Symmetric Linear Dynamical Systems are Learnable from Few Observations

Technical Summary: Symmetric Linear Dynamical Systems are Learnable from Few Observations

Problem Statement
The paper addresses the problem of learning the parameters of an $N$-dimensional stochastic linear dynamical system from a single trajectory of length $T$. The system evolves according to the equation $x_t = A x_{t-1} + \xi_{t-1}$, where $A \in \mathbb{R}^{N \times N}$ is a symmetric dynamic state matrix with a spectral radius $\rho(A) < 1$ (ensuring stability), and $\xi_t$ represents i.i.d. zero-mean Gaussian noise. The objective is to estimate the elements of $A$ under two regimes:

1. Full Observations: The entire state vector $x_t$ is observed.
2. Partial Observations: Only a subset of the state vector $x_{t,o}$ is observed, while the remaining components are hidden (latent).

The authors focus on recovering the matrix elements with guarantees on the maximum element-wise error ($\ell_\infty$ norm), a metric critical for applications like structure discovery and anomaly detection, rather than just the aggregate $\ell_2$ norm.

Methodology
The authors propose a new estimator based on the method of moments. Unlike classical approaches that rely on Ordinary Least Squares (OLS) or $\ell_1$-regularized methods (LASSO), this estimator does not require problem-specific regularization or assumptions about sparsity.

The core of the method exploits the temporal correlations inherent in the system. By iterating the system dynamics, the authors derive a relationship where the time-shifted covariance between states encodes powers of the matrix $A$:
$$E[x_t x_s^\top] = A^{s-t} E[x_t x_t^\top]$$

They define a family of estimators $\hat{S}_m(T)$ to recover powers of the dynamic matrix ($A^m$) directly from trajectory data:
$$\hat{S}_m(T) := \frac{1}{T-m} \sum_{t=0}^{T-m-1} x_t x_{t+m}^\top - \frac{1}{T-m-2} \sum_{t=0}^{T-m-3} x_t x_{t+m+2}^\top$$
This specific construction is designed to be asymptotically unbiased. The estimator is then applied to the observed coordinates in the partial observation setting to recover sub-matrices of the system dynamics.

Key Contributions and Theoretical Results

1. Logarithmic Sample Complexity for Full Observations:
   The primary theoretical contribution (Theorem 2) establishes that the proposed estimator recovers the matrix $A$ (and its powers) to a finite element-wise max-norm error using only $T = O(\log N)$ observations. This holds for both sparse and dense matrices. Consequently, the total number of data points ($O(N \log N)$) is significantly smaller than the number of unknown parameters ($O(N^2)$), a regime where traditional least-squares methods typically fail without strong sparsity assumptions.

2. Partial Observations and Sub-matrix Recovery:
   In the partially observed setting, the authors show that the sub-matrix $B$ (representing the internal dynamics of observed nodes) can also be recovered with $O(\log N)$ samples.

   • Recovering the coupling terms $CC^\top$ requires $O(N^2 \log N)$ samples.
   • Recovering the term $CEC^\top$ (involving hidden dynamics) requires $O(N^4 \log N)$ samples.
     The paper notes that recovering the individual matrices $C$ and $E$ is impossible up to a linear transformation due to the intrinsic symmetry of the problem (invariance under orthogonal transformations of the hidden space), but the derived block estimators for $CC^\top$ and $CEC^\top$ are identifiable.

3. No Regularization Required:
   A significant methodological contribution is that the estimator achieves these bounds without relying on $\ell_1$-regularization or incoherence assumptions, which are typically required for support recovery in sparse settings.

Experimental Results
The authors validate their theoretical bounds through numerical experiments on both fully and partially observed systems:

• Fully Observed: Experiments on both sparse graphs (degree 2) and dense structured graphs (star graphs) confirm that the proposed estimator achieves $T = O(\log N)$ scaling. In contrast, standard Least Squares (LS) exhibits linear scaling $O(N)$, and $\ell_1$-regularized LS (LS+L1) achieves logarithmic scaling only in the sparse case, failing in the dense case.
• Partially Observed: Experiments confirm that the observed sub-matrix $B$ scales as $O(\log N)$. While the theoretical bound for recovering coupling terms ($CC^\top$, $CEC^\top$) suggests higher polynomial dependencies, empirical results show a more favorable scaling of $O(N \log N)$ for these terms, outperforming the conservative theoretical predictions.

Significance and Claims
The paper claims that its new estimator provides a robust, regularization-free solution for learning symmetric linear dynamical systems from limited data. Its significance lies in:

• Efficiency: Achieving recovery with $T = O(\log N)$ observations, making it feasible to learn high-dimensional systems from short trajectories.
• Generality: The method works effectively for both sparse and dense matrices, overcoming the limitations of sparsity-dependent methods like LASSO.
• Application Relevance: The focus on element-wise error guarantees makes the method particularly suitable for structure discovery (identifying non-zero entries) and anomaly detection, where individual matrix elements are more critical than aggregate norms.

The authors conclude that while the current work focuses on stable, symmetric systems, the framework offers a promising direction for future exploration into marginally stable and non-symmetric systems.