System Identification

Suppose you hand a robotics team thirty minutes of logged data from a prototype drone: commanded motor thrusts, accelerometer readings, gyro readings, timestamps. Nobody handed them a physics engine. Yet they need something that answers questions like "if we pulse the left rotor now, how will the roll angle respond over the next half second?" System identification is the discipline of building that something: starting from recorded trajectories of inputs and outputs, and recovering a predictive dynamical model of the system that produced them. The team's logs are passive records of what already happened; identification turns them into an active instrument, a simulator they can query about futures that have not occurred yet. That transformation, from a fixed record to an answerable counterfactual, is the entire point of the enterprise, and every technique in this chapter is a different route to it.

This sits in a specific place relative to the last chapter. In Bayesian Filtering and Belief States, we took the dynamics $f$ and observation function $g$ as given and inferred the hidden state as data streamed in. System identification flips the problem: the dynamics themselves are unknown, and we infer them offline from batches of data. You can think of filtering as running the world model forward with uncertainty, while identification is building the world model in the first place. In fact the two problems are duals in a precise sense: filtering assumes the model and estimates the trajectory state; identification assumes (noisily measured) trajectory evidence and estimates the model. Classical estimation theory couples them tightly: for a linear state-space model with Gaussian noises, the Kalman predictor supplies the conditional predictions evaluated inside the prediction-error method. Every model-based planner, controller, and agent discussed later in this book is downstream of this step. If Chapter 1 gave us MDPs and POMDPs as the grammar of sequential decision making, this chapter is about estimating the transition function that makes those formalisms usable. A planner written against a wrong transition function is not slow or suboptimal; it is confidently wrong, and this chapter exists largely to explain why that happens and how often to expect it.

We will work through four connected ideas. First, identifiability: which models can the data even distinguish, and what makes a dataset informative at all (excitation). Second, parametric identification: fitting a model with a fixed functional form and a finite parameter vector, culminating in the prediction-error method. Third, nonparametric and subspace methods: estimating dynamics directly as impulse responses or by realizing state-space matrices from data geometry. Fourth, neural system identification: replacing hand-picked function classes with learned networks, which connects classical practice directly to the modern world models of Part V and Part VI. These four ideas form roughly a historical arc as well as a conceptual one: identifiability theory and classical parametric methods matured through the 1960s and 1970s; nonparametric estimation and realization methods developed across the following decades; canonical subspace-identification algorithms emerged in the 1990s; and neural approaches have dominated the last decade. But the arc is not a story of progress leaving older ideas behind. As you will see, each generation inherits the same structural constraints, excitation requirements, and evaluation pitfalls, only wearing new clothes.

Notice that this definition quietly embeds a judgment call: "well enough for the intended use." There is no single universal criterion of model quality. A model good enough for simulating three steps ahead may be useless for planning twenty-step rollouts; a model matched to one operating regime may be dangerous in another. Throughout the chapter we will distinguish predictive fidelity (matching future outputs), representation quality (whether latent variables carry meaning), uncertainty quality (whether confidence intervals mean what they say), and decision usefulness (whether planners built on the model succeed). Collapsing these into one notion of "fit" is among the most common failure modes in applied work.

Before any estimation, we fix vocabulary. Throughout this chapter, a discrete-time system generates trajectories. Depending on how lucky we are with sensors, we see one of two things:

• Full-state data: at every time step $t$ we observe the state $x_t \in \mathbb{R}^n$ together with the applied action or input $u_t \in \mathbb{R}^m$. Many simulated settings and instrumented physical setups provide this. Full-state data is rare in the wild but pedagogically precious: with it, identifying dynamics becomes ordinary regression, and all the statistical machinery you already know applies directly.
• Input-output data: we observe only measurements $y_t \in \mathbb{R}^p$ produced by an observation map, $y_t = g_\star(x_t, u_t) + v_t$, never the state itself. This is the standard POMDP-flavored situation from Part II, Chapter 1. Nearly everything interesting in practice lives here, because real sensors measure projections of reality, not reality. Your drone's accelerometer measures acceleration, not position and velocity; its gyro measures angular rates, not attitude.

The underlying generative story is a stochastic discrete-time system:

where $f_\star$ and $g_\star$ are the unknown true maps, $w_t$ is process noise (disturbances inside the dynamics), and $v_t$ is measurement noise (sensor corruption). Including $u_t$ in $g_\star$ permits direct feedthrough, such as an accelerometer responding immediately to a motor command; examples without feedthrough simply ignore that argument. This distinction between process noise and measurement noise matters downstream. Process noise means the world does something unpredictable at each step; even the true model cannot forecast those disturbances, so there is an irreducible floor on rollout accuracy no estimator can beat. Measurement noise means our view is blurred but the underlying evolution may be perfectly predictable; a clever estimator can average the blur away. Confusing the two leads to characteristic mistakes: treating sensor jitter as if physics were random wastes modeling effort and inflates uncertainty where none exists, while treating turbulence as sensor error produces overconfident models that fail exactly when conditions get rough.

We will state our standing assumptions explicitly, because every identification guarantee rests on them:

• Markov assumption: the next state depends on the past only through the current state, so $(x_t, u_t)$ suffices to predict $x_{t+1}$. When the assumption fails, because the system has hidden memory, delays, or slowly varying unmodeled parameters, the standard remedy is to augment the state until Markovianity approximately holds, by stacking lagged observations into the effective state. Most autoregressive models are exactly this trick wearing simple clothes.
• Stationarity: the maps $f_\star$, $g_\star$ and noise statistics do not change over the recording window. A drone whose battery drains or whose payload shifts violates this. Battery drain changes the gain from command to thrust; payload shift changes the inertia. Both turn what looks like one long dataset into a sequence of slightly different systems, and any single fitted model will be a compromise nobody asked for.
• Noise conventions: in the clean theory and simulations below, $w_t$ and $v_t$ are zero mean and temporally uncorrelated unless a section says otherwise; Gaussianity is optional except where maximum-likelihood or efficiency statements invoke it explicitly. For open-loop least squares, current innovations must also be independent of, or at least orthogonal to, past regressors and actions. The worked drone example later uses process noise without measurement noise; the ARX discussion then shows why adding sensor noise changes the plain-OLS error structure.
• Causal exogeneity (for many results): the regressors must be uncorrelated with the current disturbance. Closed-loop data can break this condition when feedback carries serially correlated or contemporaneous noise into the input, entangling signal and error. Feedback alone is not enough to imply bias: with white innovations and regressors built only from past variables, the current disturbance can remain orthogonal to the regressors. We return to this distinction in the excitation section.
• Correct-enough model class: the estimator can represent $f_\star$, or at least the projection of $f_\star$ that matters for our use. The qualifier "projection" matters. A model that omits high-frequency flexibility nobody queries might still be perfect for its consumer; a model missing a mode the planner will hit is broken regardless of its training loss.

