This article provides a comprehensive explanation of subspace identification methods — a class of algorithms that estimate state-space models directly from input-output data using linear algebra. It covers the fundamental idea behind subspace methods, the three major algorithms (N4SID, MOESP, CVA), model order selection, and the connection to the author's research on periodically time-varying and multirate system identification. Related articles, research papers, and MATLAB links are placed at the bottom.

Author: Hiroshi Okajima, Associate Professor, Kumamoto University, Japan — 20 years of control engineering research

For the overview of system identification methods, see the hub article: System Identification: From Data to Dynamical Models — A Comprehensive Guide

Why Subspace Identification Matters

In model-based control engineering, obtaining an accurate state-space model of the plant is essential for controller and observer design. Subspace identification methods provide a powerful, non-iterative approach to this problem: they estimate the state-space matrices $A, B, C, D$ directly from measured input-output data using singular value decomposition (SVD) and least squares — without requiring iterative optimization.

Subspace identification is widely used because:

It handles multivariable (MIMO) systems as naturally as single-input single-output (SISO) systems, with no additional complexity in the algorithm.
It avoids the local minima and convergence issues that plague iterative methods such as the prediction error method (PEM).
It provides a natural mechanism for model order selection through the singular values of the data matrices.
It produces state-space models directly, which are the natural starting point for state feedback, observer design, and model predictive control.
It serves as the foundation for advanced identification methods, including the cyclic reformulation-based approaches for periodically time-varying (LPTV) and multirate systems developed in the author's research.

The three major subspace identification algorithms — N4SID, MOESP, and CVA — share the same core idea but differ in how they weight the data. This article explains the common framework and the distinguishing features of each method.

Problem Setting

Consider a discrete-time linear time-invariant (LTI) system in innovations form:

$\displaystyle x(k+1) = Ax(k) + Bu(k) + Ke(k)$

$\displaystyle y(k) = Cx(k) + Du(k) + e(k)$

where $x(k) \in \mathbb{R}^{n}$ is the state, $u(k) \in \mathbb{R}^{m}$ is the input, $y(k) \in \mathbb{R}^{l}$ is the output, $e(k) \in \mathbb{R}^{l}$ is a zero-mean white noise innovation sequence, and $K$ is the Kalman gain. The system order $n$ is assumed to be unknown.

Goal: Given measured input-output data ${u(k), y(k)}_{k=0}^{N-1}$ , estimate the system matrices $A, B, C, D$ (and possibly $K$ ) and the system order $n$ .

The Core Idea: From Data to State Space via SVD

Block Hankel Matrices

The first step in any subspace method is to arrange the input-output data into block Hankel matrices. Given a block row parameter $i$ (a user-chosen integer, typically $i \gg n$ ), define the output block Hankel matrix:

$\displaystyle Y_{0|2i-1} = \begin{pmatrix} y(0) & y(1) & \cdots & y(j-1) \cr y(1) & y(2) & \cdots & y(j) \cr \vdots & \vdots & \ddots & \vdots \cr y(2i-1) & y(2i) & \cdots & y(2i+j-2) \end{pmatrix}$

where $j = N - 2i + 1$ is the number of columns. The input block Hankel matrix $U_{0|2i-1}$ is defined similarly. These matrices are then split into "past" and "future" blocks:

$\displaystyle Y_p = Y_{0|i-1}, \quad Y_f = Y_{i|2i-1}$

$\displaystyle U_p = U_{0|i-1}, \quad U_f = U_{i|2i-1}$

The subscripts "p" and "f" denote the past and future portions, respectively.

The Extended Observability Matrix

The key theoretical fact underlying all subspace methods is the following relationship. Under appropriate conditions, the column space of a certain weighted matrix derived from the future output data (after removing the input contribution) coincides with the column space of the extended observability matrix:

$\displaystyle \mathcal{O}_i = \begin{pmatrix} C \cr CA \cr CA^2 \cr \vdots \cr CA^{i-1} \end{pmatrix} \in \mathbb{R}^{il \times n}$

