The Same Thing Twice

Simulate first, agree afterward

A cross-domain parallel is cheap to assert and easy to dress up. The only version worth keeping is one that survives a measurement it could have failed. So every claim below gets at least two routes to the same number — one analytic, one simulated — written so that neither can see the other's answer. Where they disagree, the disagreement is the finding.

The seven investigations are below. Most produced a surprise, and in every case the surprise was a measurement error with a diagnosable cause rather than a broken theorem. One of those errors recurred four times in three disguises before I understood it; that is collected near the end, because it turned out to be the most transferable part of the whole exercise.

Rigidity is a counting problem before it is a geometry problem

CAD constraint solvers typically throw the whole system at Newton–Raphson and discover over-constraint when the Jacobian goes singular. Combinatorial rigidity theory answered this in 1970, in a literature that barely overlaps.

Laman's theorem characterises generic rigidity in the plane by edge counting alone, and the pebble game decides it in near-linear time on the constraint graph — before any numerics run, and with the ability to name which existing constraints make a new one redundant.

That is the known result that diagonal-bracing every cell of a grid is wasteful, reproduced from the graph with nothing but the counting rule. A hinge — two rigid triangles sharing one joint — comes back as two overlapping rigid components, which is what a hinge is.

I audited this tool, and it failed

Laman's theorem describes generic configurations, and CAD sketches are the opposite of generic — axis-aligned, mirror-symmetric, laid out on round numbers. I wrote that down as a caveat and moved on, which is not the same as measuring it. So: 18 CAD-like sketches, each paired with a generic twin carrying the identical constraint graph at random coordinates. The pebble game sees only the graph, so its verdict is the same for both; any divergence is caused by geometry alone.

Each was checked against the rank of the true constraint Jacobian at its actual coordinates. The audit runs the shipped engine, not a reimplementation of it.

Regular coordinates are not the problem; degeneracy is. A 4×4 braced lattice at exact integer coordinates is answered correctly. What breaks the method is mirror symmetry, collinearity and concentricity — and a collinear chain is reported well-constrained while actually carrying three internal degrees of freedom. Perturb any of those coordinates by 10⁻⁶ and the combinatorial answer becomes right again, which is the signature of a genuine genericity failure.

And one failure was not genericity at all. This page used to say that horizontal-plus-vertical on a point pair reads redundant "though it is geometrically fine" — attributing it to special coordinates. That was wrong. The generic twin appeared to fail identically — but that twin used random coordinates that did not satisfy the horizontal and vertical constraints, so it was not a realisation of the sketch at all, and the argument built on it was unsound. Redone properly (points sampled on the constraint variety, including the square and a 10⁻⁴ sliver), the rectangle has true rank 8 at every sample while the pebble game reports 6. The conclusion stands; the evidence for it is now the right kind. The real cause: the (2,3) pebble game is a theorem about distance constraints, while direction constraints obey (2,2)-sparsity, so mixing them needs a mixed matroid. It also errs in the opposite direction, under-counting independence rather than over-counting it. A caveat pointing the wrong way is worse than no caveat.

The engine now returns a warnings array saying both of these in place, and the bench shows them. The method still earns its keep on ordinary braced structures — but "fast pre-filter, not a proof" was too gentle, and the tool was not being honest about which of its answers to distrust.

Which constraint should I delete?

When the solver says a sketch is over-constrained, the number of constraints to remove is not up for negotiation: every maximal independent set of a matroid has the same size, so it is always constraints minus rank. What is open is which ones. The engine answered that by accident of drawing order — it flags whatever you drew last, and reversing the order changed the flagged set in every one of 60 test sketches.

I expected the explanation it gives for each flagged bar to be too broad, since it lists every constraint inside the failed dense block. It is not. The failed search reaches the minimal tight block containing both endpoints, a tight block carries exactly 2|V|−3 edges so nothing else can sit inside it, and that block plus the rejected bar is the circuit. On 145 redundant bars the engine's list matched the true circuit, computed independently from the Jacobian, every time. A prediction that failed in the useful direction.

The fix needs no new machinery. Run the pebble game with the constraints you least want to lose inserted first; the ones it rejects are a minimum-cost repair. That is Kruskal's algorithm — sort, then keep whatever stays independent — the same procedure that built the spanning tree in investigation 05, with "acyclic" swapped for "Laman-sparse". The greedy algorithm is optimal on exactly the structures that are matroids, and both of these are.