Two distinctions from Part I, Chapter 3 deserve repeating here because they organize everything that follows. First, a model can score excellently on one-step-ahead prediction and still be useless for closed-loop decision making, where the model's own outputs are fed back recursively. The reason is compounding: at step one the model sees ground truth, but at step two it sees its own imperfect output, at step three a doubly corrupted estimate, and small per-step errors accumulate geometrically in systems with unstable or lightly damped modes. We will measure this gap numerically later in the chapter, and it will recur in essentially every evaluation we perform. Second, identifiability is not the same as accuracy: an unidentified model can accidentally predict well on the data distribution that generated it while answering counterfactual questions (what happens if we act differently?) catastrophically badly. Planners ask exactly those counterfactual questions, so we care about identification depth, not surface fit. A planning algorithm is, structurally, a machine for proposing unusual action sequences; the unusual ones are precisely where a shallow-fit model has had no chance to be tested, and precisely where it is consulted most earnestly.

This section covers the question that precedes every fitting procedure: does the dataset contain enough information to pin down the model at all? It is tempting to skip straight to algorithms, since fitting code runs happily on garbage data without complaint. But the concepts here determine whether that run meant anything, and they explain phenomena that otherwise look mysterious, like estimators returning wildly different parameters on different days despite nominally similar datasets.

Identifiability is a property of the triple (system, model class, dataset). Keep all three components in mind, because failures of each kind have different cures. Structural problems belong to the class and cannot be fixed by data collection; practical problems belong to the dataset and usually can be. The model class defines a mapping from parameters $\theta$ to input-output behavior. The data identifies $\theta$ if no other parameter value in the class produces the same distribution of trajectories on experiments like the ones we ran.

The distinction between the two definitions is worth internalizing via their respective symptoms. Structural non-identifiability manifests as flat directions in the log-likelihood that stay flat no matter how much data arrives: infinitely many optima, exactly tied. Practical non-identifiability manifests as directions that are technically distinguishable but only barely. Their variance still shrinks with sample count when Fisher information is positive, but the prefactor can be so large that improvement is negligible at realistic dataset sizes. The Fisher-information framing makes the connection quantitative: practical identifiability asks whether the log-likelihood's curvature matrix has full rank and reasonable condition number near the optimum.

Structural non-identifiability often hides in symmetry. Take any linear state-space realization $(A, B, C)$. For any invertible $T$, the transformed pair $(TAT^{-1}, TB)$ produces identical input-output behavior when combined with $CT^{-1}$. The verification is a two-line computation: $(CT^{-1})(TAT^{-1})^k(TB) = CA^kB$, so every impulse-response coefficient, hence every input-output trajectory, survives the change of coordinates untouched. Input-output data therefore can never recover $A$ itself, only $A$ up to a similarity transformation. This is not a defect of your optimizer; it is a fact about information. No quantity computable from input-output records distinguishes the original from the transformed system, in the same way that a photograph of a sculpture cannot reveal which direction you call north. Two attitudes are available toward this ambiguity, and both appear in practice. Subspace methods (later in this chapter) embrace it and return whatever basis they like, while planners happily accept it because plans depend on predicted measurements, not on the internal coordinate system. Only when latents must correspond to interpretable physical quantities, as in scientific modeling, does the ambiguity become a substantive obstacle requiring extra information to break.

Practical non-identifiability is dataset-dependent and far more common in practice. Consider a linearized aircraft model and a recording made entirely during cruise at one airspeed and altitude. Hundreds of thousands of samples, yet the aerodynamic derivatives governing stall behavior are untouched by that data: different stall-coefficient values produce indistinguishable cruise recordings. The sample count was never the binding constraint. The binding constraint was behavioral coverage. If anything, the abundance of cruise samples makes the problem worse psychologically, because a large dataset creates a false sense of security; the numbers were many, but the questions they answer were few and repetitive. This is the first appearance of a theme that will dominate the end of this chapter: what a dataset certifies is determined by what it varies, not by how much of it there is.

Excitation quantifies behavioral coverage for the input channel. When the input sequence probes enough independent directions, the empirical autocovariance of appropriately stacked inputs becomes positive definite, so the associated regression has a unique estimate and, under suitable limit conditions, a unique probability limit. Persistent excitation alone does not make that limit equal the truth. Consistency additionally requires a correctly specified model, exogeneity or regressor-orthogonal disturbances, and the usual stability and regularity conditions.

The "consistently, forever" clause deserves emphasis because it is where naive experiments fail. A startup chirp or initial burst of randomness gives you full-rank covariance once, at the beginning, but as the record grows those early samples dilute toward zero weight in the empirical average. Persistent excitation requires that richness continue appearing throughout the record, so that no matter how long you wait, every parameter direction keeps getting fresh evidence. Methods exist for tracking time-varying systems under much weaker conditions (forgetting factors effectively downweight old data), but for stationary offline identification on infinite horizons, persistence is the right requirement.

What does sufficient richness look like concretely? A few canonical examples calibrate intuition, ordered from least to most informative:

• Zero input identifies nothing about the input channel. Recording a passive system tells you about its autonomous dynamics and nothing about how it responds to commands. You recover eigenvalues but no $B$ matrix; the system hums along its own modes and never reveals its opinions about your instructions.
• A constant input is minimally informative. You learn the response at a single operating point; the local sensitivity of behavior to input level, which is what $B$ measures, stays invisible. The system settles into one equilibrium, and equilibrium data pins down ratios and steady-state gains while leaving the dynamic sensitivities undetermined.
• A single sinusoid is rich of a very low order. The classical result is that a sum of $m$ distinct sinusoids is persistently exciting of order roughly $2m$, so you need several frequencies, not several cycles, to probe higher-order dynamics. Repeating the same frequency for hours adds confidence at that frequency but no new directions in the regressor space.
• Filtered random noise is rich of effectively arbitrary order, which is why dithering with small random perturbations is the default cure for degenerate excitation in industrial practice. Random signals populate all directions statistically, so magnitude becomes the only design variable: enough dither to inform, not so much as to disrupt operations or violate actuator limits.

