System Identification › Subspace Identification

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.

Compare the recovered eigenvalues with the true ones from earlier: identical to numerical precision. Also glance at the singular-value printout: with exact data, two singular values carry all the energy and the rest sit at machine zero, confirming the rank-two structure asserted by the factorization argument before any code ran. The dynamics were read off a rank-two factorization of a response table, with no gradient descent anywhere. The row slicing above is valid because this example has one output; a multi-output implementation must remove whole output-sized block rows, and similarly read whole input/output blocks for $B$ and $C$. Historically, realization ideas fed directly into ERA for aerospace modal testing and later into broader subspace-identification methods; process-industry practice also drew heavily on prediction-error methods rather than on Ho-Kalman alone.

The order-detection step of the recipe deserves its own picture: the singular-value spectrum of the Hankel matrix should show a sharp cliff at the true system order.

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.