The bench now has a Suggested repair panel in its details sheet. Change what you prefer to keep and the set of bars changes; the count never does. Apply it and the sketch goes from over-braced to rigid. It refuses to apply when direction constraints are present, for the reason in the audit above: wrong matroid, wrong answer.

One eigenvalue, four textbooks

Mixing time in Markov chains, expansion in graph theory, the first non-rigid vibration mode in structural mechanics, and iterative-solver convergence in numerical PDE. Taught as four unrelated facts; λ₂ of the Laplacian is all of them.

λ₂ was computed once per graph. Each of the four quantities was then measured by its own simulation: a lazy random walk released on a node, brute-force search over all 2ⁿ⁻¹ cuts, a Verlet integration of masses and springs read through an FFT, and Richardson iteration with optimal damping.

The true sparsest cut lands inside the Cheeger window every time, pressed against the lower bound exactly as a bottlenecked graph should be.

Where it stops being true: every graph above is regular, and that is load-bearing. On an irregular graph the four come apart — mixing and expansion follow the normalized Laplacian, unit-mass vibration follows the combinatorial one.

The repair, and what it says the object really is

That split suggested the unifying object was never the graph but the quadratic form xᵀLx under a choice of inner product — so giving each node a mass equal to its degree, solving Lv = λDv instead of Lv = λv, should put the four back together on any graph. Tested on four deliberately lopsided graphs, degrees running from 1 to 7:

The tidiest part came from the column I wasn't watching. Plain Richardson tracks the combinatorial condition number; Jacobi-preconditioned Richardson tracks the generalized one. Changing a solver's preconditioner and changing a structure's mass matrix are the same act — both are a choice of inner product, and both move the answer from one λ₂ to the other.

A multimeter that counts spanning trees

Effective resistance in circuit theory, commute time in Markov chains, and edge frequency in uniform spanning trees. Three fields, and the prediction is specific enough that it cannot survive by coincidence.

For a graph of 1 Ω resistors, the resistance measured across an edge should equal the fraction of all spanning trees that contain that edge. Four independent routes were run: a linear solve of Ohm's law, Monte-Carlo random walks counting steps and nothing else, brute-force enumeration of every spanning tree, and the Laplacian pseudoinverse.

Edge 2–3 sits in exactly 81 of 108 spanning trees. A meter across it reads 0.75 Ω. You can answer a combinatorial enumeration question with a breadboard.

Unlike λ₂, this needed no regularity — it holds on degrees 2 through 4 with no patch. Resistance is built from the Laplacian pseudoinverse rather than its spectrum, so the degree weighting is absorbed into the 2m factor instead of having to be cancelled.

π, from a grid of identical resistors

On the infinite square lattice the resistance between diagonal neighbours should be 2/π. Nothing transcendental is present in the setup: every resistor is 1 Ω and every node is an integer lattice point.

The convergence order came out p = 1.993 — boundary corrections fall like 1/L², measured from the computed values rather than assumed.

And the adjacent case is forced by conservation

Foster's theorem says the effective resistances sum to exactly n − 1. On a lattice where every edge is equivalent by symmetry, with coordination number z and therefore nz/2 edges, every edge must carry R = 2/z. That prediction was written down before any of these were computed.

Each lattice was generated from coordinates and bonded purely by distance, with its interior coordination number asserted before anything was solved. The sharp test is the last pair: kagome and square are different lattices that share z = 4, and they land on the same number to seven decimal places.

A 1956 algorithm computes a 2005 invariant

The death times in the 0-dimensional persistence barcode, the merge heights of a single-linkage dendrogram, and the edge weights of a minimum spanning tree. Topological data analysis, statistical clustering, and Kruskal.

There is a trap in testing this. The usual way to compute H₀ persistence is union-find over sorted edges — which is Kruskal, so comparing those two proves nothing. So here H₀ is computed the long way: sort every simplex by filtration value, build the GF(2) boundary matrix, and reduce it with column operations until each column's lowest nonzero row is unique. That shares no structure with greedy edge selection or with agglomerative merging.

Four routes, on three point clouds: Prim growing one tree by nearest outside vertex; naive agglomerative clustering that never builds a graph at all; boundary-matrix reduction; and a sweep that rebuilds the graph from scratch at each radius and counts components by breadth-first search.

Not agreement to machine precision — the same floating-point numbers, on all three clouds. Four computations that share no code produce one list.