Two practical extensions matter, which we return to below.

To see the excitation hierarchy concretely, we can compute the empirical second-moment matrix of stacked inputs $\bar u_t$ for each canonical input type and record how many well-conditioned directions it actually covers.

The two extensions are these. First, for nonlinear model classes, persistent excitation generalizes to the requirement that the empirical covariance of relevant feature embeddings, e.g. $[\phi(x_t), u_t]$ stacking for the regressor you intend to fit, is nonsingular; excitation is always relative to a function class. This point generalizes all the way up the expressivity ladder: exciting for a linear model is not automatically exciting for a sine-of-position term, because feature embeddings define different notion of "directions" that need covering. The deeper the class, the harder thorough coverage becomes, and eventually coverage degrades from a checkable property to an act of faith, which is precisely why extrapolation risk dominates in the neural setting later.

And for closed-loop data, excitation becomes subtle: if the controller is a deterministic feedback law $u_t = \pi(y_t)$, the input may inherit most of its variation from disturbances, and some parameter directions can become entangled with the assumed noise model. Serially correlated or contemporaneous disturbances can then correlate the regressors with the current error and invalidate clean least-squares guarantees. With white innovations and regressors formed only from past variables, however, closed-loop operation does not by itself force that bias; lack of independent excitation can still make parameters difficult or impossible to distinguish. A common remedy is to add an externally generated dither signal, restoring excitation independently of the disturbances. This is also the provenance trap behind naive behavioral cloning from logged data: every record was generated by some policy, possibly a human, so its coverage silently shapes what either a policy learner or a learned model can certify.

Since excitation is a property of the data, we sometimes get to choose it. Experiment design asks: subject to constraints (the robot cannot smash itself, the patient cannot receive arbitrary drug doses, the flight envelope is bounded), which input profiles maximize information about the parameters we care about? Common criteria minimize $\operatorname{tr}(\mathcal I^{-1})$ (A-optimality), maximize $\log\det \mathcal I$ (D-optimality), or maximize $\lambda_{\min}(\mathcal I)$ (E-optimality), where $\mathcal I$ is the Fisher information matrix. Maximizing $\operatorname{tr}(\mathcal I)$ is another possible scalar objective, but it is not A-optimality in general because it can reward directions that are already well excited. Each choice expresses a different preference about which parameter combinations deserve precision. Guidelines worth carrying forward:

• Vary amplitude, not just timing: commands spanning the meaningful dynamic range reveal sensitivity everywhere you plan to operate. Small-amplitude excitation linearizes your knowledge around mild behavior; the model then knows nothing about whether gains change at aggressive settings.
• Cover the bandwidth that matters: fast modes need fast excitation, slow integrator-like modes need long records, and identification for aggressive control needs excitation where the control bandwidth lives. A signal rich at low frequencies leaves fast transducers of behavior completely dark, and vice versa.
• Prefer multiple short, diverse experiments to one long monotonous one when drift threatens stationarity. Diversity across experiments buys behavioral coverage while limiting each experiment's exposure to slowly accumulating drift.
• Cheap multimodal information beats expensive unimodal volume: sometimes adding a sensor helps identifiability more than doubling the record length. We connect this to multimodal sensing in Part III, Chapter 1. Intuitively, a new sensor changes the observation map so that previously hidden state components leave fingerprints; more data through the same narrow keyhole just sharpens the same partial picture.

Parametric methods commit to a function class with finitely many parameters and optimize those parameters against data. The commitment is a bargain: you accept misspecification risk in exchange for statistical efficiency, interpretability, and computational tractability. When the class contains (or nearly contains) the truth, parametric methods extract an extraordinary amount from modest data; when it does not, they converge efficiently to the wrong answer. This section builds the two workhorse approaches: regression for linear systems and the prediction-error method as the general template.

Start with the friendliest case, friendly enough that we can see the entire estimation logic laid bare before complications arrive. Suppose we have full-state data $\{x_t, u_t\}_{t=0}^{T-1}$ and posit the linear model class:

where:

• $A \in \mathbb{R}^{n \times n}$: the unknown state-transition matrix, propagating today's state forward under the system's internal dynamics
• $B \in \mathbb{R}^{n \times m}$: the unknown input matrix, converting commanded actions into state change
• $\varepsilon_t \in \mathbb{R}^n$: the process-noise residual at time $t$, the part of $x_{t+1}$ the deterministic part cannot account for

The crucial observation is that for a linear model class, each row of the unknown matrices appears in one state equation independently, so what looks like a dynamical-systems problem decomposes into ordinary multiple regression, one row at a time. Stack the unknowns row-wise into a matrix $\Theta = [A \;\; B]$ of shape $n \times (n+m)$, and define the regressor $\varphi_t = [x_t^\top, u_t^\top]^\top$. Every transition is a linear equation:

Arranging all $T$ transitions into matrices $\Phi \in \mathbb{R}^{(n+m)\times T}$ (regressors as columns) and $X_+ \in \mathbb{R}^{n \times T}$ (successor states as columns), ordinary least squares solves

which is exactly $n$ ordinary linear regressions sharing one regressor matrix, one per state dimension, each using the same $T$ transition observations. Read the middle expression physically: $\Phi \Phi^\top$ is the empirical second moment of the stacked state-input regressors used in this fit. Positive definiteness of this finite-record matrix is exactly the condition for a unique least-squares estimate. Without full rank, the displayed inverse does not exist and some parameter directions remain unidentified, although a pseudoinverse still returns one of many minimizing solutions. Persistent excitation is the related asymptotic richness condition on the input sequence: it is designed to keep the relevant empirical information matrices nonsingular as data accumulates, but it is stronger than full rank in one finite record and is not literally identical to the joint state-input rank condition here.

Under i.i.d. Gaussian innovations, least squares coincides with maximum likelihood and is asymptotically efficient under the usual regularity and excitation conditions. With lagged stochastic states as regressors, finite-sample estimates can retain an $O(1/T)$ dynamic-regression bias even though current innovations are uncorrelated with current regressors. Under non-Gaussian white noise with the same exogeneity conditions, least squares remains consistent and asymptotically unbiased but generally loses Gaussian efficiency. The clean distinction for this chapter is therefore consistency versus persistent misspecification bias, not exact finite-sample unbiasedness.