Since $\mathcal{O}_{i}$ has $n$ linearly independent columns (when the system is observable), the system order $n$ can be determined from the rank of this matrix — which in practice is estimated by examining the singular values of the data matrix.

The SVD Step

The core computational step is a singular value decomposition of an appropriately weighted and projected data matrix. The generic form is:

$\displaystyle W_1 \cdot \mathcal{O}_i \cdot X_f \cdot W_2 \approx U_s \Sigma_s V_s^T$

where $W_{1}$ and $W_{2}$ are weighting matrices (which differ between N4SID, MOESP, and CVA), and $X_{f}$ represents the future state sequence. The key observation is:

The matrix $\Sigma_{s}$ contains the significant singular values — their number equals the system order $n$ .
The left singular vectors $U_{s}$ span the column space of $\mathcal{O}_{i}$ , from which $C$ and $A$ can be extracted.

Once $\mathcal{O}_{i}$ is estimated, the system matrices are recovered:

$C$ is read off as the first $l$ rows of $\mathcal{O}_{i}$ .
$A$ is obtained from the shift structure of $\mathcal{O}_{i}$ : the submatrix formed by the first $(i-1)l$ rows, when left-multiplied by $A$ , yields the submatrix formed by the last $(i-1)l$ rows.
$B$ and $D$ are then estimated by least squares from the input-output data and the estimated state sequence.

The Three Major Algorithms

N4SID (Numerical Algorithms for Subspace State Space System IDentification)

N4SID, developed by Van Overschee and De Moor (1994), uses oblique projections of the future outputs onto the past data, along the direction of the future inputs. The oblique projection extracts the state sequence from the data by removing the contribution of future inputs while retaining the information from past inputs and outputs.

In the framework of the weighted SVD described above, the simplest variant of N4SID applies no additional weighting beyond the oblique projection itself:

$\displaystyle W_1 = I, \quad W_2 = I$

However, the original N4SID framework by Van Overschee and De Moor actually encompasses multiple weighting choices — the identity weighting is one specific case. MATLAB's n4sid function with the default 'auto' setting selects the weighting adaptively based on the data characteristics. The essential feature that distinguishes N4SID from the other methods is the use of the oblique projection, not the specific choice of $W_{1}$ and $W_{2}$ .

Key characteristics:

The oblique projection provides a direct estimate of the state sequence, from which the system matrices are obtained by least squares.
N4SID is numerically robust and widely implemented.
It is the default subspace method in MATLAB's System Identification Toolbox (n4sid function).

Reference: P. Van Overschee and B. De Moor, Subspace Identification for Linear Systems: Theory, Implementation, Applications, Kluwer Academic Publishers, 1996.

MOESP (Multivariable Output-Error State Space)

MOESP, developed by Verhaegen (1994), uses orthogonal projections (QR decomposition) to remove the input contribution from the output data, and then applies SVD to the projected output to estimate the extended observability matrix.

The key difference from N4SID is the projection strategy: MOESP first projects the output data onto the orthogonal complement of the row space of the future input data. This effectively removes the direct input-to-output feedthrough and reveals the state-related component.

In the weighted SVD framework, the MOESP weighting can be understood as follows: $W_{1} = I$ , while $W_{2}$ is implicitly determined by the orthogonal projection that removes the future input subspace. Specifically, let $\Pi_{U_{f}}^{\perp}$ denote the projection onto the orthogonal complement of the row space of $U_{f}$ . Then the SVD is applied to $Y_{f} \Pi_{U_{f}}^{\perp}$ (or an equivalent quantity obtained via QR factorization), which corresponds to setting $W_{2} = \Pi_{U_{f}}^{\perp}$ .

Key characteristics:

MOESP focuses on the deterministic part of the system (input-output dynamics), making it particularly effective when the noise model is of secondary interest.
It uses an instrumental variable interpretation, which can provide better statistical properties in certain noise conditions.
In MATLAB, MOESP is available through the n4sid function by setting the N4Weight option to 'MOESP' (via n4sidOptions).