The practical reading is in the clustered case. Three Gaussian blobs give bars of 0.810, 0.864, 0.866, 1.040 and 2.221: a noise floor, then two that stand clear of it. Two long finite bars plus the one infinite bar is the barcode saying three clusters — and it is the identical gap a statistician looks for when deciding where to cut a dendrogram. Choosing the number of clusters and reading a barcode are one act, performed by people who cite different papers for it.

Backpropagation was already a hundred years old

The adjoint method in computational physics, backpropagation in machine learning, and Pontryagin's costate equation in optimal control are three derivations of reverse-mode gradients. The literatures barely overlap. The question is not only whether they agree, but which two are the same algorithm and which is merely the continuum limit of the others.

A damped pendulum under a time-varying torque, integrated by explicit Euler, with a loss over the trajectory. The gradient of that loss with respect to the damping coefficient was then computed four ways, sharing no code: central finite differences as ground truth; a tape-based reverse-mode autodiff engine written from scratch and run over the unrolled simulation; a hand-derived backward recursion for the Euler scheme; and the costate ODE integrated backward from its terminal condition.

Backpropagation and the discrete adjoint do not merely agree — they agree to 1.6×10⁻¹⁵, and across every one of the two thousand control gradients the largest disagreement is 1.4×10⁻¹⁷. That is not two methods landing on one answer. That is one algorithm, written down twice by people who did not know about each other.

The continuous adjoint is the interesting one, because it is wrong at any finite step size — and correctly so. It differentiates the differential equation; the other two differentiate the program that approximates it. The gap is first order in the step:

Halving the step halves the error, converging on a ratio of exactly 2. This is the distinction between discretize-then-optimize and optimize-then-discretize, and it is not academic: an optimizer fed the continuous gradient is descending a surface slightly different from the one it is being scored on.

Worth recording: this is the only investigation here whose prediction survived on the first run, with no anomaly to chase. Written down in advance, confirmed, nothing to debug. Reported because the other three were not like that, and a notebook that only shows the clean runs is lying by selection.

One theorem, four names, and one that would not convert

Max-flow equals min-cut. The value of a zero-sum game equals the value of its dual game. A linear program's optimum equals its dual's. Either a system has a solution or Farkas hands you a certificate that it has none. Investigation 02 met the first of these; this asks whether the other three are the same object, not just the same slogan.

These are textbook theorems, so agreement of values proves little. The tests worth running are the ones where a certificate in one language has to be translated into another and still check out. All on 60 random directed networks (5–7 nodes, 5–19 edges, integer capacities, maximum flow 1–18), each solved four independent ways: Edmonds–Karp augmenting paths, brute-force enumeration of every s–t cut, the flow LP, and the cut LP written out by hand.

The Farkas certificate is a cut, literally. For the system “send d units from s to t within capacities,” I did not derive the certificate by hand. I assembled the constraint matrix, asked a generic solver for a vector y with yA ≥ 0 and y·b < 0, and read what came back. Normalising the source multiplier to 1, the conditions on y are exactly the cut LP’s constraints (a potential on each node, an edge variable at least the potential drop across it), and y·b < 0 says that cut costs less than d. Thresholding the potentials produced a genuine cut of capacity below d in all 60 cases. When the demand was exactly the max flow, the solver found no certificate in all 60 — the “exactly one of the two” in Farkas, observed.

A prediction that could have failed. The usual proof that the cut LP has integral optima averages over thresholds: the LP value equals the mean capacity of the level sets, and no cut beats the minimum, so if the potential is optimal then every level set must itself be a minimum cut. That is checkable on a deliberately fractional optimum. On the 15 networks with two or more distinct minimum cuts I took the 50/50 mix of two of them (a potential taking the value ½ on the nodes they disagree about), confirmed it is dual-optimal, and tested the level sets at thresholds 0.25, 0.5 and 0.75. All 15 were minimum cuts. Fractional dual optima are not noise to be rounded away; they are convex combinations of cuts that all cut the same.

The game, and the version that failed

The standard way to see any LP as a game is Dantzig’s embedding into a symmetric game, with the LP optimum recovered as x/t from the equilibrium strategy. I ran it, solving the game by multiplicative-weights self-play so that nothing shared code with the LP solver. It did not work. Over 12 networks the recovered flow value was off by 30% at the median (range 1% to 60%), and the recovered flow broke a capacity or conservation constraint by 0.79 units at the median. The cause is visible in the output: the normaliser t, the weight on the last row, settled near 0.04, so dividing by it magnified an ordinary 1/√T equilibrium error by 25×. The embedding is correct and its conditioning is dreadful, and a reduction that is exact on paper can be unusable as an algorithm.