When only inputs and outputs are available, we posit an autoregressive model, the classical ARX form ("auto-regressive with exogenous input"). For a second-order system with scalar output,

where:

• $\alpha_1, \alpha_2$: coefficients on the lagged outputs; together they form the characteristic polynomial $\lambda^2 - \alpha_1 \lambda - \alpha_2$, whose roots are the poles we will shortly compare against the true system eigenvalues
• $\beta_1, \beta_2$: coefficients on the lagged inputs
• $e_t$: the one-step-ahead prediction error. Plain least squares is consistent only when this composite error is orthogonal to the lagged-output and lagged-input regressors; process noise or measurement noise entering those lagged outputs can violate that condition. Whiteness is an additional condition used by standard uncertainty and efficiency calculations.

and again least squares applies with regressor $[y_t, y_{t-1}, u_t, u_{t-1}]$. The philosophical move here is elegant: instead of estimating a hidden state and then the dynamics on it, we promote the lagged output history itself to be the state. The model implicitly treats its own past predictions as state, so ARX is a one-step predictor even when the physical system was written in state-space form. For minimal observable systems these two descriptions contain identical input-output information; they differ only in coordinates and in what happens to noise, which is exactly the caveat coming next. Fitting the ARX denominator coefficients and comparing their roots to the true system eigenvalues reveals something important: the recovery happens in input-output coordinates, a different (but equivalent, for minimal systems) description than the internal state.

Notice what changed, because it quietly invalidates the clean theory. With full-state data, the regression noise is the current process disturbance and the current regressor is built from variables available before that disturbance. With output data, eliminating the hidden state produces a composite equation error containing earlier process disturbances; sensor noise can also enter the lagged outputs. Either source can therefore appear in both today's regressor and today's equation error, biasing plain ARX least squares. In the drone example below, process noise alone is enough: the eliminated equation contains $e_t=w_{t,1}-0.75w_{t-1,1}+0.35w_{t-1,2}$, while $y_t$ contains $w_{t-1,1}$. Such correlation can persist with more data, so the estimator may converge to a shifted value. Correcting that bias leads to the general template.

The prediction-error method (PEM) generalizes linear regression into a principle for essentially any model class. Rather than regressing outputs on raw lagged quantities and hoping the implicit noise bookkeeping works out, PEM asks each candidate model to do its own optimal prediction, noise structure included. Define a one-step-ahead predictor associated with candidate model $\theta$: $\hat y_{t+1|t}(\theta)$, the conditional expectation of the next output given everything observed through time $t$, assuming the world behaves like $\theta$. The definition matters: this conditional expectation uses only data available through time $t$, so when past measurements are noisy, the optimal predictor weighs recent observations less than it would trust clean states. For a linear state-space model with Gaussian noises, computing this predictor is exactly the Kalman filter recursion from our previous chapter, run with $\theta$'s matrices. Here the two chapters meet operationally: the filter, previously the consumer of a known model, becomes an internal subroutine of the estimator searching for the model. PEM then optimizes:

typically with quadratic loss $\ell(e) = \|e\|^2$, whereupon PEM reduces to negative log-likelihood maximization for Gaussian output noise. Three properties justify why PEM became the classical centerpiece:

• Statistical grounding. With structural containment, identifiability, informative data, and standard stability and regularity assumptions, a correctly specified maximum-likelihood PEM is consistent and asymptotically efficient. Under misspecification it instead converges to the best approximation selected by its loss, which need not be the physical system.
• Noise-model honesty. Unlike bare ARX regression, PEM lets you declare whether disturbances corrupt the state dynamics (equation error), only the sensor (output error), or something in between. A wrong disturbance model can shift the pseudo-true parameters and forfeits the clean maximum-likelihood efficiency guarantee, even when some parameter estimates remain consistent under weaker conditions. When a network minimizes squared one-step error against noisy observations, it adopts a Gaussian one-step-prediction scoring rule on noisy targets; that loss alone does not specify a classical output-error model structure.
• Model-quality framing. PEM scores exactly the quantity a filter consumes: predictive likelihood of new data. Models good for filtering are judged on their own terms. This alignment between the training objective and the downstream consumer is rarer than it sounds; later we will see planning consumes models quite differently, and the mismatch will matter.

Its costs: PEM optimization is nonlinear in $\theta$ (each loss evaluation runs a filter pass), initialization-sensitive, and locally convergent at best. Each gradient or search step requires filtering the entire record, and the loss surface over multi-dimensional parameter spaces can harbor local minima corresponding to qualitatively wrong dynamics that nonetheless predict the data mediocre-but-decently. Multi-start strategies, careful initialization from cheaper methods, and structure-exploiting implementations mitigate but do not eliminate this friction. This computational friction motivated the family we visit next.

Consider a crude lateral dynamics model of a quadrotor hovering, linearized about hover and discretized at 50 Hz:

where $x_1$ is the lateral velocity proxy, $x_2$ is its rate, $u$ is the differential motor command, and $y_t = x_{t,1}$ records that only the first component of the state reaches the sensor. The complex-conjugate eigenvalues of this $A$ have modulus about $0.853$, describing a moderately damped side-to-side oscillation: errors rotate as they decay rather than disappearing in one step. We simulate it under a spectrally rich shaped-noise input, then fit (a) the full-state linear regression and (b) a second-order ARX model on the scalar output, checking that both recover the dynamics, one up to nothing (states observed, exact coordinates) and one up to input-output equivalence. We then evaluate the fitted full-state model recursively on overlapping starts from the held-out tail, under the recorded exogenous action sequence. Comparing its horizon-30 error with the true-model floor isolates the excess caused by estimation imperfection; the code reports that gap explicitly rather than asking you to infer it from the plot.

Let's code this. Everything numerical in this chapter uses NumPy first, PyTorch only in the neural section.

The one-step curve is nearly flat with horizon because it is really the same measurement repeated: each evaluation restarts from truth, discards the model's previous guess entirely, and asks a fresh single question. Its flatness across "horizons" is therefore an artifact of the protocol, not a property of the model. In this seeded run the recursive curve climbs steadily: at horizon ten the model has consumed nine of its own guesses, each contributing error that the light damping preserves. The dashed true-model reference is nonzero because deterministic mean rollout cannot reproduce the realized process disturbances. The printed horizon-30 difference separates the fitted model's error from that run-specific reference; its positive value is an empirical observation for this seed, not a build invariant or a theorem about every finite record. The printed count also makes clear that the curve averages overlapping starts from the held-out tail. This experiment evaluates recursive use of the fitted full-state model under recorded controls, not an ARX rollout or a planner-generated action sequence. Part XI returns to these evaluation distinctions systematically.

