The architecture in context
The system we are building
The state-space experiment represents memory through diagonal complex recurrences. Rather than choosing every pole randomly or learning all poles by gradient descent, one arm estimates them from delayed output snapshots using dynamic mode decomposition. The resulting recurrent features feed a ridge readout trained against the observed outputs.
Who does what in the stack
- PyTorch complex tensors
- Implement diagonal recurrent states.
- SVD / eigvals
- Identify a reduced dynamics operator.
- Ridge solve
- Fits output channels on a fixed recurrent representation.
The project connects output-based pole estimation to the same feature/readout interface used by generic and gradient-trained controls. Complex arithmetic represents oscillation and decay together. The fitted poles are inferred from training observations, not supplied as an oracle in this arm.
From module map to executable structure
Inside Complex diagonal state-space model
Default 64 complex states, input 1/output 2, train length 96/test 192;1,500 training and 400 test sequences.
| Layer or branch | Output shape | Implementation detail |
|---|---|---|
| Input projection | B × N complex | Fixed projection P maps d input channels to N modes. |
| Diagonal recurrence | B × N complex | h[t]=p*h[t−1]+u[t]P; one multiplier per mode. |
| Real feature vector | B × 2N | Concatenate real and imaginary state components. |
| Ridge readout | B × dout | Fit linear output coefficients on training states. |
| Alternative pole estimators | N complex poles | Fixed pole families, learned poles or Hankel/DMD inference are separate arms. |
The diagonal transition replaces a dense N×N state multiply with N elementwise products. Complex conjugate structure or real/imaginary expansion represents damped oscillation. A state model can be cheap to advance and still be poorly identified; pole radius, input projection and observability determine whether the output is useful.
The equation and the update
Default learned-pole arm 600 Adam steps at 5e-3; three seeds. DMD and fixed-pole arms do not use that optimizer budget. Known-system poles in a matched arm supply privileged information.
Implementation card / no invented benchmarks
Capacity, budget and execution evidence
- Parameters / retained state
- State is N complex values per sequence; readout width 2N before any bias convention. Pole-learning parameters are distinct from fixed projection and readout storage.
- Duration and hardware evidence
- No generic speedup is inferred from diagonal structure; training, identification and rollout need separate timings.
- Source coordinates
- E40 lines 65–155 and configuration 159–166; E41 metrics
- Current reproduction context
- Current workstation, supplied by the author: Apple M4, 128 GB unified RAM, 40 GPU cores and 16 CPU cores. This is context for prospective reproduction, not attribution of every archived run. Python and framework versions are not fully locked for these historical sources; declarations, when available, are identified separately.
Counts above are calculated from the stated layer shapes unless identified as saved measurements. They exclude optimizer state and nontrainable buffers. No archived training was rerun for this revision.
What these design choices change
Increasing N adds dynamical modes but may make the readout ill-conditioned. Longer Hankel windows change the identification problem, not just batch size. The source projects inferred eigenvalues inside radius .999; this avoids explosive autonomous modes but can bias nearly undamped dynamics.
Reproduction and measurement protocol
Check the impulse response of a single real pole and a conjugate pair against the recurrence. Reset state at sequence boundaries. Keep noisy observation targets separate from the clean reference used for evaluation, and report matched-pole arms as oracle-assisted rather than fully learned.
For a new run, save the resolved Python/framework versions, backend, dtype, seed, input shapes, batch size and exact source revision. Start with one batch and one update. Log training steps separately from epochs or environment steps. Do not equate the configured maximum with a completed budget or convergence.
Measure initialization/compilation, data preparation, warmed forward pass, training updates and evaluation separately. Synchronize accelerator work around timed regions using the chosen framework’s supported mechanism. Report peak process memory and framework allocation separately; parameter bytes exclude activations, gradients, optimizer state and input buffers. On a shared machine, begin with a single CPU worker and a small batch rather than claiming all available resources.
A closer look at the implementation
The code that carries the idea
The snippet forms paired Hankel matrices, truncates the SVD and computes a reduced transition matrix. Its eigenvalues become discrete-time poles. Radial projection caps magnitudes at 0.999 while preserving phase, preventing an estimated unstable pole from exploding during longer rollouts.
def infer_poles_dmd(y, N, delay=None):
"""Closed-form system-ID of the dynamics operator from observed output (Hankel delay-embedding + DMD). NO oracle, NO
gradient. Returns up to N discrete poles (complex). This is the 'model the substrate' / OSNR operator-inference arm."""
# average a few sequences' Hankel into one identification (use the first training seqs)
yy = y[:64].reshape(y.shape[0] if y.shape[0] < 64 else 64, y.shape[1], -1).mean(-1) # (b,T) collapse channels
b, T = yy.shape
L = delay or max(2 * N, 16)
L = min(L, T - 2)
# build Hankel snapshot pairs across all sequences
X0, X1 = [], []
for i in range(b):
H = torch.stack([yy[i, k:k + L] for k in range(T - L)], 1) # (L, T-L)
X0.append(H[:, :-1]); X1.append(H[:, 1:])
X0 = torch.cat(X0, 1); X1 = torch.cat(X1, 1) # (L, m)
U, S, Vh = torch.linalg.svd(X0, full_matrices=False)
r = min(N, (S > 1e-8 * S[0]).sum().item(), L)
Ur, Sr, Vr = U[:, :r], S[:r], Vh[:r, :].conj().T
Atil = Ur.conj().T @ X1 @ Vr @ torch.diag(1.0 / Sr) # reduced DMD operator
ev = torch.linalg.eigvals(Atil) # discrete poles
mag = ev.abs() # ENFORCE STABILITY (physical prior): |pole| <= 0.999.
ev = ev * (torch.clamp(mag, max=0.999) / (mag + 1e-12)) # project unstable modes onto the unit disk, keep phaseVerbatim archive excerpt from closed_form_neat_osnr_ssm.py. Context-dependent historical code, not a standalone runnable program. Comments retain their original wording; the article distinguishes implemented behavior from stale or overbroad comments.
The boundary that matters
Projection changes the estimate; it is a modeling prior, not proof that the identified system is correct. Collapsing output channels before identification can discard modes, and repeated-pole padding needs a shape test when too few eigenvalues are retained. Stability alone does not guarantee accurate extrapolation.
