Preferential Subspace Identification (PSID)

• PSID is a linear system identification approach that constructs state-space models from primary signals to predict secondary signals by extracting latent dynamics.
• It employs a two-stage SVD procedure to extract shared and residual latent states, ensuring robust parameter estimation and optimal prediction.
• Extended with filtering and smoothing algorithms, PSID achieves minimal MSE in secondary signal estimation, notably in applications like neural decoding.

Preferential Subspace Identification (PSID) is a linear system identification methodology designed for multivariate time series, particularly where a primary signal (such as neural activity) is used to predict a secondary signal (such as behavior). PSID constructs a state-space model that extracts latent dynamical structure from a "primary" observation sequence to achieve optimal prediction of a "secondary" sequence. Originally developed for optimal prediction given past information ("innovations form"), recent advances have extended PSID to incorporate optimal filtering (using concurrent data) and forward–backward smoothing (using all available data), leveraging both primary and secondary signals for enhanced dynamic state estimation and decoding (Sani et al., 21 Jul 2025).

1. State-Space Model and Core PSID Algorithm

In PSID, the underlying dynamics are represented by a forward stochastic state-space model: $$\begin{cases} x_{k+1} = A x_k + w_k,\ y_k = C_y x_k + v_k,\ z_k = C_z x_k + \eta_k, \end{cases}$$ where $y_k \in \R^{n_y}$ (primary), $z_k \in \R^{n_z}$ (secondary), and $x_k \in \R^{n_x}$ is the latent state. The process, observation, and secondary-signal noises $(w_k, v_k, \eta_k)$ are zero-mean Gaussian, with covariance structure: $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$ and $A, C_y, C_z$ are system, primary observation, and secondary observation matrices, respectively.

The original PSID aims to estimate a model in the Kalman predictor ("innovation") form, determining (a) the predictor for the dynamics of $y_k$ and (b) the low-dimensional shared subspace enabling optimal prediction of $z_k$ from the history of $y_k$. This is achieved via a two-stage subspace identification process:

• Stage 1 (shared subspace extraction): SVD of the projection of the future block Hankel matrix of $y_k \in \R^{n_y}$0 onto the row-space of the past block Hankel matrix of $y_k \in \R^{n_y}$1, selecting the top $y_k \in \R^{n_y}$2 modes and reconstructing the corresponding latent states.
• Stage 2 (residual latent extraction): Remove shared subspace influence from future blocks of primary data, perform SVD on reprojected components to extract the remaining $y_k \in \R^{n_y}$3 latent dimensions.
• System parameter estimation: Stack estimated latent trajectories, then estimate $y_k \in \R^{n_y}$4 (Kalman predictor gain) using ordinary least squares. Specifically,

$y_k \in \R^{n_y}$5

2. Filter Gain Identification and Reduced-Rank Regression

While the predictor-form gain $y_k \in \R^{n_y}$6 is uniquely determined by the observed data, the Kalman filter update gain $y_k \in \R^{n_y}$7 generally is not when the state and observation noises are correlated ($y_k \in \R^{n_y}$8). With access to the secondary signal during training, PSID uniquely identifies $y_k \in \R^{n_y}$9 by directly minimizing the filtering mean-squared error (MSE) on $z_k \in \R^{n_z}$0: $z_k \in \R^{n_z}$1 The problem is formulated as a reduced-rank regression (RRR): $z_k \in \R^{n_z}$2 where $z_k \in \R^{n_z}$3, $z_k \in \R^{n_z}$4, $z_k \in \R^{n_z}$5. The optimal $z_k \in \R^{n_z}$6 is computed, and under observability, $z_k \in \R^{n_z}$7 itself can be recovered if the concatenated observability matrix of $z_k \in \R^{n_z}$8 ($z_k \in \R^{n_z}$9) has a left-inverse.

A key property is that access to the secondary signal resolves non-uniqueness, uniquely determining $x_k \in \R^{n_x}$0 within the equivalence class of observable stochastic-form models, as formalized by Faurré’s theorem. In the limit of large sample size, the estimated $x_k \in \R^{n_x}$1 converges to the true value (Sani et al., 21 Jul 2025).

3. PSID with Filtering: Algorithmic Workflow

The PSID + filtering workflow involves the following steps:

1. Learn innovation-form model: Apply original PSID (SVD-based) to estimate $x_k \in \R^{n_x}$2.
2. Compute innovations: Calculate one-step-ahead innovations $x_k \in \R^{n_x}$3 and predicted states.
3. Form regression data matrices: $x_k \in \R^{n_x}$4, $x_k \in \R^{n_x}$5.
4. Reduced-rank regression: Solve $x_k \in \R^{n_x}$6, yielding $x_k \in \R^{n_x}$7.
5. (Optional) Gain recovery: If observability matrix $x_k \in \R^{n_x}$8 is invertible, recover $x_k \in \R^{n_x}$9 directly.
6. Filter application: For new data, apply

$(w_k, v_k, \eta_k)$0

This extension yields state estimates and filtered predictions of the secondary signal with optimal MSE performance (Sani et al., 21 Jul 2025).

4. Forward–Backward PSID Smoothing

Building on classical two-filter Kalman smoothers, PSID is extended to a forward–backward smoothing algorithm targeting the secondary signal:

• Forward pass: Run PSID + filtering on the entire sequence $(w_k, v_k, \eta_k)$1, yielding filtered secondary estimates $(w_k, v_k, \eta_k)$2 and residuals $(w_k, v_k, \eta_k)$3.
• Backward pass: Reverse and process $(w_k, v_k, \eta_k)$4, $(w_k, v_k, \eta_k)$5 using PSID + filtering, producing backward-filtered residuals $(w_k, v_k, \eta_k)$6, then reverse time to obtain $(w_k, v_k, \eta_k)$7.
• Smoothing combination: The final smoothed secondary estimate is given by

$(w_k, v_k, \eta_k)$8

which achieves the minimum MSE for the secondary signal over all available data, matching the theoretical optimum from a Rauch–Tung–Striebel (RTS) smoother (Sani et al., 21 Jul 2025).

5. Identifiability and Theoretical Guarantees

PSID’s identifiability properties differ from classical single-signal subspace identification:

• Single-signal case ($(w_k, v_k, \eta_k)$9): The filter gain $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$0 is non-identifiable; filtering yields the observation itself.
• Two-signal case ($\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$1 primary, $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$2 secondary): The joint statistics of $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$3 uniquely determine the optimal filter $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$4 minimizing the filtering MSE for $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$5. This secondary signal removes the stochastic-form equivariance present in the single-signal case.
• Parameter identifiability: The predictor-form parameters $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$6 are generically identifiable up to similarity transforms. The RRR step identifies the external characteristic $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$7; under observability conditions, $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$8 is also identified.
• Consistency: As the number of samples $\E\left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right] \left[ \begin{smallmatrix} w_k \ v_k \end{smallmatrix} \right]^T = \begin{pmatrix} Q & S \ S^T & R \end{pmatrix}, \quad \E[\eta_k \eta_k^T]=R_z$9, least-squares and SVD estimates converge (in Frobenius norm) to ground truth system parameters; the RRR solution for $A, C_y, C_z$0 converges as well (Sani et al., 21 Jul 2025).

6. Empirical Validation and Performance

PSID + filtering and smoothing methods were empirically validated on simulated state-space models, with random choices of model dimensions $A, C_y, C_z$1 and training samples up to $A, C_y, C_z$2:

• Parameter recovery: Model parameters $A, C_y, C_z$3 were recovered with normalized Frobenius error below 1% for identifiable parameters. As expected, non-identifiable $A, C_y, C_z$4 did not converge in the single-signal case.
• Decomposition performance: On held-out data, one-step-ahead prediction $A, C_y, C_z$5, filtered estimate $A, C_y, C_z$6, and smoothed estimate $A, C_y, C_z$7 achieved coefficient of determination ($A, C_y, C_z$8) values matching the performance of predictions from the true underlying model, including ideal filter and RTS smoother (Sani et al., 21 Jul 2025).

A plausible implication is that PSID + filtering provides a principled and practical means of identifying dynamic models and constructing optimal decoders for two-signal systems, particularly in the analysis of neural and behavioral data streams.

7. Significance and Applications

PSID and its extensions address the broader challenge of inferring dynamical latent representations from multivariate time series where joint prediction, filtering, and smoothing of secondary variables are essential—for example, neural decoding of behavior. The unique identifiability and empirical optimality guarantees position PSID as a rigorous tool for probing dynamic interactions and constructing evidence-based decoders in neuroengineering and related disciplines (Sani et al., 21 Jul 2025).