Nonparametric methods refuse, up front, to fix a finite parameterization. Instead of asking "which of our candidate transfer functions is it?", they ask "what can the data tell us about the response function directly?" The question being dodged is the riskiest one in parametric practice: order selection, choosing how many poles, lags, or basis functions the model possesses. Guess too low and the model cannot represent the system regardless of optimization effort; guess too high and variance explodes with little improvement in fit. Nonparametric methods postpone or dissolve that commitment entirely, letting the data decide the effective complexity.

The oldest and most useful instance is the impulse response: the output sequence caused by a momentary unit input kick, $G(k) = C A^k B$ for a linear system in state-space form. With the convention $x_{t+1}=Ax_t+Bu_t$, $y_t=Cx_t$, zero initial state, and no direct feedthrough, the input-driven output is $y_{t+1}=\sum_{k=0}^{\infty}G(k)u_{t-k}$; equivalently, an unshifted impulse response places $CB$ at lag one and a direct-feedthrough term $D$ at lag zero. Absolute summability of the impulse-response weights is the input-output, or BIBO, stability condition that permits finite truncation with a bounded tail. For a minimal realization this agrees with internal asymptotic stability, but a nonminimal realization can hide unstable unobservable or uncontrollable modes. In practice the discarded tail is not literally zero: it is accepted only after its bound or estimated contribution falls below a declared numerical or noise tolerance. Visible unstable modes do not decay and therefore defeat this truncation argument. Estimating $G(k)$ from data is regression under linearity and a justified finite-horizon approximation: stack lagged inputs into the regressor and let least squares fill in the weights. Two properties follow immediately. Each tap is estimated independently, so a bad estimate at one lag says nothing about neighbors, and no finite state-space order has been imposed. The frequency response is the Fourier transform of the impulse response, the object engineers have long probed by driving systems with swept sinusoids and reading amplitudes and phases off instruments. Frequency-domain plotting carries real diagnostic power: resonant peaks mark lightly damped modes, phase crossings flag stability margins, and rolls in magnitude reveal bandwidth. Impulse- and frequency-domain estimators remain industrial defaults precisely because they answer "what does the system do?" without first demanding "what kind of system is it?"

The modern incarnation regularizes the impulse response toward smoothness. Under open-loop exogeneity, raw finite-impulse-response least squares is unbiased when the chosen truncation contains the true response, or when omitted-tail effects are orthogonal to the included regressors; otherwise truncation itself can introduce bias. Even in the favorable case, long-lag estimates are noisy because the true weights are small while the estimation difficulty is unchanged. Stable-kernel Gaussian-process methods address this dilemma by choosing a covariance construction whose induced impulse-response estimates satisfy the intended stability and decay conditions; the guarantee comes from that kernel construction, not from Gaussian-process regression in general. This gives a principled bias-variance dial: the kernel encodes qualitative prior knowledge without committing to a transfer-function order. The kernel choice occupies the role that model-class choice played in parametric work, but its consequences are typically softer than fixing an incorrect finite order.

Where nonparametrics shine and struggle deserves honest accounting on both sides:

• They excel as a first look, as the highest-fidelity available estimate when no reliable structure is known, and for validating parametric models afterward. A parametric model whose predicted frequency response deviates visibly from a kernel-smoothed nonparametric estimate has been caught misbehaving without the validator ever having committed to the parametric assumptions.
• They struggle for nonlinear systems, where convolution ceases to describe behavior, and for control design, where most synthesis theory wants structured matrices, not ten thousand FIR taps. Convolution is a statement about superposition, and nonlinear systems break superposition at the root; extending nonparametric ideas to the nonlinear world generally forces you to reintroduce some structural commitment, at which point much of the appeal evaporates.

It helps to hold parametric and nonparametric identification as two ends of the same bias-variance tradeoff rather than rival camps. Parametric estimators inject strong assumptions (linearity, order, form of the nonlinearity) and win dramatically when assumptions hold, collapsing variance; nonparametric ones hedge, paying variance for protection against misspecification. Modern practice mixes freely, for example estimating an impulse response nonparametrically to decide the model order, then fitting a parametric state-space model of that order and validating against the nonparametric estimate. The same dialectic reappears between physical priors and deep networks later in this chapter, with networks playing the flexible high-variance extreme and grey-box structure playing the biased-but-efficient one.

Subspace methods occupy a sweet spot: they are noniterative (no nonlinear optimization), work directly from input-output data, and return state-space matrices in a similarity-transform sense. The importance of "noniterative" is easy to undersell. PEM-style search needs initialization, can stagnate in poor local minima, and offers no certificate that another restart would not find something better; a direct linear-algebraic route runs in one deterministic pass, making it attractive for large industrial problems and for bootstrapping the initialization that iterative methods still require. Their intellectual root is realization theory: Ho and Kalman's 1966 realization algorithm showed that if you know a system's impulse response exactly, the dynamics follow from linear algebra on a Hankel matrix, with no search involved. Realization theory answered a foundational question: given input-output behavior, does an internal state-space description exist, how big can it be, and is it unique? The answers, existence for any realizable response, minimality equivalent to simultaneous controllability and observability, uniqueness only up to similarity transform, match the identifiability discussion above exactly.

The core construction is worth seeing in detail because it demystifies what "learning dynamics" can mean when done geometrically. Rather than searching for parameters that make predictions match, subspace methods exploit rigid internal structure that the data matrix must exhibit if the underlying system is linear, and read the dynamics off that structure like pulling a fingerprint off a surface.

1. Estimate (or, in simulation, compute) the first $2s-1$ strictly proper Markov parameters in the shifted sequence $G(k)=CA^kB$, $k = 0, \dots, 2s-2$, for an $s$-by-$s$ block Hankel matrix. Here $G(0)=CB$ denotes the response one state transition after the impulse. In the unshifted impulse-response convention used earlier, the lag-zero term is the direct feedthrough $D$ (zero for a strictly proper system), and $CB$ is the lag-one term. Naming the shift explicitly keeps the two conventions consistent. Everything that follows derives from the repetitive multiplication structure baked into this formula.
2. Arrange them into the Hankel matrix $H$, whose $(i,j)$ entry is $G(i+j)$:

Because each block factorizes through $A^i$, the Hankel matrix has rank at most $n$: each entry can be written as $G(i+j) = (C A^i)(A^j B)$, splitting cleanly into a left factor determined only by the row index $i$ and a right factor determined only by the column index $j$. Stacking the left factors row-wise gives the extended observability matrix $\mathcal{O}_s$ (blocks $C A^i$, $i = 0, \dots, s-1$), and stacking the right factors column-wise gives the controllability matrix $\mathcal{C}_s$ (blocks $A^j B$, $j = 0, \dots, s-1$), so $H = \mathcal{O}_s \, \mathcal{C}_s$: however many rows and columns the table has, its rank cannot exceed the inner dimension $n$. This is the central insight. An arbitrarily large table of response numbers carries only $n$ independent dimensions, arranged as a product of an observability object times a controllability object; the dimension of the hidden state is literally visible as the rank deficiency of an ordinary matrix. All of the dynamics live in a low-rank bilinear factorization of this data matrix. 3. Compute the truncated SVD, $H \approx U_n \Sigma_n V_n^\top$. The rank-$n$ truncation both discovers the system order (singular value profile) and hands you balanced factors: $\hat{\mathcal{O}} = U_n \Sigma_n^{1/2}$ acts as an extended observability matrix, $\hat{\mathcal{C}} = \Sigma_n^{1/2} V_n^\top$ as its controllability counterpart, with the shift structure $\hat{\mathcal{O}}_{\text{down}} \, \hat{A} \approx \hat{\mathcal{O}}_{\text{up}}$, where $\hat{\mathcal{O}}_{\text{down}}$ means $\hat{\mathcal{O}}$ with its bottom block row discarded and $\hat{\mathcal{O}}_{\text{up}}$ means $\hat{\mathcal{O}}$ with its top block row discarded; multiplication by $\hat{A}$ advances every power, turning $C A^i$ into $C A^{i+1}$. Note the double duty performed here: order detection and factor extraction come from a single decomposition, with no separate hypothesis-testing machinery needed for the former. 4. Recover $\hat A$ by solving the overdetermined shift invariance: drop the last block row of $\hat{\mathcal{O}}$ and the first to get $\hat{\mathcal{O}}_1, \hat{\mathcal{O}}_2$, then $\hat A = \hat{\mathcal{O}}_1^\dagger \hat{\mathcal{O}}_2$; read $\hat B$ from the first block column of $\hat{\mathcal{C}}$ and $\hat C$ from the first block row of $\hat{\mathcal{O}}$. The pseudoinverse is appropriate precisely because $\hat{\mathcal{O}}_1$ is tall and overdetermined: many rows of shifted observability data vote on the single $\hat A$ that best explains their agreement, which is also what confers noise robustness on this step. The word balanced refers to the particular similarity transform chosen by the SVD, one that scales states equally by controllability and observability, numerically desirable and otherwise arbitrary, another reminder that coordinates are conventions. 5. Handle noise. Real data contain measurement noise, so stochastic subspace variants separate input effects and noise before the SVD. Classic MOESP formulations use LQ/QR factorizations and orthogonal projections relative to the future-input subspace, while N4SID is commonly expressed through oblique projections of future outputs onto past data along future inputs. Both replace raw, tap-wise impulse-response estimates with structured projections of stacked measured windows; the conceptual skeleton remains unchanged. Persistently exciting inputs are essential because the projected Hankel matrices must retain the rank needed to reveal the system order.

Let us verify the recipe numerically, pretending the impulse response came from clean measurement.

Relatedly, eigenvalue problems on structured data matrices are reappearing inside modern research on learning Koopman operators and linear encodings of nonlinear systems, topics we take up when discussing architectures in Part V. The lineage is direct: lift a nonlinear system into coordinates where it acts linearly, and suddenly Ho-Kalman-style algebra replaces iterative optimization, trading expressive power for the same geometric cleanliness enjoyed here.

Everything so far leaned on structure we rarely enjoy: global linearity, low order, scalar outputs. The interesting worlds are nonlinear, high-dimensional, and partially observed. Neural system identification replaces the hand-designed model class with a learned function, typically producing what modern readers would simply call a world model's dynamics module. This equivalence is worth stating bluntly: if you take the classical identification pipeline, swap the linear predictor for a network, the lagged-history feature map for a learned encoder, and the squared one-step loss for a probabilistic likelihood, you have reinvented the core architecture of contemporary world modeling. The main changes are scale, generality of the observation spaces, and engineering discipline around training. The conceptual issues (excitation, compounding, and noise-model specification) carry over, but their classical remedies require adaptation rather than transferring unchanged.

The dominant architecture keeps a recurrent latent state $z_t$ instead of the true state $x_t$, and learns three maps jointly from raw trajectories. The division of labor mirrors the classical pipeline exactly, which is why classical vocabulary remains useful for diagnosing failures:

• Encoder $q(z_t \mid o_{\le t}, u_{<t})$: compresses the observation-and-action history into a latent summary, playing the role of the belief state from our filtering chapter, but parameterized by a recurrent network or transformer rather than a Bayesian recursion. Where the Bayesian belief tracked an exact posterior distribution, the encoder learns whatever compression suffices to support accurate prediction, which may discard details irrelevant to future outcomes.
• Transition model $\hat z_{t+1} = \hat f_\theta(z_t, u_t)$: the actual dynamics hypothesis, trained to predict its own next latent state (often as a distribution, following the probabilistic styling of our filtering discussion). This component plays precisely the role $f_\star$ played in our setup section; everything the classical theory said about the transition function's responsibilities applies here.
• Decoder $\hat o_t = \hat g_\phi(z_t)$: ties latents back to measurable observations so the whole thing is supervised by prediction of real data. When the predictive likelihood is expressed through this decoder, latents that do not support accurate observational predictions cannot lower the predictive loss, however tidy their internal dynamics. An explicit reconstruction loss strengthens that grounding but is an optional auxiliary term in the objective below.

At the conceptual level, training targets a one-step or few-step predictive log-likelihood over transitions drawn from the data-collecting policy, such as $\min_{\theta,\phi} \sum_t -\log p(o_{t+1} \mid o_{\le t}, u_{\le t})$, optionally with auxiliary terms, reconstruction losses, and regularization. Latent-variable implementations usually optimize a variational bound or sampled approximation rather than this exact marginal likelihood.