Reference: M. Verhaegen, "Identification of the deterministic part of MIMO state space models given in innovations form from input-output data," Automatica, Vol. 30, pp. 61–74, 1994.

CVA (Canonical Variate Analysis)

CVA, developed by Larimore (1983, 1990), uses canonical correlations between past and future data to determine the system order and estimate the state. The weighting matrices are chosen to maximize the correlation between past and future, which provides an optimal statistical basis for distinguishing signal from noise.

The weighting corresponds to:

$\displaystyle W_1 = \Sigma_{ff}^{-1/2}, \quad W_2 = \Sigma_{pp}^{-1/2}$

where $\Sigma_{ff}$ is the sample covariance matrix of the future output data block $Y_{f}$ , and $\Sigma_{pp}$ is the sample covariance matrix of the past data block (combining $Y_{p}$ and $U_{p}$ ). These weightings normalize the data so that the SVD produces canonical correlations between past and future — values between 0 and 1 that directly indicate the statistical significance of each state component.

Key characteristics:

CVA provides the statistically optimal weighting for determining the system order — the canonical correlations directly indicate the relevance of each state component.
Model order selection via CVA is based on the canonical correlation coefficients, which have well-defined statistical distributions under certain conditions, enabling formal hypothesis testing.
CVA is particularly effective for stochastic systems where the noise plays a significant role.
In MATLAB, CVA is available by setting N4Weight to 'CVA' in n4sidOptions. MATLAB also uses CVA by default for frequency-domain data.

Reference: W.E. Larimore, "Canonical variate analysis in identification, filtering and adaptive control," Proc. 29th IEEE Conference on Decision and Control, pp. 596–604, 1990.

Comparison of the Three Methods

Feature	N4SID	MOESP	CVA
Projection	Oblique	Orthogonal (QR)	Canonical correlation
Weighting $W_{1}, W_{2}$	Identity (simplest variant)	Orthogonal complement of input subspace	Covariance-based normalization
Order selection	Singular values	Singular values	Canonical correlations
Noise handling	General	Focuses on deterministic part	Statistically optimal
Best suited for	General purpose	Deterministic-dominant systems	Stochastic systems
MATLAB option	`'auto'` (default)	`N4Weight = 'MOESP'`	`N4Weight = 'CVA'`

In practice, the differences between the three methods are often small for well-conditioned, moderately-sized datasets. MATLAB's n4sid function offers an 'all' option that runs all three methods and returns the one with the best fit.

Model Order Selection

A significant advantage of subspace methods is the built-in mechanism for model order selection. After the SVD step:

$\displaystyle \Sigma = \mathrm{diag}(\sigma_1, \sigma_2, \ldots, \sigma_{il})$

with $\sigma_{1} \geq \sigma_{2} \geq \cdots \geq \sigma_{il} \geq 0$ , the system order $n$ is determined by identifying a gap in the singular value spectrum: the first $n$ singular values are significantly larger than the remaining ones, which represent noise.

In practice, order selection methods include:

Visual inspection of the singular value bar chart — look for a clear drop-off after the $n$ -th singular value.
Information criteria such as the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC).
In MATLAB, specifying a range of orders (n4sid(data, 1:10)) generates a bar chart of singular values with an automatic order recommendation.

For example, in MATLAB:

% Load or create input-output data
data = iddata(y, u, Ts);

% Estimate with automatic order selection (range 1 to 10)
sys = n4sid(data, 1:10);

% The function displays a singular value plot
% and recommends the optimal order

MATLAB Implementation

MATLAB's System Identification Toolbox provides the n4sid function for subspace identification. Below is a basic workflow.

Basic Usage

%% Data preparation
Ts = 0.01;             % Sampling period
data = iddata(y, u, Ts);  % Create iddata object from output y and input u

%% Subspace identification with specified order
n = 3;                 % Desired model order
sys = n4sid(data, n);  % Identify a 3rd-order state-space model