A different game works. Let one player choose an s–t path and the other choose an edge, with payoff 1/c_e if the edge lies on the path. The path player minimises. Its value should be 1/(max flow). Then the path player’s optimal mixed strategy, scaled by the max flow, should be a maximum flow written as a path decomposition, and the edge player’s strategy should pile onto the minimum cut. With the game value from an LP: |value − 1/F| ≤ 6×10⁻¹⁷ on all 60. With 30,000 rounds of self-play, the bracket around the value contained 1/F in 60 of 60 (median width 3.5% of the value, worst 22%), the scaled path strategy loaded edges to 1.004× capacity at the median (worst 1.11×), and the edge player put a median 97% of its weight on edges belonging to a minimum cut (worst 66%).

So the picture is: the path player plays the flow, the edge player plays the cut, and the value is their shared price. One saddle point, three vocabularies. The earlier failure was about conditioning, not correctness: the structured game recovered both the flow and the cut from the same kind of self-play.

At scale. At a few thousand nodes the strategy sets cannot be enumerated, so the path player became a best response (one shortest-path computation per round, with edge lengths q_e/c_e) against multiplicative-weights edge play. Two kinds of network, each checked against an exact max flow. On plain random sparse digraphs (n up to 3,000) the result was uninformative: the minimum cut was just the star around the source or sink, which any method finds. So I planted a bottleneck: two dense random halves, 25 fat edges out of s and into t (star capacity 250 each), joined by six edges of capacity 3, so the true cut is 18 and sits in the middle of the graph. On n = 200, 1,000 and 3,000 (up to 21,000 edges, 4,000–6,000 rounds) the sweep over level sets of the distance potential from the averaged edge strategy returned a cut of capacity exactly 18 in all three, and in all three it was the planted split, n/2 nodes on each side. The averaged path strategy scaled by the max flow loaded edges to 1.000× capacity. The surprise is what did not work: the plain lower bound on the game value, the best response to the averaged edge strategy, stayed loose (0.27–0.53 of the true value at the end, and worse as n grew) and only 27–53% of the edge weight sat on the crossing edges. The raw strategy is a poor certificate of the value, and its level sets are an excellent certificate of the cut. That is the Investigation 06 level-set theorem doing real work in a regime where nothing else did. Three instances is a demonstration, not a rate: I did not measure a 1/√T law, the planted network is symmetric (six equal edges make the flow side unusually easy), and the failed Dantzig embedding was not revisited.

What this does and does not show. Nothing here is new mathematics; every identity is a theorem and was used as a target. What was measured is that the translations hold on concrete instances — certificate to cut, fractional optimum to level sets, equilibrium to flow and cut — and that the one translation I tried from the textbook’s reduction was numerically unusable. Networks were small enough to enumerate every cut, and the self-play bounds are loose; this says nothing about large instances. Code and the per-instance results are in the project.

Buying memory with arithmetic, at the exact exchange rate

Investigation 04 ended with a cost: reverse-mode differentiation through a chain of l steps must remember the forward trajectory. Store all l states and the backward pass is free; store almost nothing and you recompute constantly. The question is what the best possible exchange rate is, and whether a known answer survives being checked three independent ways.

The model. Keep the initial state plus s further checkpoint slots, and one working state. The backward pass needs states l−1, l−2, …, 0 in that order. Moving forward one step costs 1; loading a stored state, discarding one, and using the working state in a backward step are free. Splitting at the best interior checkpoint k gives a recurrence, T(l,s) = min_k [ k + T(l−k, s−1) + T(k, s) ], with T(1,s) = 0 and T(l,0) = l(l−1)/2 (no spare slot means recomputing from the start every time).

A convention error, caught by the test. I first checked the recurrence against an exhaustive search of the equivalent pebble game on a line and they disagreed in 35 of 45 cases. The recurrence was fine; my search charged the in-flight state and the initial state differently. Once both used the same bookkeeping, the search and the recurrence agreed on every one of 45 cases (l = 2–10, s = 0–4). This is the same species of error as the rest of this page — two honest calculations that disagree because they are counting different things — and, as before, it showed up as a clean, plausible, wrong table.

The closed form. The binomial formula from the literature (Griewank–Walther, the “revolve” schedule) fits exactly, but only with the right offset: with r the least integer such that C(s+1+r, r) ≥ l,