where:

• $o_t$: the raw observation received at time $t$, whatever the sensors actually stream
• $o_{\le t}$: the entire observation history up to and including time $t$
• $p(o_{t+1} \mid o_{\le t}, u_{\le t})$: the model's predictive distribution over the next observation given past observations and the action history through the latest command
• $\theta, \phi$: trainable parameters of the encoder and transition model ($\theta$) and of the decoder ($\phi$). Architectural choices (recurrent networks, gated units, attention, state-space layers like S4/Mamba-style selective SSMs) supply inductive biases about memory and timescale; we survey these in Part III, Chapters 3 and 4, and the full lineage (emergent latent dynamics models such as E2C, PlaNet, Dreamer, and successors) belongs to Part V and Part VIII.

Reading the objective against the classical checklist makes several things visible at once. The loss is one-step (at best few-step), so it grades the model exactly where Figure @fig-rollout-compounding warned us grading is cheapest and decisions demand the most. The noise model is implicitly "outputs are contaminated, dynamics may be treated as flexible," inheriting ARX-style assumptions without anyone declaring them. And nothing in the objective penalizes confidently wrong predictions in regions the data-collecting policy never visited, because no data lives there to score. None of these observations amounts to a flaw in the recipe; each flags a place where additional machinery must be bolted on deliberately rather than inherited for free.

Neural identification also admits hybrid forms closer to classical spirit. Physics-informed and grey-box models hard-code trusted structure, estimate unknown physical parameters, and may learn a residual for whatever the trusted equations omit. Classical operator inference fits reduced linear or polynomial operators by least squares on prescribed lifted features; related Koopman-autoencoder variants learn the lifting itself. Implicit training on multistep rollout errors explicitly buys closed-loop quality at gradient cost. The grey-box strategy is especially attractive in engineering contexts where some physics is trustworthy and some is not: encoding conservation laws removes entire families of physically impossible predictions from consideration, shrinking what the data must teach and protecting the parts of the model that matter most for safety.

Four recurring failure surfaces deserve explicit treatment. Three have direct classical analogues, while optimization through long learned recurrences is a modern training issue that interacts with the older questions of excitation and conditioning:

• Compounding rollout error. Multistep autoregressive prediction feeds each guess back in, so small per-step bias can compound rapidly, especially through unstable modes. Mitigations include multistep objectives, scheduled teacher forcing, and per-step spectral or Lipschitz constraints. These constraints bound local amplification, but the realized rollout rate depends on products of Jacobians—the associated Lyapunov exponents or joint spectral behavior—not on the eigenvalues of any single Jacobian alone.
• Distribution shift and model exploitation. Data collected under one policy trains a model accurate near that policy's trajectory tube. A planner then queries the model off-tube, where it is silently wrong, acts foolishly, and steers the next data collection somewhere even stranger. A decision procedure can actively seek out the region where model validity lapses, because optimization concentrates wherever errors happen to favor unrealistically optimistic outcomes. Model-based reinforcement learning (Part II, Chapter 6) turns this cycle into an algorithm; the epistemic lesson belongs in every practitioner's head regardless.
• Similarity-transform ambiguity. Latent coordinates are identified only up to transformations preserving behavior, so latents need not align with interpretable physical quantities. If you followed the linear identifiability discussion, you will recognize this instantly as the structural ambiguity from decades ago, generalized to arbitrary learned encodings, and with the same fundamental status: it reflects what input-output information can determine, not a defect of any particular training run. For planning this is harmless; for scientific modeling it is a reproducibility hazard, since two laboratories fitting the same data can obtain latents differing by an arbitrary reparameterization yet making identical predictions, defeating naive comparisons of learned variables across studies.
• Optimization pathologies through time. Backpropagating through hundreds of recurrence steps meets vanishing or exploding gradients, while rare events contribute little loss mass when the data under-excites them. Classical excitation discipline helps diagnose the coverage problem, but it does not by itself solve long-horizon optimization. A nominally excellent model can therefore remain blind to exactly the transitions a cautious planner most wants predicted accurately.

To make the stakes concrete, we study a forced damped pendulum, smooth and innocent-looking, but governed by $\sin$ of position. The pendulum is a useful adversary for a linear predictor on the native variables $(\theta,\omega,u)$: no affine map on those variables can reproduce the sine term exactly over a nondegenerate interval. A lifted linear model containing $\sin(\theta)$ could represent the Euler step, so the limitation belongs to the chosen features rather than to every conceivable linear representation. We fit a linear model and a small feedforward network on identical lagged feature definitions and rollout protocols, using separate independently generated training and validation trajectories.

First the environment and data. The state is $(\theta, \omega)$: angular position and velocity; the action is applied torque.

In this seeded run both curves increase overall with horizon, and the structure-blind linear regression stays ahead across all tested horizons. That run-specific outcome illustrates how a misspecified but well-conditioned estimator can beat a more expressive network trained only on a modest one-step dataset. The logarithmic axis makes multiplicative differences readable without implying that it magnifies them; the printed and plotted values determine the size of the gap. The oscillatory mismatch likely includes phase error, but this experiment does not decompose the discrepancy by cause. Notice too what the exercise did not require: reward signals, interaction, or a planner. Neural identification here is pure supervised sequence learning. The gap between "good sequence model" and "good basis for deciding" is precisely the space of the coming chapters: Classical Planning and Optimal Control will consume models like these, and Model-Based Reinforcement Learning will interrogate them adversarially.

No identified model outruns its data. This deserves more than a bullet list, so let us dwell on the deeper issues.

A central failure mode is confident extrapolation into regions the data never visited. On-distribution fit diagnostics need not flag it because they contain no evidence about those regions. For example, a low-order lift-curve fit around cruise can produce a smooth but physically invalid extrapolation at a 60-degree angle of attack. When such a model drives a planner, the planner may select actions that look favorable precisely where the model is least reliable, then collect new data in that poorly supported region. Aggressive optimization can therefore preferentially select favorable model errors. This is not inevitable for every imperfect model and optimizer, but it is a practical reason to regularize optimism, maintain uncertainty over plausible models, keep planners within validated envelopes, and collect data where the current model admits ignorance. Ensemble disagreement is one inexpensive warning signal, although shared data and inductive biases can make the members confidently agree for the same wrong reason.