%% Extract state-space matrices
[A, B, C, D] = ssdata(sys);

Choosing the Subspace Algorithm

%% Use MOESP weighting
opt = n4sidOptions('N4Weight', 'MOESP');
sys_moesp = n4sid(data, n, opt);

%% Use CVA weighting
opt = n4sidOptions('N4Weight', 'CVA');
sys_cva = n4sid(data, n, opt);

%% Try all methods and pick the best
opt = n4sidOptions('N4Weight', 'all');
sys_best = n4sid(data, n, opt);

Model Validation

%% Compare model output with measured data
compare(data, sys);

%% Check residuals (should be white noise)
resid(data, sys);

%% Bode plot of identified model
bode(sys);

Notes on MATLAB's n4sid Function

By default, n4sid estimates a discrete-time model with the same sample time as the data. To estimate a continuous-time model, set 'Ts' to 0.
The 'Feedthrough' option controls whether the $D$ matrix is estimated (default: not estimated, i.e., $D = 0$ ).
The 'DisturbanceModel' option controls whether the noise model (Kalman gain $K$ ) is estimated. Setting it to 'none' yields $K = 0$ .
The 'Focus' option in n4sidOptions determines whether the estimation minimizes prediction error ('prediction') or simulation error ('simulation').

Connection to Parametric Methods

Subspace identification and classical parametric methods (ARX, ARMAX, Output-Error, Box-Jenkins) represent two fundamentally different approaches to system identification.

Parametric methods assume a specific polynomial model structure (e.g., $A(q)y(k) = B(q)u(k) + C(q)e(k)$ ) and estimate the polynomial coefficients by minimizing a prediction error criterion. They require the user to specify the model orders and may involve iterative optimization.

Subspace methods bypass the polynomial structure entirely. They work directly with the state-space representation and use SVD rather than iterative optimization. The model order is determined from the data rather than specified a priori.

Aspect	Parametric (ARX/ARMAX)	Subspace (N4SID/MOESP/CVA)
Model structure	Polynomial (transfer function)	State-space
Order selection	User-specified	From singular values
Optimization	Iterative (PEM) or least squares (ARX)	Non-iterative (SVD + LS)
Local minima risk	Yes (for PEM)	No
MIMO handling	Complex	Natural
Noise model	Integral part of model structure	Optional

In MATLAB, ssest combines both approaches: it uses n4sid for initialization and then refines the estimate iteratively using a prediction error method.

For the detailed treatment of parametric methods and PEM, see: Classical Parametric System Identification: ARX, ARMAX, and the Prediction Error Method

Connection to Kernel-Based Identification

A more recent alternative to both parametric and subspace methods is kernel-based (regularized) identification, which estimates the impulse response directly in a reproducing kernel Hilbert space. Kernel methods avoid model order selection entirely by using regularization to control complexity.

Compared to subspace methods:

Subspace methods are best when a finite-dimensional state-space model is desired (e.g., for observer or controller design).
Kernel methods are best when the model order is uncertain or very high, or when a smooth nonparametric model is sufficient.
Both approaches can serve as complementary tools: subspace identification for model-based control design, and kernel methods for initial exploration or validation.

Connection to the Author's Research

Subspace identification is the key building block in the author's research on system identification for periodically time-varying (LPTV) and multirate systems.

Cyclic Reformulation-Based Identification for LPTV Systems

In the cyclic reformulation approach, an LPTV system with period $M$ is transformed into an equivalent higher-dimensional LTI system by constructing cycled signals. The cycled system is then identified using a standard subspace method such as N4SID. A subsequent state coordinate transformation recovers the original LPTV parameters from the identified cycled system matrices.

The fact that subspace methods produce state-space models directly (rather than polynomial models) is essential for this approach, because the cyclic reformulation structure is naturally expressed in state-space form.

For the detailed method, see the blog article: Cyclic Reformulation-Based System Identification for Periodically Time-Varying Systems