Running it for real. Counting is not the same as differentiating, so I executed the optimal schedule on an actual nonlinear chain, x_i+1 = tanh(W_ix_i + b_i) in 8 dimensions, and compared the gradient to the store-everything version. The gradients agreed to the last bit (max error 0.0) in all five runs; the number of extra forward steps equalled the recurrence’s prediction exactly (15, 69, 172, 374, 255 for (l,s) = (10,2), (30,3), (64,4), (100,3), (100,6)), and peak storage was s+2 states each time. (My first version of that run came in below the prediction, because its base case quietly stored a whole segment; the counter was right and my memory accounting was not. Fixed to recompute, then it matched.)

Reading the curve. The cost is polynomial in l for fixed memory (roughly l^1+1/s) but becomes only logarithmic overhead once the slot count grows like log l — a few dozen states make a million-step chain affordable at a handful of recomputed passes. The uniform √l scheme looks better in the table, but it is not the same machine: it stores a whole segment while replaying it, so it spends about 2√l states to get a single extra pass. The binomial schedule is the optimum for the model above, not for every model.

What this does and does not show. The schedule and formula are known results; what was done is verify them from scratch three ways and see the gradients come out right. Not tested: chains whose steps differ in cost or state size, branching graphs, or disk-like two-level memory. The link to reversible computing (Bennett’s pebbling) is the same game with the extra constraint that nothing may be erased for free, and was not examined.

The bug that kept coming back

One error appeared four times across these investigations, in three different costumes, before I understood what it was. It never looked like an error. It looked like a clean, confident number that happened to be wrong.

Any measurement that reads a mode off a simulation can fail to excite that mode, and will then report the next one up. There is no warning, no ragged output, no failed assertion — just a plausible value, usually matching a neighbouring eigenvalue to five or six digits.

• A probe on a nodal point. A random 3-regular graph reported its first vibration mode 59% high. Node 0 sat almost exactly on a node of the Fiedler mode — amplitude −0.018 — so the sensor was blind to the mode being measured and returned √λ₄ instead. A single accelerometer in the wrong place does exactly this on real hardware.
• The same thing in the JavaScript port, because the port watched one fixed random projection rather than every node. The identical symptom, reintroduced while translating a fix I had already made.
• A random walk starting on a near-nodal vertex (overlap 0.0039 with the Fiedler mode). It decayed at the λ₃ rate and matched 1 − λ₃/2 to five digits.
• A random initial condition nearly orthogonal to a mode — amplitude 0.0067 where typical is about 1. Pure chance, and it reappears whenever the random stream shifts, which makes it look like flakiness rather than a systematic flaw.

The durable fix is not a better peak detector. It is several independent excitations: sum the spectra over three random initial conditions, take the slowest decay over four random walk starts. A mode one trial misses, another catches. Single-shot excitation was the bug the whole time — probing at a single point was only its first costume.

Two smaller mistakes worth keeping

• The profiled guess. Rigidity analysis took 578 ms at 400 joints. The DFS inner loop looked like the culprit, so it got optimised; the time did not move. The real cost was a Set rebuilt inside a filter predicate, O(m) times per redundancy. Hoisting it gave 51 ms. Optimising the thing that looks slow is not profiling.
• The assumed convergence order. Richardson extrapolation on the resistor grid made the answer three orders of magnitude worse, because it assumed error ∝ 1/L when the data said 1/L². Measuring p from three grid sizes instead of assuming it took π from five digits to eight.

Threads not yet pulled

• Mixed matroid, done. Sketches mixing bars with horizontal/vertical constraints are now analysed by the rank of the constraint Jacobian at a generic point on the constraint variety (one random x per vertical-connected class, one random y per horizontal-connected class), rather than by pebble counting. It is still a matroid, so the same greedy repair applies unchanged — only the independence oracle differs. Checked against a separate Python implementation on 120 random mixed sketches at a different sample: statuses, ranks, internal DOF and repair sets all agree (0 mismatches), and the old pebble path misreported the rank on 80 of them. Caveat: the Python check shares the same formula for surviving trivial motions, so internal DOF is confirmed numerically, not derived independently. Parallel, perpendicular, angle and non-point objects remain unhandled.
• Duality at scale, the rest. Investigation 06 now has a thousands-of-nodes check on planted bottlenecks. Still open: whether the self-play error follows a 1/√T law (not measured), asymmetric bottlenecks where the flow side is not trivially easy, and whether a better normalisation rescues the Dantzig embedding.
• Checkpointing beyond a chain. Investigation 07 solved the uniform-cost line. Non-uniform step costs, branching graphs (the general pebble game), and the reversible-computing variant where erasure is not free are untested.