The evaluation question is likewise subtler than "low prediction loss." Responsible practice reports, at minimum: one-step metrics (did the model class capture local dynamics?); recursive rollout metrics under realistic action distributions (can anything plan with it?); behavior under held-out excitation regimes and held-out initial conditions (is it identifiable-level good, or merely fitted?); closed-loop performance with a controller designed against it (the ultimate, circular-in-a-good-way test); and calibration of predictive uncertainty (does the 95% interval contain truth about 95% of the time?). Reading the list as a sequence of escalating claims clarifies its logic: one-step competence certifies local fidelity; rollout competence certifies compositional fidelity over the tested horizon; held-out excitation certifies that accuracy survived a change of behavioral regime; closed-loop success certifies usefulness under the feedback the model induces; calibration certifies honesty about its own limitations. Each level licenses decisions the previous one does not. Comparison against the irreducible floor from process noise keeps expectations honest. And because high-dimensional perceptual models can look excellent while missing task-critical abstractions entirely, task-grounded evaluations, measuring whether plans formed from the model succeed, belong beside perplexity-style scores, a theme elaborated when we treat evaluation formally in Part XI. The separation between predictive fidelity and decision usefulness reaches its sharpest form here: a video-prediction model may score superbly on pixel likelihood while a planner relying on it misses obstacles, because the lost abstractions never carried much pixel mass.

On impact: system identification has a long industrial record in flight control, process control, automotive powertrains, and spacecraft attitude systems. The methods and validation practices vary by domain, so those examples are evidence of broad use rather than a claim that every deployed controller follows one recipe. Neural and latent approaches expand the observation spaces and parameter counts that identification can address, while inheriting the same need to state the operating regime, excitation coverage, and noise assumptions under which a model was validated.

Three additional caveats round out the picture:

• Closed-loop coupling. Identification under feedback quietly changes what the data can certify, as noted above; whenever you retrofit a model onto logs whose recorded actions came from a human or algorithmic policy reacting to observations, you are doing closed-loop identification whether or not you say so. Driving logs, interaction traces, and other action-annotated behavioral records often have this policy-dependent structure, so the caveat applies more widely than the word "control" suggests.
• Nonstationarity and adaptation. Systems drift, wear, and reconfigure. Online and adaptive identification extends these estimators with forgetting factors and recursive updates, blurring into the continual-learning questions of Part VI. Recursive least squares with forgetting, for instance, runs the same regression incrementally while exponentially discounting old evidence, converting a stationary batch estimator into one that tracks slow drift, at the price of re-admitting the noise sensitivity that averaging used to suppress.
• Partial observability. When sensors see projections of hidden state, identification of reduced-order models risks conflating modeling error with latent-variable mismatch; resolving the two is the business of Part III. A poorly predicting reduced model raises two confounded diagnoses: the reduction kept too few dimensions, or the prediction mechanism is wrong within the kept ones, and distinguishing them requires tools from the representation-learning chapters ahead.

The key parameters for the identification implementations are:

• proc_std: Standard deviation of injected process noise (0.01) in the drone trajectory, including its held-out tail. It sets a per-step expected squared-error contribution; the irreducible multi-step rollout error is horizon dependent because new disturbances accumulate along the trajectory.
• r, c: Row and column counts of the Hankel matrix (both 8). For a minimal order-$n$ deterministic system, $r \ge n$ and $c \ge n$ support a full-rank Hankel factorization under sufficient excitation, while recovering the state-transition operator from a shifted Hankel requires enough adjacent blocks. Here 8 is comfortably above the required depth and deliberately generous for a numerical-rank check; it is not a minimal or "twice the order" prescription.
• sv_threshold: Relative singular-value cutoff (1e-8) declaring numerical rank and thereby the estimated system order.
• Excitation gains (0.8, 0.6): Shaping-filter coefficients coloring the simulation input; their choice governs whether persistent excitation holds.
• ARX lags: Number of lagged outputs and inputs entering the regressor (two each here); this fixes the asserted model order.
• DT: Integration timestep of the pendulum environment (0.05), tying discrete transitions to continuous dynamics.
• hidden_dims: Widths of the feedforward dynamics network (6 -> 64 -> 64 -> 2 with Tanh units).
• lr: Adam learning rate (1e-3) used for the network fit.
• epochs: Number of full-batch training iterations (400).

• System identification learns the transition model from trajectories, complementing filtering, which infers state given a known model. Together they form the estimation backbone of Part II.
• Identifiability is a joint property of system, class, and data. Structural ambiguities like similarity transforms are unavoidable facts; practical unidentifiability from under-excited datasets is a solvable experiment-design problem.
• Persistent excitation formalizes input richness. With a correctly specified identifiable model, suitable exogeneity, stability and regularity conditions, sufficiently rich inputs support convergence to a unique probability limit. Under feedback, independently rich dither can restore missing excitation, but dither alone does not repair model misspecification or endogenous noise.
• Prediction-error methods generalize least squares to arbitrary model classes by optimizing predictive likelihood. Consistency with the true system requires structural containment, identifiability, informative data, and suitable regularity conditions; correct noise specification additionally supports the maximum-likelihood interpretation and asymptotic efficiency.
• Nonparametric methods estimate impulse and frequency responses directly, trading structure for robustness, and modern kernel methods regularize them with stability-aware priors.
• Subspace methods realize dynamics from data geometry, reading state-space matrices out of a low-rank Hankel factorization with no iterative optimization; rank discovery reveals system order en route.
• Neural identification replaces function-class choice with learned networks, buying expressive power for nonlinear partially observed worlds while paying in compounding rollout error, exploitation risk, and lost interpretability.
• Grading differs from using: one-step loss, recursive rollout error, off-distribution probes, closed-loop success, and calibration each certify a different claim about a world model.

You now possess the tools to turn trajectories into models, and the vocabulary to ask the uncomfortable question: what has this dataset certified about them? The next chapter, Classical Planning and Optimal Control, puts those models to work. Given a roughly correct, non-exploitable, usefully excited identified model, we will ask how to extract decisions from it: dynamic programming over states, trajectory optimization through continuous spaces, and the stability and optimality bounds that classical control attaches to model-derived policies. Those qualifiers are not throwaway hedging; each corresponds to a concrete risk cataloged in this chapter. Knowing which qualifier the evidence supports is often the difference between a planner that ships and a planner that surprises.

Ready to test your understanding? Take this quick quiz to reinforce what you've learned about system identification.