Paper: H. Okajima, Y. Fujimoto, H. Oku and H. Kondo, Cyclic Reformulation-Based System Identification for Periodically Time-Varying Systems, IEEE Access (2025)

Multirate System Identification

Multirate systems — where different sensors operate at different sampling rates — are a special case of LPTV systems. The identification algorithm again uses subspace identification (N4SID) as Step 2, applied to the cycled signals constructed from the multirate data. The subspace method identifies the cycled LTI system, and then the original time-invariant plant parameters are extracted via a state coordinate transformation.

For the detailed algorithm and numerical results, see the blog article: System Identification Under Multirate Sensing Environments

Paper: H. Okajima, R. Furukawa and N. Matsunaga, System Identification Under Multirate Sensing Environments, Journal of Robotics and Mechatronics, Vol. 37, No. 5, pp. 1102–1112 (2025) (Open Access)

The Role of Subspace Methods in the Identification Pipeline

The relationship between subspace identification and the author's cyclic reformulation approach can be summarized as follows:

Step 1 (Cyclic reformulation): Multirate/LPTV data → Cycled signals → LTI cycled system

Step 2 (Subspace identification): Apply N4SID to the cycled signals → Identified cycled model

Step 3 (Coordinate transformation): Recover original LPTV/multirate parameters

Subspace identification serves as the engine in Step 2, while the cyclic reformulation (Step 1) and coordinate transformation (Step 3) handle the multirate/LPTV-specific structure.

Key References

1. Van Overschee and B. De Moor, Subspace Identification for Linear Systems: Theory, Implementation, Applications, Kluwer Academic Publishers, 1996.
1. Verhaegen, "Identification of the deterministic part of MIMO state space models given in innovations form from input-output data," Automatica, Vol. 30, pp. 61–74, 1994.
W.E. Larimore, "Canonical variate analysis in identification, filtering and adaptive control," Proc. 29th IEEE CDC, pp. 596–604, 1990.
1. Ljung, System Identification: Theory for the User, 2nd ed., Prentice Hall, 1999.
MATLAB System Identification Toolbox documentation: n4sid

MATLAB Code

Code Ocean (Multi-Rate System): Multi-Rate System Code — Includes subspace identification as a key step
MATLAB File Exchange (Multi-Rate Observer): State Estimation under Multi-Rate Sensing: IEEE ACCESS 2023
GitHub (LMI basics): MATLAB Fundamental Control LMI
GitHub (System identification basics): MATLAB system identification

Self-Introduction

Hiroshi Okajima — Associate Professor, Graduate School of Science and Technology, Kumamoto University. Member of SICE, ISCIE, and IEEE.

www.control-theory.com / Blog / YouTube / GitHub

If you found this article helpful, please consider bookmarking or sharing it.

SubspaceIdentification #N4SID #MOESP #CVA #SystemIdentification #StateSpaceModel #MATLAB #ControlEngineering #MIMO #SingularValueDecomposition

Why Subspace Identification Matters

Problem Setting

The Core Idea: From Data to State Space via SVD

Block Hankel Matrices

The Extended Observability Matrix

The SVD Step

The Three Major Algorithms

N4SID (Numerical Algorithms for Subspace State Space System IDentification)

MOESP (Multivariable Output-Error State Space)

CVA (Canonical Variate Analysis)

Comparison of the Three Methods

Model Order Selection

MATLAB Implementation

Basic Usage

Choosing the Subspace Algorithm

Model Validation

Notes on MATLAB's n4sid Function

Connection to Parametric Methods

Connection to Kernel-Based Identification

Connection to the Author's Research

Cyclic Reformulation-Based Identification for LPTV Systems

Multirate System Identification

The Role of Subspace Methods in the Identification Pipeline

Related Articles and Videos

Blog Articles (blog.control-theory.com)

Research Web Pages (www.control-theory.com)

Video

Key References

MATLAB Code

SubspaceIdentification #N4SID #MOESP #CVA #SystemIdentification #StateSpaceModel #MATLAB #ControlEngineering #MIMO #SingularValueDecomposition