# Learning Low-dimensional Latent Dynamics from High-dimensional Observations:
Non-asymptotics and Lower Bounds

Yuyang Zhang  Shahriar Talebi  Na Li

###### Abstract

In this paper, we focus on learning a linear time-invariant (LTI) model with
low-dimensional latent variables but high-dimensional observations. We provide
an algorithm that recovers the high-dimensional features, i.e. column space of
the observer, embeds the data into low dimensions and learns the low-
dimensional model parameters. Our algorithm enjoys a sample complexity
guarantee of order $\tilde{\mathcal{O}}(n/\epsilon^{2})$, where $n$ is the
observation dimension. We further establish a fundamental lower bound
indicating this complexity bound is optimal up to logarithmic factors and
dimension-independent constants. We show that this inevitable linear factor of
$n$ is due to the learning error of the observer’s column space in the
presence of high-dimensional noise. Extending our results, we consider a meta-
learning problem inspired by various real-world applications, where the
observer column space can be collectively learned from datasets of multiple
LTI systems. An end-to-end algorithm is then proposed, facilitating learning
LTI systems from a meta-dataset which breaks the sample complexity lower bound
in certain scenarios.

Machine Learning, ICML

##  1 Introduction

Analyzing high-dimensional time series data is essential for numerous real-
world applications in finance (Mudassir et al., 2020), economics (Maliar &
Maliar, 2015; Masini et al., 2023) and biology (Churchland et al., 2012;
Hajnal et al., 2023; Xia et al., 2021; Gallego et al., 2020; Stringer et al.,
2019). High-dimensional time series observations often find succinct
representation through a set of low-dimensional latent variables. In this
paper, the focus is to learn the low-dimensional dynamics capturing the very
essence of the time series, which becomes useful in various down-stream tasks
like prediction and inference (Churchland et al., 2012; Mudassir et al., 2020;
Pandarinath et al., 2018).

Popular techniques for such analysis range from linear models, such as linear
time-invariant (LTI) systems (Sikander & Prasad, 2015; Bui-Thanh et al., 2008;
Hespanha, 2018) and auto-regressive models (Dong et al., 2022; Poloni &
Sbrana, 2019; Qin, 2022), to more complex nonlinear models, exemplified by
recurrent neural networks (Yu et al., 2021; Sussillo & Barak, 2013; Medsker &
Jain, 1999). Herein, we focus on LTI models as they are interpretable and
require much less computational power in many real-world applications (Gallego
et al., 2020; Churchland et al., 2012; Dong et al., 2022).

Specifically, we learn partially observed LTI systems in the following form

$x_{t+1}=Ax_{t}+Bu_{t}+w_{t},\quad y_{t}=Cx_{t}+\eta_{t},$ (1)

where $y_{t}\in\mathbb{R}^{n}$ are _high-dimensional_ observations,
$x_{t}\in\mathbb{R}^{r}$ are low-dimensional latent variables with $n\gg r$,
$u_{t}\in\mathbb{R}^{m}$ are inputs and $w_{t},\eta_{t}$ are process and
observation noise, respectively. Matrices $A$ and $B$ are the system
parameters for the dynamics and inputs, respectively, and $C$ can be seen as
the observer, mapping latent variables to _noisy observations_. Recently, a
contemporary perspective has revitalized well-established system
identification (SYSID) algorithms for similar setups (Shirani Faradonbeh et
al., 2018; Sarkar & Rakhlin, 2019; Oymak & Ozay, 2019; Sarkar et al., 2021;
Zheng & Li, 2021). Without any structural assumptions, they have established
finite time convergence results with a rate of
$\tilde{\mathcal{O}}(\sqrt{n\cdot\text{poly}\left(r,m\right)/T})$, leading to
a sample complexity of
$\tilde{\mathcal{O}}(n\cdot\text{poly}\left(r,m\right)/\epsilon^{2})$ for an
$\epsilon$-well approximation of model parameters. Such results are not
satisfying for LTI systems with a large observation dimension $n$. One
alternative way to directly applying the SYSID algorithm is to initially
project the high-dimensional data to lower dimensions before conducting
further analysis (Saul et al., ; Churchland et al., 2012), either for
efficiency or for interpretability. Moreover, these projections may often
consist of “meta information” that are shared across different tasks. This is
consistent with many scenarios in meta-learning settings where we deal with
more complex datasets collected from similar observers but possibly with
different latent dynamics. The above facts motivate us to investigate the
provable efficiency of this type of projection method and its generalization
ability.

Contributions. In this paper, we study learning LTI systems with _high-
dimensional noisy observations_. We adopt a two-stage procedure to first
extract high-dimensional features, i.e. column space of the observer. We
subsequently perform standard SYSID on the resulting low-dim data and recover
the original model parameters. We establish sample complexity for this bifold
procedure that scales as
$\tilde{\mathcal{O}}([n+\text{poly}\left(r,m\right)]/\epsilon^{2})$ which
further reduces to only
$\tilde{\mathcal{O}}(\text{poly}\left(r,m\right)/\epsilon^{2})$ in the absence
of observation noise. This then naturally leads to the following question:

_Can this linear dependency on observer dimension $n$ be possibly improved by
some carefully designed algorithm?_

We show that the linear dependency on $n$ is unavoidable in the presence of
observation noise–as shown in our lower bound result Theorem 4.1, requiring at
least $\mathcal{O}(n/\epsilon^{2})$ samples from our learning problem.
Unfortunately, observation noises always exist for real-world applications, as
the sensors always involve uncertainties. The lower bound also indicates that
the sample complexity of the proposed algorithm is optimal up to some
logarithmic factors and dimension-independent constants.

Additionally, our pragmatic solution–namely, the separation of learning the
high-dimensional features (i.e. the observer column space) and learning the
rest of the model parameters– can be extended to meta-learning setups by
collecting metadata from a set of dynamical systems that share the same
observer model. By adopting a “leave-one-out” strategy for statistical
consistency, we show that as long as the metadata is _collectively_ large in
order of $\tilde{\mathcal{O}}(n/\epsilon^{2})$, we can successfully obtain an
$\epsilon$-well approximation of all the systems parameters involved in the
meta-data (Theorem 5.1). We finally note that such metadata with the same
observer is common in real-world applications (Hajnal et al., 2023; Xia et
al., 2021; Marks & Goard, 2021; Gallego et al., 2020). One example is in
neuroscience where neuron activities are measured over a long period by the
same set of electrodes or imaging devices. Although the subject may perform
different tasks or demonstrate behavior corresponding to different latent
dynamics, the observer model (i.e. the electrodes or the imaging device)
remains the same, resulting in the metadata considered here.

To summarize our contributions: Firstly, we provide Column Space Projection-
based SYSID (Col-SYSID) (Algorithm 1 in Section 3) that learns LTI systems
through high-dimensional observations with complexity
$\tilde{\mathcal{O}}(n/\epsilon^{2})$ (Theorem 3.4 in Section 3). As compared
with existing algorithms, ours reduces the polynomial dependency on the latent
dimension $r$ and the input dimension $m$. We then establish a sample
complexity lower bound for this problem, indicating the optimality of the
above result (Theorem 4.1 in Section 4). Our algorithmic idea is further
extended to a meta-learning setting, where we have access to datasets from
multiple similar systems. We provide an end-to-end framework for handling this
meta-dataset and learning all included systems (Algorithm 3 in Section 5).
With the help of the meta-dataset, we break the sample complexity lower bound
in certain scenarios (Theorem 5.1 in Section 5).

Before proceeding, we set the following notations throughout the paper.

Notations: Without further explanation, let $\delta$ be any probability in
$(0,1/e)$. We use $\text{poly}\left(\cdot\right)$ to denote polynomial and
logarithmic dependences. We use $\lesssim$ and $\gtrsim$ to hide absolute
constants. We use $\tilde{\mathcal{O}}(\cdot)$ to hide all problem-related
constants, absolute constants and logrithmic factors. Let $[N]$ be the set of
integers $\\{1,\dots,N\\}$. Let $\mathcal{M}_{[N]}$ denote the set
$\\{\mathcal{M}_{n}\\}_{n\in[N]}$ and let $\mathcal{M}_{-n}$ denote the set
$\\{\mathcal{M}_{n^{\prime}}\\}_{n^{\prime}\in[N]\setminus\\{n\\}}$ when the
full set $[N]$ is clear from the context. For any matrix
$M\in\mathbb{R}^{m\times n}$, let
$\sigma_{1}(M)\geq\sigma_{2}(M)\geq\cdots\geq\sigma_{\max\\{m,n\\}}(M)$ be its
singular values, and let $\sigma_{\min}(M)$ be the minimum non-zero singular
value. Let $\left\|M\right\|$ denote its operator norm, and $M{}^{\intercal}$
denote its transpose. Let $\text{col}(M)$ be the column space of $M$ and let
$\Phi_{M}$ denote any orthonormal matrix whose columns form a basis of
$\text{col}(M)$. For any orthonormal matrix $\Phi$, we use $\Phi^{\perp}$ to
denote the matrix such that $\begin{bmatrix}\Phi&\Phi^{\perp}\end{bmatrix}$ is
unitary, and we refer to $\Phi^{\perp}$ as the orthogonal complement of
$\Phi$. For any positive semi-definite matrix $\Sigma\in\mathbb{R}^{n\times
n}$, we slightly abuse the notation and let $\mathcal{N}(0,\Sigma)$ denote the
distribution of $\Sigma^{1/2}x$, where $x\sim\mathcal{N}(0,I_{n})$ follows the
standard Gaussian distribution.

###  1.1 Other Related Works

##### Other Linear Models with Low-Dimensional Dynamics.

Recently, there are lines of research on autoregressive models with low-
dimensional representations (Qin, 2022; Dong et al., 2022; Poloni & Sbrana,
2019). Indeed, we may also be able to reconstruct low-dimensional dynamics
from the model with the learned representations. However, sample complexity
results are not currently available for such models. Another related line of
research is on dynamic factor models (Hallin et al., 2023; Anderson et al.,
2022; Breitung & Eickmeier, 2006). Anderson et al. (2022) only provides
asymptotic convergence guarantees. Finite time convergence is indeed developed
in Hallin et al. (2023). However, their results can not imply dependence on
the system dimensions and are only about the convergence rates. Moreover,
certain assumptions might not be easily verified in our setting.

##### General SYSID Algorithms.

The most relevant sample complexity results are provided in recent SYSID
literature for learning partially observed LTI systems (Zheng & Li, 2021; Lee,
2022; Djehiche & Mazhar, 2022; Sarkar et al., 2021; Oymak & Ozay, 2019). The
last four papers focus on systems with stable latent dynamics, while the last
paper extends to unstable latent dynamics. Unfortunately, all algorithms have
suboptimal sample complexities. There also exist a line of work providing
sample complexity for learning fully observed systems (Sarkar & Rakhlin, 2019;
Shirani Faradonbeh et al., 2018), whose analysis techniques is applicable to
unstable dynamics. This paper focuses on stable latent dynamics, and we leave
for future directions to extend the results to unstable latent dynamics.

##### Existing SYSID Lower Bounds.

There exist literature providing lower bounds for learning partially observed
LTI systems. However, their results are not directly comparable to ours in our
setting. The most relevant papers are Sun et al. (2023); Fattahi (2021). The
former provides a lower bound depending on $r/n$, which is tailored for
systems with $r\gg n$. The latter provides a lower bound proportional to
$1/\epsilon^{4}$. The lower bound follows a slower rate than ours, though
being logarithmic in the system dimensions. Sun et al. (2022) provides lower
bounds defined for noise-free settings. Mao et al. (2021) provides lower
bounds for systems with $C=I$. Bakshi et al. (2023) provides lower bounds for
learning almost uncontrollable and unobservable systems. Other related works
include Djehiche & Mazhar (2021); Jedra & Proutiere (2023); Simchowitz et al.
(2018), which develop lower bound for fully observed systems.

##### General Subspace Learning Algorithms.

To learn the observer column space, several existing literature may be
related. Vaswani et al. (2018); Balzano et al. (2018); Candès et al. (2011);
Blanchard et al. (2007) focus on analysis for PCA subspace learning. However,
they assume that the dataset is i.i.d. sampled, which is not the case for
dynamic systems. Dynamic factor model techniques are also related (Hallin et
al., 2023). As discussed previously, the results can not imply the dependence
on the system dimensions. Other ideas include Deng et al. (2020), whose sample
complexity remains an open question.

##  2 Preliminaries and The Problem Setup

Consider linear dynamical systems in the form of Equation 1 with $x_{0}=0$
(for simplicity), process noises
$w_{t}\overset{i.i.d}{\sim}\mathcal{N}(0,\Sigma_{w})$ and isotropic
observation noises
$\eta_{t}\overset{i.i.d}{\sim}\mathcal{N}(0,\sigma_{\eta}^{2}I)$. Here
$\Sigma_{w}$ and $\sigma_{\eta}^{2}I$ are positive semi-definite matrices. We
denote such systems by
$\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\sigma_{\eta}^{2}I)$. In standard system
identification setup, given the input-output trajectory data from the system,
the objective is to learn system parameters up to the well-known similarity
transformation class, i.e., to learn tuple
$(\widehat{A},\widehat{B},\widehat{C})$ such that for some invertible matrix
$S$,111Here, $S$ represents a change of basis for the latent variables
resulting in a modification of the system representations from $(A,B,C)$ to
$(S^{-1}AS,S^{-1}B,CS)$–yet describing exactly the same input-output behavior.
See Section 4.4 in (Hespanha, 2018).

$\begin{split}\widehat{A}=S^{-1}AS,\quad\widehat{B}=S^{-1}B,\quad\widehat{C}=CS.\end{split}$

Furthermore, in order to ensure this learning problem is well-posed, we assume
$(A,B)$ is controllable and $(A,C)$ is observable. This is often referred to
as a minimal realization (a state-space description of minimal size that
explains the given input-output data (Schutter, 2000)). The corresponding
system $\mathcal{M}$ is then called a minimal system.

#####  High-dimensional System Identification Problem (HD-SYSID):

Consider learning minimal system
$\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\sigma_{\eta}^{2}I)$ with high-
dimensional observations and low-dimensional latent states and inputs; namely,
$r,m\ll n$. Here, covariances $\Sigma_{w}$, $\sigma_{\eta}^{2}I$ are positive
semi-definite covariances. For $k=1,2$, we choose input trajectories
$\mathcal{U}_{k}=\\{u_{k,t}\\}_{t=0}^{T_{k}-1}$ and get corresponding
observations $\mathcal{Y}_{k}=\\{y_{k,t}\\}_{t=0}^{T_{k}}$. Here, every input
$u_{k,t}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\Sigma_{u})$ is sampled
independently with positive definite covariance $\Sigma_{u}\succ 0$.222Our
approach can be easily extended to a weaker assumption on $\Sigma_{u}$:
$(A,(\Sigma_{w}+B\Sigma_{u}B{}^{\intercal})^{1/2})$ being controllable. Here
we assume $\Sigma_{u}\succ 0$ to keep the results clean and interpretable.
With the two datasets $\mathcal{D}_{1}=\mathcal{Y}_{1}\cup\mathcal{U}_{1}$ and
$\mathcal{D}_{2}=\mathcal{Y}_{2}\cup\mathcal{U}_{2}$, our objective is to
output approximate system matrices
$\left(\widehat{A},\widehat{B},\widehat{C}\right)$ such that with high
probability,

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\leq\epsilon\end{split}$

for some invertible matrix $S$.

In HD-SYSID problem, we specifically consider systems with $r,m\ll n$ and
isotropic observation noise. We select inputs with positive definite
covariance $\Sigma_{u}$ so that the system is fully excited and therefore
learnable. The problem setup can be extended to the case where we have access
to $K\geq 2$ independent data trajectories. Our proposed approach and
theoretical analysis in the following sections can also be readily adapted.

##  3 SYSID with Column Space Projection

Before diving into the details, we first state the necessary assumption and
definition for HD-SYSID.

######  Assumption 3.1.

There exist constants $\psi_{A}\geq 1$ and $\rho_{A}\in(0,1)$, which are
independent of system dimensions, such that

$\begin{split}\left\|A^{i}\right\|\leq\psi_{A}\left(\rho_{A}\right)^{i-1},\quad\forall
i\in\mathbb{N}.\end{split}$

The existence of $\psi_{A}$ and $\rho_{A}$ is a standard assumption (Oymak &
Ozay, 2019) and is guaranteed by Lemma 6 in Talebi et al. (2023) as long as
the spectral radius of $A$ is smaller than $1$. However, the constants may be
dependent on system dimensions in the literature. Here we assume $\psi_{A}$
and $\rho_{A}$ to be independent of system dimensions.

######  Definition 3.2 ( System Identification Oracle (Sys-Oracle)).

Consider a system $\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\Sigma_{\eta})$ with a
minimal realization $(A,B,C)$ and arbitrary dimensions $(r,m,n)$. Assume
$\mathcal{M}$ satisfy Assumption 3.1. Given an input-output trajectory data
sequence $\mathcal{D}=\mathcal{Y}\cup\mathcal{U}$ where
$\mathcal{U}=\\{u_{t}\\}_{t=0}^{T-1}$ with
$u_{t}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\Sigma_{u})$ and
$\mathcal{Y}=\\{y_{t}\\}_{t=0}^{T-1}$ is the corresponding output, the Sys-
Oracle outputs approximation
$\left(\widehat{A},\widehat{B},\widehat{C}\right)$ such that with probability
at least $1-\delta$,

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\\\
{}&\quad\leq\frac{\text{poly}\left(r,n,m\right)}{\sqrt{T}}\cdot\kappa_{2},\end{split}$

for any $T\geq\kappa_{1}\cdot\text{poly}\left(r,n,m\right)$ and some
invertible matrix $S$. Here
$\kappa_{1}=\kappa_{1}\left(\mathcal{M},\mathcal{U},\delta\right),\kappa_{2}=\kappa_{2}\left(\mathcal{M},\mathcal{U},\delta\right)$
are problem-related constants independent of system dimensions modulo
logarithmic factors.

The Sys-Oracle defined above is applicable to all minimal systems with
arbitrary dimensions $(r,n,m)$, which does not necessarily lie in the regime
where $n\gg r,m$. It represents standard system identification algorithms,
including the celebrated Ho-Kalman algorithm (Sarkar et al., 2021; Oymak &
Ozay, 2019). The constraints of such existing algorithms are captured by
Assumption 3.1. These algorithms also require the inputs and noises to be
sampled independently and fully exciting, which is captured by the definition
of system $\mathcal{M}$ and inputs $\mathcal{U}$ in Sys-Oracle. In Appendix D,
we will provide an example of such an oracle for the completeness of this
paper.

###  3.1 The Algorithm

Algorithm 1 Column Space Projection-based SYSID (Col-SYSID)

1: Inputs: Data $\mathcal{Y}_{1},\mathcal{D}_{2}$; Noise covariance
$\Sigma_{\eta}$;Subroutines Col-Approx, Sys-Oracle;

2: Approximate the observer column space:
$\widehat{\Phi}_{C}\leftarrow\texttt{Col-Approx}\left(\mathcal{Y}_{1}\right)$

3: Project dataset
$\mathcal{D}_{2}=\\{y_{2,t}\\}_{t=0}^{T_{2}}\cup\\{u_{2,t}\\}_{t=0}^{T_{2}-1}$
onto the column space:
$\tilde{\mathcal{D}}_{2}\leftarrow\big{\\{}\widehat{\Phi}_{C}{}^{\intercal}y_{2,t}\big{\\}}_{t=0}^{T_{2}}\cup\left\\{u_{2,t}\right\\}_{t=0}^{T_{2}-1}$

4: Identify low-dimensional parameters:
$\widehat{A},\widehat{B},\tilde{C}\leftarrow\acs{oracle}\left(\tilde{\mathcal{D}}_{2}\right)$

5: Recover the high-dimensional observer:
$\widehat{C}\leftarrow\widehat{\Phi}_{C}\tilde{C}$

6: Outputs: $\left(\widehat{A},\widehat{B},\widehat{C}\right)$

Our proposed algorithm consists of two components. Firstly in Line 2, using
the first data trajectory, we approximate the column space of the high-
dimensional observer $C$ and get matrix
$\widehat{\Phi}_{C}\in\mathbb{R}^{n\times\text{rank}(C)}$. The columns of
$\widehat{\Phi}_{C}$ form an orthonormal basis of the space. This
approximation is accomplished by the Column Space Approximation Subroutine
(Col-Approx), which essentially calculates the covariance of the observations
$\\{y_{1,t}\\}_{t=0}^{T_{1}}$ and extracts the eigenspace of the large
eigenvalues. The details of the subroutine is postponed to Section 3.3.

Secondly, in Lines 3-5, we learn the system parameters with the second data
trajectory. With the help of $\widehat{\Phi}_{C}$, we can project the high-
dimensional observations onto lower dimensions, i.e. we let
$\tilde{y}_{2,t}=\widehat{\Phi}_{C}{}^{\intercal}y_{2,t}$. This projected
sequence of observations satisfies the following dynamics

$\begin{split}x_{t+1}={}&Ax_{t}+Bu_{t}+w_{t},\\\
\tilde{y}_{2,t}={}&\widehat{\Phi}_{C}{}^{\intercal}Cx_{t}+\widehat{\Phi}_{C}{}^{\intercal}\eta_{t}.\end{split}$

that is generated by an equivalent system $\widehat{\mathcal{M}}$ denoted by

$\widehat{\mathcal{M}}=(r,\text{rank}(\widehat{\Phi}_{C}),m,A,B,\widehat{\Phi}_{C}{}^{\intercal}C,\Sigma_{w},\sigma_{\eta}^{2}I).$
(2)

######  Remark 3.3 (Why two trajectories?).

. The reason for us to switch to the second trajectory for parameter learning
is to ensure the independence between $\widehat{\Phi}_{C}$ and the second
trajectory, which will ensure that $\widehat{\Phi}_{C}{}^{\intercal}\eta_{t}$
is still independent of other system variables. This is critical for the
application of Sys-Oracle. Note that if the algorithm is equipped with any
Sys-Oracle that handles dependent noise, one can easily simplify Algorithm 1
to a single trajectory. Developing such Sys-Oracle is an interesting future
direction and ideas from (Tian et al., 2023; Simchowitz et al., 2018) may be
related.∎

Lastly, we then feed this low-dimensional dataset to the Sys-Oracle for
learning the corresponding low-dimensional parameters
$\left(\widehat{A},\widehat{B},\tilde{C}\right)\approx\left(A,B,\widehat{\Phi}_{C}{}^{\intercal}C\right)$.
Finally, we recover $C$ from $\tilde{C}$ through $\widehat{\Phi}_{C}$. If
$\widehat{\Phi}_{C}$ approximates the column space of $C$ well, then
$\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}$ should also be close to
$\Phi_{C}\Phi_{C}{}^{\intercal}$, which is the projection onto the column
space of $C$. Therefore, it is intuitive to expect that
$\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}C\approx\Phi_{C}\Phi_{C}{}^{\intercal}C=C$.

###  3.2 Sample Complexity of the Algorithm

We now provide the following sample complexity result for Algorithm 1.

######  Theorem 3.4.

Consider system $\mathcal{M}$ and datasets
$\mathcal{D}_{1}=\mathcal{U}_{1}\cup\mathcal{Y}_{1},\mathcal{D}_{2}=\mathcal{U}_{2}\cup\mathcal{Y}_{2}$
(with lengths $T_{1}$, $T_{2}$ respectively) in HD-SYSID. Suppose
$\mathcal{M}$ satisfies Assumption 3.1. If

$\begin{split}T_{1}\gtrsim\kappa_{3}\cdot n^{2}r^{2},\quad
T_{2}\geq\kappa_{1}\cdot\text{poly}\left(r,m\right),\\\ \end{split}$

then ($\widehat{A},\widehat{B},\widehat{C}$) from Algorithm 1 satisfy the
following for some invertible matrix $S$ with probability at least $1-\delta$

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\\\
{}&\lesssim\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|+\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,m\right)}{T_{2}}}.\end{split}$

Here
$\kappa_{1}=\kappa_{1}(\widehat{\mathcal{M}},\mathcal{U}_{2},\delta),\kappa_{2}=\kappa_{2}(\widehat{\mathcal{M}},\mathcal{U}_{2},\delta),\kappa_{3}=\kappa_{3}\left(\mathcal{M},\mathcal{U}_{[2]},\delta\right)$
and $\kappa_{4}=\kappa_{4}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)$ are
all problem-related constants independent of system dimensions modulo
logarithmic factors with $\widehat{\mathcal{M}}$ defined in Equation 2. Also,
$\kappa_{1}$ and $\kappa_{2}$ are defined in Definition 3.2, while the
definitions of $\kappa_{3},\kappa_{4}$ are summarized in Theorem A.1 in the
appendix.

The above error bound consists of two terms. The first term, also the
dominating term, is the error of observer column space approximation and the
second term is due to learning of the rest of the system. Together, the error
bound directly translates to a sample complexity of
$\tilde{\mathcal{O}}([n+\text{poly}\left(r,m\right)]/\epsilon^{2})$. This
complexity is equivalent to $\tilde{\mathcal{O}}(n/\epsilon^{2})$ in our
setting with $r,m\ll n$.

######  Remark 3.5.

Here the first term indicates an interesting insight: the norm of our
approximated observation matrix may affect the performance of the algorithm.
Therefore, in practice, one might want to choose Sys-Oracle that outputs
$\tilde{C}$ with a reasonable norm (line 4 of Algorithm 1). This ensures
$\left\|\widehat{C}\right\|$ is not too large because
$\left\|\widehat{C}\right\|\leq\left\|\widehat{\Phi}_{C}\right\|\left\|\tilde{C}\right\|=\left\|\tilde{C}\right\|.$
∎

For more concrete results, we instantiate Algorithm 1 with Ho-Kalman algorithm
(Sarkar et al., 2021; Oymak & Ozay, 2019), Algorithm 4 in Appendix D as Sys-
Oracle. Then the above theorem leads to the following corollary.

######  Corollary 3.6.

Consider the special minimal system
$\mathcal{M}=(r,n,m,A,B,C,\sigma_{w}^{2}I,\sigma_{\eta}^{2}I)$ and datasets
$\mathcal{D}_{1}=\mathcal{U}_{1}\cup\mathcal{Y}_{1},\mathcal{D}_{2}=\mathcal{U}_{2}\cup\mathcal{Y}_{2}$
(with length $T_{1},T_{2}$ respectively), where the inputs are sampled
independently from $\mathcal{N}(0,\sigma_{u}^{2}I)$. Suppose $\mathcal{M}$
satisfies Assumption 3.1. If

$\begin{split}&T_{1}\gtrsim\tilde{\kappa}_{3}\cdot n^{2}r^{2},\quad
T_{2}\geq\tilde{\kappa}_{1}\cdot r^{3}(r+m),\end{split}$

then $(\widehat{A},\widehat{B},\widehat{C})$ from Algorithm 1 satisfy the
following for some invertible matrix $S$ with probability at least $1-\delta$

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\\\
{}&\lesssim\tilde{\kappa}_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|+\tilde{\kappa}_{2}\cdot\sqrt{\frac{r^{5}(r+m)}{T_{2}}}.\end{split}$

Here
$\tilde{\kappa}_{1},\tilde{\kappa}_{2},\tilde{\kappa}_{3},\tilde{\kappa}_{4}$
are all problem-related constants only dependent of system $\mathcal{M}$ and
inputs $\mathcal{U}_{1}\cup\mathcal{U}_{2}$ (as compared to $\kappa_{1}$ and
$\kappa_{2}$ in Theorem 3.4 that depend on $\widehat{\mathcal{M}}$). They are
also independent of system dimensions modulo log factors. Detailed definitions
of the constants are listed in Corollary D.2.

Based on the above corollary, we now compare the error of our algorithm
instantiated with the example oracle, denoted by $\Delta_{1}$, with the error
of the example oracle when it is directly applied to the dataset, denoted by
$\Delta_{2}$. For a fair comparison, we consider $\Delta_{2}$ as the error
with datasets generated by inputs
$\mathcal{U}=\\{u_{t}\\}_{t=0}^{T_{1}+T_{2}-1}$ of length $T_{1}+T_{2}$. Here
the inputs are sampled independently from $\mathcal{N}(0,\sigma_{u}^{2}I)$.

The error of our algorithm is upper bounded by

$\Delta_{1}\leq\tilde{\mathcal{O}}\left(\sqrt{\frac{n}{T_{1}}}+\sqrt{\frac{r^{5}(r+m)}{T_{2}}}\right).$

The above error directly translates to the sample complexity of
$\tilde{\mathcal{O}}([n+r^{5}(r+m)]/\epsilon^{2})=\tilde{\mathcal{O}}(n/\epsilon^{2})$
when $r,m\ll n$. On the other hand, from Corollary D.1

$\Delta_{2}\leq\tilde{\mathcal{O}}\left(\sqrt{\frac{r^{5}(r+n+m)}{T_{1}+T_{2}}}\right).$

And this gives a sample complexity of order
$\tilde{\mathcal{O}}([nr^{5}+r^{5}(r+m)]/\epsilon^{2})=\tilde{\mathcal{O}}(nr^{5}/\epsilon^{2})$
when $r,m\ll n$. Therefore, in the regime where $r,m\ll n$, it is clear that
our algorithm enjoys a better sample complexity as compared to directly
applying the example oracle.

###  3.3 Proof Ideas of Theorem 3.4

Although conceptually simple, justifying this tight sample complexity result
is hard. The difficulty lies in ensuring the accuracy of the observer column
space approximation $\widehat{\Phi}_{C}$, which is learned mainly by PCA on
the given data (detailed in the following algorithm).

Algorithm 2 Col-Approx

1: Input: Observations $\mathcal{Y}=\\{y_{t}\\}_{t=0}^{T}$;

2: Calculating data covariance
$\Sigma_{y}\leftarrow\sum_{t=0}^{T}y_{t}y_{t}{}^{\intercal}$

3: Estimating observer rank
$\widehat{r}_{c}\leftarrow\mathop{\arg\max}_{i}\left(\sigma_{i}(\Sigma_{y})-\sigma_{i+1}(\Sigma_{y})>T^{3/4}\right)$

4: Estimating column space  $\widehat{\Phi}_{C}\leftarrow\text{ first
$~{}\widehat{r}_{c}~{}$ eigenvectors of }~{}\Sigma_{y}$

5: Output: $\widehat{\Phi}_{C}$

Existing analysis of PCA (Candès et al., 2011; Vaswani et al., 2018; Chen &
Storey, 2015), matrix factorization (Gribonval et al., 2015) and other
subspace learning techniques (Tripuraneni et al., 2021; Zhang et al., 2023) do
not apply to our setting. This is because these methods assume i.i.d. data
samples, while our dataset $\mathcal{D}_{1}$ consists of correlated data
points generated by a dynamical system. We need to apply martingale tools
(Abbasi-yadkori et al., 2011; Sarkar & Rakhlin, 2019) for our analysis.

Moreover, the noises in the dynamical systems accumulate across time steps,
leading to large noise terms in the analysis. Directly applying the
aforementioned tools leads to loose sample complexity bounds. Since we are
learning a subspace embedded in $\mathbb{R}^{n}$, it is reasonable to expect a
sample complexity of $\mathcal{O}(n/\epsilon^{2})$ for an $\epsilon$-good
approximation of the subspace 333Other intuitions for this complexity comes
from standard covariance concentration results for i.i.d. data (Corollary 2.1
of (Rudelson, 1999)). In these results, $\mathcal{O}(n/\epsilon^{2})$ samples
are required to learn an $\epsilon$-good covariance, which directly translates
to a good column space estimation of the covariance. . We achieve this result
by developing our own subspace perturbation lemma specifically tailored to
this setting. The above ideas on analyzing Col-Approx translate to the
following lemma.

######  Lemma 3.7.

Consider system $\mathcal{M}$ and dataset
$\mathcal{D}_{1}=\mathcal{U}_{1}\cup\mathcal{Y}_{1}$ in HD-SYSID. Suppose
$\mathcal{M}$ satisfies 3.1. If

$\begin{split}T_{1}\gtrsim\kappa_{5}\cdot n^{2}r^{2},\end{split}$

then $\widehat{\Phi}_{C}=\texttt{Col-Approx}(\mathcal{Y}_{1},\Sigma_{\eta})$
satisfies the following with probability at least $1-\delta$

$\begin{split}\text{rank}(\widehat{\Phi}_{C})=\text{rank}(C),\quad\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\lesssim\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}.\end{split}$

Here $\kappa_{4}=\kappa_{4}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)$
and $\kappa_{5}=\kappa_{5}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)$ are
both problem-related constants independent of system dimensions modulo
logarithmic factors (details in Lemma A.2).

Recall $\Phi_{C}$ denotes the orthonormal matrix whose columns form a basis of
$\text{col}(C)$ and $\widehat{\Phi}_{C}^{\perp}$ denote the matrix such that
$\begin{bmatrix}\widehat{\Phi}_{C}^{\perp}&\widehat{\Phi}_{C}\end{bmatrix}$ is
unitary.

###### Proof Sketch of Lemma 3.7.

$\widehat{\Phi}_{C}$ is constructed from the eigenvectors of matrix
$\Sigma_{y}=\sum_{t=0}^{T_{1}}y_{1,t}y_{1,t}{}^{\intercal}$. We decompose it
as follows

$\begin{split}{}&\sum_{t=0}^{T_{1}}y_{1,t}y_{1,t}{}^{\intercal}=\underbrace{C\left(\sum_{t=0}^{T_{1}}x_{1,t}x_{1,t}{}^{\intercal}\right)C{}^{\intercal}}_{\Sigma_{C}}+\sigma_{\eta}^{2}(T_{1}+1)I\\\
+{}&\underbrace{\sum_{t=0}^{T_{1}}\left(\eta_{1,t}\eta_{1,t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)}_{\Delta_{2}}+\underbrace{\sum_{t=0}^{T_{1}}\left(Cx_{1,t}\eta_{1,t}{}^{\intercal}+\eta_{1,t}x_{1,t}{}^{\intercal}C{}^{\intercal}\right)}_{\Delta_{1}+\Delta_{1}{}^{\intercal}}.\end{split}$

Here $\Sigma_{C}$ contains the information on the observer column space
$\Phi_{C}$, while $\Delta_{1}$ and $\Delta_{2}$ are noise terms. The rest of
the proof is decomposed into three steps. In the first step, we upper bound
the noise terms $\Delta_{1}$, $\Delta_{2}$ and lower bound the latent state
covariance $\sum_{t=0}^{T_{1}}x_{1,t}x_{1,t}$. In the second step, we show
that the eigenvectors in the column space of $C$ has eigenvalues much larger
then the other eigenvalues. This gap eigenvalue gap enables us to identify the
dimension of the column space of $C$. Finally, in the third step, we apply our
own subspace perturbation result to bound the accuracy of the approximated
column space.

Step 1. With system $(A,B)$ being controllable, the inputs injected fully
excite the latent dynamics, leading to the following lower bound on
$\sum_{t=0}^{T_{1}}x_{1,t}x_{1,t}$:

$\begin{split}\sum_{t=0}^{T_{1}}x_{1,t}x_{1,t}\succ\tilde{\mathcal{O}}(T_{1})I.\end{split}$

The upper bound on $\Delta_{2}$, i.e.
$\left\|\Delta_{2}\right\|\leq\tilde{\mathcal{O}}(\sqrt{T_{1}})$, follows from
standard guassian concentration arguments. To upper bound $\Delta_{1}$, we
apply the martingale tools from previous works (Abbasi-yadkori et al., 2011;
Sarkar & Rakhlin, 2019) and get the following relative error bounds

$\begin{split}\left\|(\Sigma_{C}+T_{1}I)^{-\frac{1}{2}}\Delta_{1}\right\|\leq\tilde{\mathcal{O}}(\sqrt{T_{1}}).\end{split}$
(3)

Step 2. Since the perturbations are small, standard eigenvalue perturbation
bounds gives

$\begin{split}\sigma_{\text{rank}(C)}(\Sigma_{y})\geq\tilde{\mathcal{O}}(T_{1}),\quad\sigma_{\text{rank}(C)+1}(\Sigma_{y})\leq\tilde{\mathcal{O}}(\sqrt{T_{1}}).\end{split}$

This eigenvalue gap as large as $\tilde{\mathcal{O}}(T_{1})$ makes it easy to
determine the rank of $C$. Therefore, with high probability,
$\text{rank}(\widehat{\Phi}_{C})=\text{rank}(C)$.

Step 3. Finally, we bound the accuracy of the eigenspace of
$\Sigma_{y}=\Sigma_{C}+\Delta_{2}+(\Delta_{1}+\Delta_{1}{}^{\intercal})$. Here
$\Delta_{1}$ is the so-called relative perturbation, because we bound its norm
by comparing it with the information $(\Sigma_{C}+T_{1}I)^{\frac{1}{2}}$ (in
Equation 3). Since $\left\|(\Sigma_{C}+T_{1}I)^{\frac{1}{2}}\right\|$ can be
large, directly applying standard subspace perturbation results leads to sub-
optimal error bound. For this, we develop a specific relative subspace
perturbation bound, adapted from classic eigenspace perturbation results. It
ensures both $\Delta_{2}$ and the relative noise $\Delta_{1}$ will not perturb
the eigenspace of $\Sigma_{y}$ too much. This finishes the proof. ∎

With the above lemma, the proof of Theorem 3.4 can then be decomposed into two
steps. In the first step, we guarantee the feasibiliy to learn the low-
dimensional system parameters $(A,B,\widehat{\Phi}_{C}{}^{\intercal}C)$ from
the projected dataset (Line 3 in Algorithm 1). Since $\widehat{\Phi}_{C}$ is
accurate enough, the projected dataset almost preserves all the information.
This ensures that the equivalent system generating the projected dataset
($\widehat{\mathcal{M}}$ in Equation 2) is observable and controllable.
Therefore, Sys-Oracle outputs accurate estimation of the low-dimensional
system parameters. In the second step, we analyze the errors of the recovered
high-dimensional parameters
$(\widehat{A},\widehat{B},\widehat{\Phi}_{C}\tilde{C})$, where $\tilde{C}$ is
our estimation of $\widehat{\Phi}_{C}{}^{\intercal}C$. The errors on
$\widehat{A}$ and $\widehat{B}$ are automatically bounded by the definition of
Sys-Oracle. The error on the recovered high-dimensional observer
$\widehat{C}=\widehat{\Phi}_{C}\tilde{C}$ can be decomposed as follows

$\begin{split}{}\left\|CS-\widehat{C}\right\|\leq{}&\left\|CS-\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS\right\|+\left\|\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS-\widehat{\Phi}_{C}\tilde{C}\right\|\\\
={}&\left\|\widehat{\Phi}_{C}^{\perp}\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}CS\right\|+\left\|\widehat{\Phi}_{C}{}^{\intercal}CS-\tilde{C}\right\|.\end{split}$

Here the first term is the column space approximation error and the second
term is the error learning the low-dimensional observer
$\tilde{C}\approx\widehat{\Phi}_{C}{}^{\intercal}C$. Combining the above
inequality with Lemma 3.7 and Definition 3.2, along with some algebraic
arguments, proves the theorem.

##  4 Lower Bounds

Now that Algorithm 1 accomplishes HD-SYSID with sample complexity
$\tilde{\mathcal{O}}(n/\epsilon^{2})$, one natural question is: _Is this
linear dependency on observer dimension $n$ unavoidable?_ In this section, we
provide the following theorem to show that this dependency is necessary.

While we present this lower bound for learning LTI systems with a single
input-output data trajectory, it can be extended to multiple trajectories. The
more general version is deferred to the appendix (Theorem B.2).

######  Theorem 4.1.

Suppose $n\geq 2$, and choose positive scalars $\delta\leq\frac{1}{2}$ and
$\epsilon\leq 0.6$. Consider the class of minimal systems
$\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\Sigma_{\eta})$ with different
$A$,$B$,$C$ matrices. All parameters except $A,B,C$ are fixed and known.
Moreover, $r<n$ and $\Sigma_{\eta}$ is positive definite. Let
$\mathcal{D}=\\{y_{t}\\}_{t=0}^{T}\cup\\{u_{t}\\}_{t=0}^{T-1}\cup\text{\\{all
known parameters\\}}$ denote the associated single trajectory dataset. Here
the input $u_{t}$ satisfy: 1). $u_{t}$ is independently sampled; 2).
$\mathbb{E}(u_{t})=0$. Consider any estimator $\widehat{f}$ mapping
$\mathcal{D}$ to
$(\widehat{A}(\mathcal{D}),\widehat{B}(\mathcal{D}),\widehat{C}(\mathcal{D}))\in\mathbb{R}^{r\times
r}\times\mathbb{R}^{r\times m}\times\mathbb{R}^{n\times r}$. If

$T<\frac{\phi_{\eta}(1-2\delta)\log
1.4}{50\left(\psi_{w}+\psi_{u}\right)}\cdot\frac{n+\log\frac{1}{2\delta}}{\epsilon^{2}},$

there exists a system
$\mathcal{M}_{0}=(r,n,m,A_{0},B_{0},C_{0},\Sigma_{w},\Sigma_{\eta})$ with
dataset $\mathcal{D}$ such that

$\begin{split}\mathbb{P}{}&\left\\{\left\|C_{0}B_{0}-\widehat{C}\left(\mathcal{D}\right)\widehat{B}\left(\mathcal{D}\right)\right\|\geq\epsilon\right\\}\geq\delta.\end{split}$

Here $\mathbb{P}$ denotes the distribution of $\mathcal{D}$ generated by
system $\mathcal{M}_{0}$. Related constants are defined as follows

$\begin{split}\phi_{\eta}&\coloneqq\sigma_{\min}(\Sigma_{\eta}),\quad\psi_{w}\coloneqq\sigma_{1}(\Sigma_{w}),\\\
\quad\psi_{u}&\coloneqq\max_{t\in[0,T-1]}\sigma_{1}\left(\mathbb{E}(u_{t}u_{t}{}^{\intercal})\right).\end{split}$

The above theorem indicates that with less than $\mathcal{O}(n/\epsilon^{2})$
samples, any estimator has a constant probability to fail learning the product
$CB$. This error lower bound on
$\widehat{C}(\mathcal{D})\widehat{B}(\mathcal{D})$ can be translated to
corresponding errors on

$\begin{split}\max\left\\{\left\|S^{-1}B-\widehat{B}\left(\mathcal{D}\right)\right\|,\left\|CS-\widehat{C}\left(\mathcal{D}\right)\right\|\right\\}\geq\epsilon,\end{split}$

under mild conditions. This indicates that the estimator fails to learn either
$B$ or $C$ well. (The proof details are provided in Corollary B.3). Another
perspective to understand this error on
$\widehat{C}(\mathcal{D})\widehat{B}(\mathcal{D})$ is as follows. In the
celebrated Ho-Kalman algorithm for HD-SYSID, the estimation error the Hankel
operator is lower bounded by the estimation error of
$\widehat{C}(\mathcal{D})\widehat{B}(\mathcal{D})$. Therefore, our lower bound
indicates the failure of estimating the Hankel operator, which prevents the
algorithm from giving satisfying outputs.

The above results then give a $\mathcal{O}(n/\epsilon^{2})$ sample complexity
lower bound on Sys-Oracle. Therefore, the sample complexity of Algorithm 1,
i.e. $\tilde{\mathcal{O}}(n/\epsilon^{2})$, is optimal up to logarithmic and
problem-related constants.

We defer the proof of Theorem 4.1 to Appendix B. But we would like to briefly
discuss the reason behind this unavoidable dependency on $n$. Therein, we
construct a set of system instances with different observation matrices such
that their column spaces are “relatively close” but “quantifiably distinct”.
In fact, we take advantage of the fact that in high dimensional space
$\mathbb{R}^{n}$, there are too many distinct subspaces that are close to each
other. And therefore, any estimator needs correspondingly many samples, i.e.
$\mathcal{O}(n/\epsilon^{2})$, to distinguish the true column space from many
other candidates. This difficulty is indeed due to the high dimensional
observation noise.

In the following section, we introduce a more general problem pertinent to
high-dimensional “meta-datasets.” While this lower bound still has
implications even for this general problem, we show how to confine its implied
difficulty on the total length of the meta-dataset and not on its personalized
portions.

##  5 Meta SYSID

Thus far, we have introduced Algorithm 1 for HD-SYSID that provides near-
optimal dependency on system dimensions utilizing the idea of low-dimensional
embedding. Now we will investigate how this algorithmic idea facilitates
learning from so-called meta-datasets (or metadata)—multiple datasets from
different systems sharing the same observer. As shown in (Hajnal et al., 2023;
Xia et al., 2021; Marks & Goard, 2021; Gallego et al., 2020), this is a real-
world setting of significant importance.

Following our analysis in Section 3.3, we observe that the
$\tilde{\mathcal{O}}(n/\epsilon^{2})$ sample complexity arises from the error
of column space $\widehat{\Phi}_{C}$. Modulo this error, the sample complexity
reduces vastly to
$\tilde{\mathcal{O}}(\text{poly}\left(r,m\right)/\epsilon^{2})$. This
motivates us to again separate the learning of the observer column space
$\Phi_{C}$ from the learning of personalized system parameters. Utilizing
multiple datasets from different systems, we can collectively learn the
observer column space much more accurately. Subsequently, equipped with this
accurate column space approximation, there is hope to learn every single
system with much fewer samples. To formalize the above idea, we introduce the
following “ High-dimensional Meta System Identification Problem (Meta-SYSID).”

High-dimensional Meta System Identification Problem (Meta-SYSID): Consider the
identification of $K$ minimal systems
$\mathcal{M}_{k}=(r,n,m,A_{k},B_{k},C,\Sigma_{w,k},\sigma_{\eta,k}^{2}I),k\in[K]$
with the same dimensions ($r,m\ll n$) and observation matrix $C$. Here
covariances $\\{\Sigma_{w,k},\sigma_{\eta,k}^{2}I\\}_{k\in[K]}$ are positive
semi-definite matrices. For every system $\mathcal{M}_{k}$, we choose an input
sequence $\mathcal{U}_{k}=\\{u_{k,t}\\}_{t=0}^{T_{k}-1}$ and get observations
$\mathcal{Y}_{k}=\\{y_{k,t}\\}_{t=0}^{T_{k}}$. Here every input
$u_{k,t}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\Sigma_{u,k})$ is sampled
independently from both the system variables and other inputs with positive
definite covariance $\Sigma_{u,k}$. This single trajectory dataset is denoted
by $\mathcal{D}_{k}=\mathcal{Y}_{k}\cup\mathcal{U}_{k}$. With the $K$ datasets
$\bigcup_{k\in[K]}\mathcal{D}_{k}$, our objective is to learn all system
parameters well. Namely, we aim to identify
$\\{\widehat{A}_{k},\widehat{B}_{k},\widehat{C}_{k}\\}_{k=1}^{K}$ such that,
for every $k\in[K]$, the following holds for some invertible matrix $S_{k}$
with high probability

$\begin{split}\max\left\\{\left\|S_{k}^{-1}A_{k}S_{k}-\widehat{A}_{k}\right\|,\left\|S_{k}^{-1}B_{k}-\widehat{B}_{k}\right\|,\right.\\\
\left.\left\|CS_{k}-\widehat{C}_{k}\right\|\right\\}&\leq\epsilon.\end{split}$

###  5.1 The Meta-Learning Algorithm and Its Sample Complexity

Algorithm 3 Meta Space Projection-based SYSID (Meta-SYSID)

1: Inputs: Meta Datasets
$\mathcal{D}_{[K]}=\mathcal{Y}_{[K]}\cup\mathcal{U}_{[K]}$;Subroutines Col-
Approx, Sys-Oracle;

2: for $k\in[K]$ do

3: Leave one out and approximate observer column space
$\widehat{\Phi}_{C,k}\leftarrow\texttt{Col-
Approx}\left(\mathcal{Y}_{-k}\right)$

4: Project dataset
$\mathcal{D}_{k}=\\{y_{k,t}\\}_{t=0}^{T_{k}}\cup\\{u_{k,t}\\}_{t=0}^{T_{k}-1}$
onto the column space:
$\tilde{\mathcal{D}}_{k}\leftarrow\big{\\{}\widehat{\Phi}_{C,k}{}^{\intercal}y_{k,t}\big{\\}}_{t=0}^{T_{k}}\cup\left\\{u_{k,t}\right\\}_{t=0}^{T_{k}-1}$

5: Identify low-dimensional parameters:
$\widehat{A}_{k},\widehat{B}_{k},\tilde{C}_{k}\leftarrow\acs{oracle}\left(\tilde{\mathcal{D}}_{k}\right)$

6: Recover the high-dimensional observer:
$\widehat{C}_{k}\leftarrow\widehat{\Phi}_{C,k}\tilde{C}_{k}$

7: end for

8: Outputs: $\\{\widehat{A}_{k},\widehat{B}_{k},\widehat{C}_{k}\\}_{k\in[K]}$

Algorithm 3 follows similar steps as Algorithm 1. For every system $k$, we
approximate the column space, learn the low-dimensional system parameters, and
finally project them back into the high-dimension. And similarly, to ensure
the independence between $\widehat{\Phi}_{C,k}$ and $\tilde{\mathcal{D}}_{k}$,
the dataset $\tilde{\mathcal{D}}_{k}$ itself is left out in the first step.
This independence is critical for the application of the Sys-Oracle
subroutine.

Figure 1: Left: Error for Column Space Approximation (Col-Approx, Algorithm
2). Center: Error for HD-SYSID (Algorithm 1) and Standard Ho-Kalman with
$n=160$. Right: Error for Meta-SYSID (Algorithm 3) and Standard Ho-Kalman for
$n\in[80,320]$.

Now we provide the theoretical guarantee for the above algorithm.

######  Theorem 5.1.

Consider systems $\mathcal{M}_{[K]}$ and datasets
$\mathcal{D}_{[K]}=\mathcal{U}_{[K]}\cup\mathcal{Y}_{[K]}$ in Meta-SYSID.
Suppose $\mathcal{M}_{[K]}$ satisfy Assumption 3.1. Then for any fixed system
$k_{0}\in[K]$, if $\\{T_{k}\\}_{k\in[K]}$ and $T_{-k_{0}}\coloneqq\sum_{k\neq
k_{0}}T_{k}$ satisfy

$\begin{split}&T_{k_{0}}\geq\kappa_{1}\cdot\text{poly}\left(r,m\right),\quad
T_{-k_{0}}\gtrsim\kappa_{3}\cdot n^{2}r^{2}\log^{8}K,\end{split}$

then $(\widehat{A}_{k_{0}},\widehat{B}_{k_{0}},\widehat{C}_{k_{0}})$ from
Algorithm 3 satisfy the following for some invertible matrix $S$ with
probability at least $1-\delta$

$\begin{split}{}&\max\left\\{\left\|S^{-1}A_{k_{0}}S-\widehat{A}_{k_{0}}\right\|,\left\|S^{-1}B_{k_{0}}-\widehat{B}_{k_{0}}\right\|,\left\|C_{k_{0}}S-\widehat{C}_{k_{0}}\right\|\right\\}\\\
{}&\lesssim\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,m\right)}{T_{k_{0}}}}+\kappa_{4}\cdot\sqrt{\frac{n}{T_{-k_{0}}}}\left\|\widehat{C}_{k_{0}}\right\|.\end{split}$

Here
$\kappa_{1}=\kappa_{1}(\widehat{\mathcal{M}}_{k_{0}},\mathcal{U}_{k_{0}},\delta),\kappa_{2}=\kappa_{2}(\widehat{\mathcal{M}}_{k_{0}},\mathcal{U}_{k_{0}},\delta),\kappa_{3}=\kappa_{3}\left(\mathcal{M}_{[K]},\mathcal{U}_{[K]},\delta\right),\kappa_{4}=\kappa_{4}\left(\mathcal{M}_{-k_{0}},\mathcal{U}_{-k_{0}},\delta\right)$
are all problem-related constants independent of system dimensions modulo
logarithmic factors. $\kappa_{1},\kappa_{2}$ are defined in Definition 3.2,
while the definitions of $\kappa_{3},\kappa_{4}$ are summarized in Theorem
C.1. $\widehat{\mathcal{M}}_{k_{0}}$ is define as follows

$\begin{split}\widehat{\mathcal{M}}_{k_{0}}=({}&r,\text{rank}(\widehat{\Phi}_{C,k_{0}}),m,\\\
{}&A_{k_{0}},B_{k_{0}},\widehat{\Phi}_{C,k_{0}}{}^{\intercal}C,\Sigma_{w},\sigma_{\eta,k}^{2}I).\end{split}$

In the above result, we only require
$T_{k_{0}}=\tilde{\mathcal{O}}(\text{poly}\left(r,m\right)/\epsilon^{2})$ for
an $\epsilon$ accurate approximation when
$T_{-k_{0}}=\tilde{\mathcal{O}}(n/\epsilon^{2})$. Namely, as long as we have
enough metadata, the number of samples needed from single systems are vastly
reduced and is independent of the observer dimension $n$. In real-world
applications when we have large, i.e. $\tilde{\mathcal{O}}(n/\epsilon^{2})$,
dataset from a single system, the above result significantly helps for few-
shot learning of other similar systems. Also, whenever we have access to
numerous, i.e. $\tilde{\mathcal{O}}(n)$ similar systems, the samples required
from every system is also independent of $n$.

##  6 Simulations

We simulate our algorithm for a set of simple systems where $r=m=1$,
$\Sigma_{w}=0$ and $\sigma_{\eta}=1$. We use $A=0.9,B=1$ and randomly sample
$C$ with orthonormal columns. We choose inputs with covariance
$\Sigma_{u}=0.1I$. The choice of $n$ will be clear in the context. The results
are reported in Figure 1.

We first simulate Col-Approx (Algorithm 2) with $n=40,80,160,320$ and a single
trajectory data of length $T=10000$ separately (Figure 1, left). For isotropic
observation noises as in HD-SYSID, i.e. $\Sigma_{\eta}=\sigma_{\eta}^{2}I$,
the approximated column space converges to the true column space as is shown
by the solid lines. Moreover, the algorithm is simulated for non-isotropic
observation noise where the covariance is set to be
$\Phi_{C}\Phi_{C}{}^{\intercal}$. As shown in by the dotted curves, the
algorithm fail to converge. This is because non-isotropic noise perturbs the
eigenspace of the observation covariance $\Sigma_{y}$ too much so that the
information on $\Phi_{C}$ is drowned.

Next, we simulate Algorithm 1 (Our Algorithm) with $T=10000$ and $n=160$
(Figure 1, center). We plot the error, i.e. $\|\widehat{C}\widehat{B}-CB\|$,
during the learning process. As a comparison, error for the standard Ho-Kalman
(Oymak & Ozay, 2019) (Naive Algorithm) is also included. It is clear that our
Algorithm 3 outperforms the standard Ho-Kalman.

Finally, we consider the Meta-SYSID setting for $n\in[80,320]$ with a meta-
dataset from $\lfloor n/40\rfloor$ randomly sampled systems. Without loss of
generality, the error $\|\widehat{C}\widehat{B}-CB\|$ for learning the first
system is plotted. For Algorithm 3 (Our Algorithm), it is clear that the
learning error doesn’t grow as the observation dimension increases. This is
because our algorithm utilizes information from the meta-dataset to learn the
observer column space, which constitutes the major challenge of system
learning. However, the standard Ho-Kalman algorithm learns every system
independently, and therefore the error is influenced by the observation
dimension.

##  7 Conclusions and Future Directions

In conclusion, our focus has been on learning a linear time-invariant (LTI)
model characterized by low-dimensional latent variables and high-dimensional
observations. The introduced Col-SYSID Algorithm serves as a solution with a
commendable complexity of $\tilde{\mathcal{O}}(n/\epsilon^{2})$. Our analysis
also delves into the fundamental limitations of this problem, establishing a
sample complexity lower bound that essentially underscores the optimality of
our proposed algorithm. Extending the scope of our results, we address a meta-
learning setting where datasets from multiple analogous systems are available.
This leads to the Meta-SYSID algorithm, an end-to-end framework adept at
managing the meta-dataset and effectively learning all included systems.

While our current work lays a solid foundation, future directions could
explore extensions to non-linear settings, or investigate adaptive approaches
to handle varying or completely unknown observation noise. Additionally,
incorporating real-world applications and practical considerations could
further enrich the utility of our results.

## Impact Statement

The goal of this work is to advance theory in high-dimension machine learning,
specifically in latent space space learning. There might be minor potential
societal consequences, but none of those are considered as necessary to be
specifically highlighted here.

## Acknowledgments

This work is supported by NSF AI institute 2112085, NSF ECCS 2328241, and NIH
R01LM014465.

## References

  * Abbasi-yadkori et al. (2011) Abbasi-yadkori, Y., Pál, D., and Szepesvári, C.  Improved Algorithms for Linear Stochastic Bandits.  In _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.  URL https://papers.nips.cc/paper_files/paper/2011/hash/e1d5be1c7f2f456670de3d53c7b54f4a-Abstract.html.
  * Anderson et al. (2022) Anderson, B. D. O., Deistler, M., and Lippi, M.  Linear System Challenges of Dynamic Factor Models.  _Econometrics_ , 10(4):35, December 2022.  ISSN 2225-1146.  doi: 10.3390/econometrics10040035.  URL https://www.mdpi.com/2225-1146/10/4/35.  Number: 4 Publisher: Multidisciplinary Digital Publishing Institute.
  * Bakshi et al. (2023) Bakshi, A., Liu, A., Moitra, A., and Yau, M.  A New Approach to Learning Linear Dynamical Systems.  In _Proceedings of the 55th Annual ACM Symposium on Theory of Computing_ , STOC 2023, pp. 335–348, New York, NY, USA, June 2023. Association for Computing Machinery.  ISBN 978-1-4503-9913-5.  doi: 10.1145/3564246.3585247.  URL https://dl.acm.org/doi/10.1145/3564246.3585247.
  * Balzano et al. (2018) Balzano, L., Chi, Y., and Lu, Y. M.  Streaming PCA and Subspace Tracking: The Missing Data Case.  _Proceedings of the IEEE_ , 106(8):1293–1310, August 2018.  ISSN 1558-2256.  doi: 10.1109/JPROC.2018.2847041.  URL https://ieeexplore.ieee.org/document/8417980.  Conference Name: Proceedings of the IEEE.
  * Blanchard et al. (2007) Blanchard, G., Bousquet, O., and Zwald, L.  Statistical properties of kernel principal component analysis.  _Machine Learning_ , 66(2):259–294, March 2007.  ISSN 1573-0565.  doi: 10.1007/s10994-006-6895-9.  URL https://doi.org/10.1007/s10994-006-6895-9.
  * Boucheron et al. (2013) Boucheron, S., Lugosi, G., and Massart, a. P.  _Concentration Inequalities: A Nonasymptotic Theory of Independence_.  Oxford University Press, Oxford, New York, March 2013.  ISBN 978-0-19-953525-5.
  * Breitung & Eickmeier (2006) Breitung, J. and Eickmeier, S.  Dynamic factor models.  _Allgemeines Statistisches Archiv_ , 90(1):27–42, March 2006.  ISSN 1614-0176.  doi: 10.1007/s10182-006-0219-z.  URL https://doi.org/10.1007/s10182-006-0219-z.
  * Bui-Thanh et al. (2008) Bui-Thanh, T., Willcox, K., and Ghattas, O.  Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space.  _SIAM Journal on Scientific Computing_ , 30(6):3270–3288, January 2008.  ISSN 1064-8275.  doi: 10.1137/070694855.  URL https://epubs.siam.org/doi/abs/10.1137/070694855.  Publisher: Society for Industrial and Applied Mathematics.
  * Candès et al. (2011) Candès, E. J., Li, X., Ma, Y., and Wright, J.  Robust principal component analysis?  _Journal of the ACM_ , 58(3):11:1–11:37, June 2011.  ISSN 0004-5411.  doi: 10.1145/1970392.1970395.  URL https://dl.acm.org/doi/10.1145/1970392.1970395.
  * Chen & Storey (2015) Chen, X. and Storey, J. D.  Consistent Estimation of Low-Dimensional Latent Structure in High-Dimensional Data, October 2015.  URL http://arxiv.org/abs/1510.03497.  arXiv:1510.03497 [stat].
  * Churchland et al. (2012) Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu, S. I., and Shenoy, K. V.  Neural population dynamics during reaching.  _Nature_ , 487(7405):51–56, July 2012.  ISSN 1476-4687.  doi: 10.1038/nature11129.  URL https://www.nature.com/articles/nature11129.  Number: 7405 Publisher: Nature Publishing Group.
  * Davis & Kahan (1970) Davis, C. and Kahan, W. M.  The Rotation of Eigenvectors by a Perturbation. III.  _SIAM Journal on Numerical Analysis_ , 7(1):1–46, 1970.  ISSN 0036-1429.  URL https://www.jstor.org/stable/2949580.  Publisher: Society for Industrial and Applied Mathematics.
  * Deng et al. (2020) Deng, H., Chen, W., Shen, Q., Ma, A. J., Yuen, P. C., and Feng, G.  Invariant subspace learning for time series data based on dynamic time warping distance.  _Pattern Recognition_ , 102:107210, June 2020.  ISSN 0031-3203.  doi: 10.1016/j.patcog.2020.107210.  URL https://www.sciencedirect.com/science/article/pii/S0031320320300169.
  * Djehiche & Mazhar (2021) Djehiche, B. and Mazhar, O.  Non asymptotic estimation lower bounds for LTI state space models with Cram\’er-Rao and van Trees, September 2021.  URL http://arxiv.org/abs/2109.08582.  arXiv:2109.08582 [math, stat].
  * Djehiche & Mazhar (2022) Djehiche, B. and Mazhar, O.  Efficient learning of hidden state LTI state space models of unknown order, February 2022.  URL http://arxiv.org/abs/2202.01625.  arXiv:2202.01625 [math, stat].
  * Dong et al. (2022) Dong, Y., Qin, S. J., and Boyd, S. P.  Extracting a low-dimensional predictable time series.  _Optimization and Engineering_ , 23(2):1189–1214, June 2022.  ISSN 1389-4420, 1573-2924.  doi: 10.1007/s11081-021-09643-x.  URL https://link.springer.com/10.1007/s11081-021-09643-x.
  * Fattahi (2021) Fattahi, S.  Learning Partially Observed Linear Dynamical Systems from Logarithmic Number of Samples.  In _Proceedings of the 3rd Conference on Learning for Dynamics and Control_ , pp. 60–72. PMLR, May 2021.  URL https://proceedings.mlr.press/v144/fattahi21a.html.  ISSN: 2640-3498.
  * Gallego et al. (2020) Gallego, J. A., Perich, M. G., Chowdhury, R. H., Solla, S. A., and Miller, L. E.  Long-term stability of cortical population dynamics underlying consistent behavior.  _Nature Neuroscience_ , 23(2):260–270, February 2020.  ISSN 1546-1726.  doi: 10.1038/s41593-019-0555-4.  URL https://www.nature.com/articles/s41593-019-0555-4.  Number: 2 Publisher: Nature Publishing Group.
  * Gribonval et al. (2015) Gribonval, R., Jenatton, R., Bach, F., Kleinsteuber, M., and Seibert, M.  Sample Complexity of Dictionary Learning and Other Matrix Factorizations.  _IEEE Transactions on Information Theory_ , 6(61):3469–3486, 2015.  ISSN 0018-9448, 1557-9654.  doi: 10.1109/TIT.2015.2424238.  URL https://www.infona.pl//resource/bwmeta1.element.ieee-art-000007088631.
  * Hajnal et al. (2023) Hajnal, M. A., Tran, D., Einstein, M., Martelo, M. V., Safaryan, K., Polack, P.-O., Golshani, P., and Orbán, G.  Continuous multiplexed population representations of task context in the mouse primary visual cortex.  _Nature Communications_ , 14(1):6687, October 2023.  ISSN 2041-1723.  doi: 10.1038/s41467-023-42441-w.  URL https://www.nature.com/articles/s41467-023-42441-w.  Number: 1 Publisher: Nature Publishing Group.
  * Hallin et al. (2023) Hallin, M., Nisol, G., and Tavakoli, S.  Factor models for high-dimensional functional time series I: Representation results.  _Journal of Time Series Analysis_ , 44(5-6):578–600, 2023.  ISSN 1467-9892.  doi: 10.1111/jtsa.12676.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/jtsa.12676.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/jtsa.12676.
  * Hespanha (2018) Hespanha, J. P.  _Linear Systems Theory: Second Edition_.  Princeton University Press, February 2018.  ISBN 978-1-4008-9008-8.
  * Jedra & Proutiere (2023) Jedra, Y. and Proutiere, A.  Finite-Time Identification of Linear Systems: Fundamental Limits and Optimal Algorithms.  _IEEE Transactions on Automatic Control_ , 68(5):2805–2820, May 2023.  ISSN 1558-2523.  doi: 10.1109/TAC.2022.3221705.  URL https://ieeexplore.ieee.org/document/9946382.  Conference Name: IEEE Transactions on Automatic Control.
  * Lee (2022) Lee, H.  Improved rates for prediction and identification of partially observed linear dynamical systems, March 2022.  URL http://arxiv.org/abs/2011.10006.  arXiv:2011.10006 [cs, eess, math, stat].
  * Li (1998) Li, R.-C.  Relative Perturbation Theory: II. Eigenspace and Singular Subspace Variations.  _SIAM Journal on Matrix Analysis and Applications_ , 20(2):471–492, June 1998.  ISSN 0895-4798.  doi: 10.1137/S0895479896298506.  URL https://epubs.siam.org/doi/abs/10.1137/S0895479896298506.  Publisher: Society for Industrial and Applied Mathematics.
  * Maliar & Maliar (2015) Maliar, L. and Maliar, S.  Merging simulation and projection approaches to solve high-dimensional problems with an application to a new Keynesian model.  _Quantitative Economics_ , 6(1):1–47, 2015.  ISSN 1759-7331.  doi: 10.3982/QE364.  URL https://onlinelibrary.wiley.com/doi/abs/10.3982/QE364.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.3982/QE364.
  * Mao et al. (2021) Mao, Y., Hovakimyan, N., Voulgaris, P., and Sha, L.  Finite-Time Model Inference From A Single Noisy Trajectory, January 2021.  URL http://arxiv.org/abs/2010.06616.  arXiv:2010.06616 [cs, eess].
  * Marks & Goard (2021) Marks, T. D. and Goard, M. J.  Stimulus-dependent representational drift in primary visual cortex.  _Nature Communications_ , 12(1):5169, August 2021.  ISSN 2041-1723.  doi: 10.1038/s41467-021-25436-3.  URL https://www.nature.com/articles/s41467-021-25436-3.  Number: 1 Publisher: Nature Publishing Group.
  * Masini et al. (2023) Masini, R. P., Medeiros, M. C., and Mendes, E. F.  Machine learning advances for time series forecasting.  _Journal of Economic Surveys_ , 37(1):76–111, 2023.  ISSN 1467-6419.  doi: 10.1111/joes.12429.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/joes.12429.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/joes.12429.
  * Medsker & Jain (1999) Medsker, L. and Jain, L. C.  _Recurrent Neural Networks: Design and Applications_.  CRC Press, December 1999.  ISBN 978-1-4200-4917-6.  Google-Books-ID: ME1SAkN0PyMC.
  * Mudassir et al. (2020) Mudassir, M., Bennbaia, S., Unal, D., and Hammoudeh, M.  Time-series forecasting of Bitcoin prices using high-dimensional features: a machine learning approach.  _Neural Computing and Applications_ , July 2020.  ISSN 1433-3058.  doi: 10.1007/s00521-020-05129-6.  URL https://doi.org/10.1007/s00521-020-05129-6.
  * Oymak & Ozay (2019) Oymak, S. and Ozay, N.  Non-asymptotic Identification of LTI Systems from a Single Trajectory, February 2019.  URL http://arxiv.org/abs/1806.05722.  arXiv:1806.05722 [cs, math, stat].
  * Pandarinath et al. (2018) Pandarinath, C., O’Shea, D. J., Collins, J., Jozefowicz, R., Stavisky, S. D., Kao, J. C., Trautmann, E. M., Kaufman, M. T., Ryu, S. I., Hochberg, L. R., Henderson, J. M., Shenoy, K. V., Abbott, L. F., and Sussillo, D.  Inferring single-trial neural population dynamics using sequential auto-encoders.  _Nature Methods_ , 15(10):805–815, October 2018.  ISSN 1548-7105.  doi: 10.1038/s41592-018-0109-9.  URL https://www.nature.com/articles/s41592-018-0109-9.  Number: 10 Publisher: Nature Publishing Group.
  * Poloni & Sbrana (2019) Poloni, F. and Sbrana, G.  Closed-form results for vector moving average models with a univariate estimation approach.  _Econometrics and Statistics_ , 10:27–52, April 2019.  ISSN 2452-3062.  doi: 10.1016/j.ecosta.2018.06.003.  URL https://www.sciencedirect.com/science/article/pii/S2452306218300327.
  * Qin (2022) Qin, S. J.  Latent vector autoregressive modeling and feature analysis of high dimensional and noisy data from dynamic systems.  _AIChE Journal_ , 68:e17703, June 2022.  ISSN 0001-1541.  doi: 10.1002/aic.17703.  URL https://ui.adsabs.harvard.edu/abs/2022AIChE..68E7703Q.  ADS Bibcode: 2022AIChE..68E7703Q.
  * Rudelson (1999) Rudelson, M.  Random Vectors in the Isotropic Position.  _Journal of Functional Analysis_ , 164(1):60–72, May 1999.  ISSN 0022-1236.  doi: 10.1006/jfan.1998.3384.  URL https://www.sciencedirect.com/science/article/pii/S0022123698933845.
  * Sarkar & Rakhlin (2019) Sarkar, T. and Rakhlin, A.  Near optimal finite time identification of arbitrary linear dynamical systems.  In _Proceedings of the 36th International Conference on Machine Learning_ , pp. 5610–5618. PMLR, May 2019.  URL https://proceedings.mlr.press/v97/sarkar19a.html.  ISSN: 2640-3498.
  * Sarkar et al. (2021) Sarkar, T., Rakhlin, A., and Dahleh, M. A.  Finite time LTI system identification.  _The Journal of Machine Learning Research_ , 22(1):26:1186–26:1246, January 2021.  ISSN 1532-4435.
  * (39) Saul, L. K., Labs, T., Ave, P., Park, F., and Roweis, S. T.  An Introduction to Locally Linear Embedding.
  * Schutter (2000) Schutter, B. D.  Minimal state-space realization in linear system theory: an overview.  _Journal of Computational and Applied Mathematics_ , 121(1):331–354, September 2000.  ISSN 0377-0427.  doi: 10.1016/S0377-0427(00)00341-1.  URL https://www.sciencedirect.com/science/article/pii/S0377042700003411.
  * Shirani Faradonbeh et al. (2018) Shirani Faradonbeh, M. K., Tewari, A., and Michailidis, G.  Finite time identification in unstable linear systems.  _Automatica_ , 96:342–353, October 2018.  ISSN 0005-1098.  doi: 10.1016/j.automatica.2018.07.008.  URL https://www.sciencedirect.com/science/article/pii/S0005109818303546.
  * Sikander & Prasad (2015) Sikander, A. and Prasad, R.  Linear time-invariant system reduction using a mixed methods approach.  _Applied Mathematical Modelling_ , 39(16):4848–4858, August 2015.  ISSN 0307-904X.  doi: 10.1016/j.apm.2015.04.014.  URL https://www.sciencedirect.com/science/article/pii/S0307904X15002504.
  * Simchowitz et al. (2018) Simchowitz, M., Mania, H., Tu, S., Jordan, M. I., and Recht, B.  Learning Without Mixing: Towards A Sharp Analysis of Linear System Identification.  In _Proceedings of the 31st Conference On Learning Theory_ , pp. 439–473. PMLR, July 2018.  URL https://proceedings.mlr.press/v75/simchowitz18a.html.  ISSN: 2640-3498.
  * Stewart (1990) Stewart, G. W. G. W.  _Matrix perturbation theory_.  Computer science and scientific computing. Academic Press, Boston, 1990.  ISBN 0-12-670230-6.
  * Stringer et al. (2019) Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D.  High-dimensional geometry of population responses in visual cortex.  _Nature_ , 571(7765):361–365, July 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-1346-5.
  * Sun et al. (2023) Sun, S., Li, J., and Mo, Y.  Finite Time Performance Analysis of MIMO Systems Identification, October 2023.  URL http://arxiv.org/abs/2310.11790.  arXiv:2310.11790 [cs, eess].
  * Sun et al. (2022) Sun, Y., Oymak, S., and Fazel, M.  Finite Sample Identification of Low-Order LTI Systems via Nuclear Norm Regularization.  _IEEE Open Journal of Control Systems_ , 1:237–254, 2022.  ISSN 2694-085X.  doi: 10.1109/OJCSYS.2022.3200015.  URL https://ieeexplore.ieee.org/document/9870857.  Conference Name: IEEE Open Journal of Control Systems.
  * Sussillo & Barak (2013) Sussillo, D. and Barak, O.  Opening the Black Box: Low-Dimensional Dynamics in High-Dimensional Recurrent Neural Networks.  _Neural Computation_ , 25(3):626–649, March 2013.  ISSN 0899-7667.  doi: 10.1162/NECO˙a˙00409.  URL https://doi.org/10.1162/NECO_a_00409.
  * Talebi et al. (2023) Talebi, S., Taghvaei, A., and Mesbahi, M.  Data-driven Optimal Filtering for Linear Systems with Unknown Noise Covariances, October 2023.  URL http://arxiv.org/abs/2305.17836.  arXiv:2305.17836 [cs, eess, math].
  * Tian et al. (2023) Tian, Y., Zhang, K., Tedrake, R., and Sra, S.  Toward Understanding Latent Model Learning in MuZero: A Case Study in Linear Quadratic Gaussian Control.  July 2023.  URL https://openreview.net/forum?id=r9YZ357Trz.
  * Tripuraneni et al. (2021) Tripuraneni, N., Jin, C., and Jordan, M.  Provable Meta-Learning of Linear Representations.  In _Proceedings of the 38th International Conference on Machine Learning_ , pp. 10434–10443. PMLR, July 2021.  URL https://proceedings.mlr.press/v139/tripuraneni21a.html.  ISSN: 2640-3498.
  * Vaswani et al. (2018) Vaswani, N., Bouwmans, T., Javed, S., and Narayanamurthy, P.  Robust Subspace Learning: Robust PCA, Robust Subspace Tracking, and Robust Subspace Recovery.  _IEEE Signal Processing Magazine_ , 35:32–55, July 2018.  ISSN 1053-5888.  doi: 10.1109/MSP.2018.2826566.  URL https://ui.adsabs.harvard.edu/abs/2018ISPM...35d..32V.  ADS Bibcode: 2018ISPM…35d..32V.
  * Vershynin (2018) Vershynin, R.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2018.  ISBN 978-1-108-41519-4.  doi: 10.1017/9781108231596.  URL https://www.cambridge.org/core/books/highdimensional-probability/797C466DA29743D2C8213493BD2D2102.
  * Xia et al. (2021) Xia, J., Marks, T. D., Goard, M. J., and Wessel, R.  Stable representation of a naturalistic movie emerges from episodic activity with gain variability.  _Nature Communications_ , 12(1):5170, August 2021.  ISSN 2041-1723.  doi: 10.1038/s41467-021-25437-2.  URL https://www.nature.com/articles/s41467-021-25437-2.  Number: 1 Publisher: Nature Publishing Group.
  * Yu et al. (2021) Yu, W., Kim, I. Y., and Mechefske, C.  Analysis of different RNN autoencoder variants for time series classification and machine prognostics.  _Mechanical Systems and Signal Processing_ , 149:107322, February 2021.  ISSN 0888-3270.  doi: 10.1016/j.ymssp.2020.107322.  URL https://www.sciencedirect.com/science/article/pii/S0888327020307081.
  * Zhang et al. (2023) Zhang, T. T. C. K., Toso, L. F., Anderson, J., and Matni, N.  Meta-Learning Operators to Optimality from Multi-Task Non-IID Data, August 2023.  URL http://arxiv.org/abs/2308.04428.  arXiv:2308.04428 [cs, eess, stat].
  * Zheng & Li (2021) Zheng, Y. and Li, N.  Non-Asymptotic Identification of Linear Dynamical Systems Using Multiple Trajectories.  _IEEE Control Systems Letters_ , 5(5):1693–1698, November 2021.  ISSN 2475-1456.  doi: 10.1109/LCSYS.2020.3042924.  URL https://ieeexplore.ieee.org/document/9284539.  Conference Name: IEEE Control Systems Letters.

##  Appendix A Upper Bounds for HD-SYSID — Proof of Theorem 3.4

Recall the setting in HD-SYSID. We consider system
$\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\sigma_{\eta}^{2}I)$ and inputs
$\mathcal{U}_{1}=\\{u_{1,t}\\}_{t=0}^{T_{1}-1}$,
$\mathcal{U}_{2}=\\{u_{2,t}\\}_{t=0}^{T_{2}-1}$ sampled independently from
$\mathcal{N}(0,\Sigma_{u})$. To simplify future analysis, we define the
following notations

$\begin{split}{}&\psi_{C}=\sigma_{1}\left(C\right),\quad\psi_{\eta}=\sigma_{1}(\sigma_{\eta}^{2}I),\quad\psi_{w}=\sigma_{1}(\Sigma_{w}+B\Sigma_{u}B{}^{\intercal}),\\\
{}&\phi_{C}=\sigma_{\min}(C),\quad\phi_{O}=\sigma_{\min}\left(\begin{bmatrix}C\\\
CA\\\ \vdots\\\
CA^{r-1}\end{bmatrix}\right),\quad\phi_{R}=\sigma_{\min}\left(\begin{bmatrix}B&AB&\dots&A^{r-1}B\end{bmatrix}\right),\quad\phi_{u}=\sigma_{\min}(\Sigma_{u}).\\\
\end{split}$ (4)

Here we assume all $\psi$’s satisfy $\psi\geq 1$, otherwise we define $\psi$
to be $\max\\{1,\sigma_{1}(\cdot)\\}$. Similarly, we assume all $\phi$’s
satisfy $\phi\leq 1$, otherwise we define $\phi$ to be
$\min\\{1,\sigma_{\min}(\cdot)\\}$.

Recall that auxiliary system $\widehat{\mathcal{M}}$ is defined as follows
with $\widehat{\Phi}_{C}$ being the approximated observer column space:

$\begin{split}x_{t+1}={}&Ax_{t}+Bu_{t}+w_{t},\\\
y_{t}={}&\widehat{\Phi}_{C}{}^{\intercal}Cx_{t}+\widehat{\Phi}_{C}{}^{\intercal}\eta_{t}.\end{split}$
(5)

Now we are ready to restate Theorem 3.4 in full details.

######  Theorem A.1 (Theorem 3.4 Restated).

Consider $\mathcal{M}$, datasets
$\mathcal{D}_{1}=\mathcal{U}_{1}\cup\mathcal{Y}_{1},\mathcal{D}_{2}=\mathcal{U}_{2}\cup\mathcal{Y}_{2}$
in HD-SYSID and constants defined above. Suppose system $\mathcal{M}$
satisfies Assumption 3.1 with constants $\psi_{A}$ and $\rho_{A}$. If $T_{1}$
and $T_{2}$ satisfy

$\begin{split}&T_{1}\gtrsim\kappa_{3}\cdot n^{2}r^{3},\quad
T_{2}\geq\kappa_{1}\cdot\text{poly}\left(r,m\right),\end{split}$ (6)

then $(\widehat{A},\widehat{B},\widehat{C})$ from Algorithm 1 satisfy the
following for some invertible matrix $S$ with probability at least $1-\delta$

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\lesssim\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|+\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,m\right)}{T_{2}}}.\end{split}$
(7)

Here $\kappa_{1}=\kappa_{1}(\widehat{\mathcal{M}},\mathcal{U}_{2},\delta)$ and
$\kappa_{2}=\kappa_{2}(\widehat{\mathcal{M}},\mathcal{U}_{2},\delta)$ are
defined in Definition 3.2.
$\kappa_{3}=\kappa_{3}\left(\mathcal{M},\mathcal{U}_{[2]},\delta\right),\quad\kappa_{4}=\kappa_{4}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)$
are detailed below. All of them are problem-related constants independent of
system dimensions modulo logarithmic factors.

$\begin{split}\kappa_{3}\left(\mathcal{M},\mathcal{U}_{[2]},\delta\right)={}&\max\left\\{\kappa_{4}^{2}\frac{\psi_{A}^{2}\psi_{C}^{2}}{(1-\rho_{A}^{2})\phi_{O}^{2}},~{}\left(\frac{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{(1-\rho_{A}^{2})\phi_{u}^{2}\phi_{C}^{4}\phi_{R}^{4}}\right)^{2}\log^{2}\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)\log^{4}(\frac{r}{\delta})\right\\},\\\
\kappa_{4}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)={}&\frac{\psi_{\eta}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\log\frac{1}{\delta}}.\end{split}$
(8)

###### Proof.

Based on LABEL:eq:3col_1, $T_{1}$ satisfies the condition of Lemma A.2. We
apply Lemma A.2 on $(\mathcal{D}_{1},\Sigma_{\eta})$ and get the following
with probability at least $1-\frac{\delta}{2}$

$\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\lesssim\kappa_{4}\sqrt{\frac{n}{T_{1}}}\coloneqq\Delta_{\Phi}.$
(9)

The system generating dataset $\tilde{\mathcal{D}}_{2}$ (line 3 in Algorithm
1), denoted by $\widehat{\mathcal{M}}$, is rewritten as in Equation 5.

##### Step 1. We first show that $\widehat{\mathcal{M}}$ is still a minimal
system.

The controllability directly comes from the fact that
$R=\begin{bmatrix}B&AB&\dots&A^{r-1}B\end{bmatrix}$ is full row rank, because
system $\mathcal{M}$ is minimal. On the other hand, for the observability
matrix $O$, we know that

$\begin{split}{}&\text{rank}(O)=\text{rank}\left(\begin{bmatrix}\widehat{\Phi}_{C}{}^{\intercal}C\\\
\widehat{\Phi}_{C}{}^{\intercal}CA\\\ \vdots\\\
\widehat{\Phi}_{C}{}^{\intercal}CA^{r-1}\end{bmatrix}\right)\\\
\geq{}&\text{rank}\left(\text{diag}\left(\widehat{\Phi}_{C},\dots,\widehat{\Phi}_{C}\right)\begin{bmatrix}\widehat{\Phi}_{C}{}^{\intercal}C\\\
\widehat{\Phi}_{C}{}^{\intercal}CA\\\ \vdots\\\
\widehat{\Phi}_{C}{}^{\intercal}CA^{r-1}\end{bmatrix}\right)=\text{rank}\left(\begin{bmatrix}\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}C\\\
\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CA\\\ \vdots\\\
\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CA^{r-1}\end{bmatrix}\right)\\\
={}&\text{rank}\left(\begin{bmatrix}C\\\ CA\\\ \vdots\\\
CA^{r-1}\end{bmatrix}-\begin{bmatrix}\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}C\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA\\\
\vdots\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{r-1}\end{bmatrix}\right).\end{split}$
(10)

Consider the second term. We first observe that the following holds for all
$i\in[r-1]$

$\begin{split}\left\|\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{i}\right\|={}&\left\|\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\Phi_{C}{}^{\intercal}CA^{i}\right\|\leq\left\|\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\left\|\Phi_{C}{}^{\intercal}CA^{i}\right\|\\\
={}&\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\left\|CA^{i}\right\|\\\
\leq{}&\Delta_{\Phi}\psi_{A}\psi_{C}\rho_{A}^{i-1}.\end{split}$ (11)

Moreover, for $i=0$, we know

$\begin{split}\left\|\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}C\right\|\leq\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\left\|C\right\|\leq\Delta_{\Phi}\psi_{C}.\end{split}$
(12)

Therefore,

$\begin{split}\left\|\begin{bmatrix}\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}C\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA\\\
\vdots\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{r-1}\end{bmatrix}\right\|={}&\sqrt{\left\|\sum_{i=0}^{r-1}A^{i}{}^{\intercal}C{}^{\intercal}\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{i}\right\|}\\\
\leq{}&\sqrt{\sum_{i=0}^{r-1}\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{i}\right\|^{2}}\\\
\leq{}&\sqrt{(\Delta_{\Phi}\psi_{C})^{2}+\sum_{i=1}^{r-1}\left(\Delta_{\Phi}\psi_{A}\psi_{C}\rho_{A}^{i-1}\right)^{2}}\\\
\leq{}&2\frac{\psi_{A}\psi_{C}}{\sqrt{1-\rho_{A}^{2}}}\Delta_{\Phi}\end{split}$
(13)

From LABEL:eq:3col_1, we know that
$T_{1}\gtrsim\kappa_{4}^{2}\frac{\psi_{A}^{2}\psi_{C}^{2}}{(1-\rho_{A}^{2})\phi_{O}^{2}}n$.
Combining Theorem 1 in (Stewart, 1990) gives

$\begin{split}{}&\sigma_{r}\left(\begin{bmatrix}C\\\ CA\\\ \vdots\\\
CA^{r-1}\end{bmatrix}-\begin{bmatrix}\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}C\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA\\\
\vdots\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{r-1}\end{bmatrix}\right)\\\
\geq{}&\sigma_{r}\left(\begin{bmatrix}C\\\ CA\\\ \vdots\\\
CA^{r-1}\end{bmatrix}\right)-\left\|\begin{bmatrix}\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}C\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA\\\
\vdots\\\
\widehat{\Phi}_{C}^{\perp}\widehat{\Phi}_{C}^{\perp}{}^{\intercal}CA^{r-1}\end{bmatrix}\right\|\\\
\geq{}&\frac{\phi_{O}}{2}>0.\end{split}$ (14)

This implies that

$\begin{split}r\geq\text{rank}(O)\geq r.\end{split}$ (15)

Namely, the system is observable. Since
$\left(A,B,\widehat{\Phi}_{C}{}^{\intercal}C\right)$ is controllable and
observable, we conclude that $\widehat{\mathcal{M}}$ is minimal.

Step 2. We now apply Sys-Oracle on this minimal system.  Recall the dynamics
of $\widehat{\mathcal{M}}$:

$\begin{split}x_{t+1}={}&Ax_{t}+Bu_{t}+w_{t},\\\
y_{t}={}&\widehat{\Phi}_{C}{}^{\intercal}Cx_{t}+\widehat{\Phi}_{C}{}^{\intercal}\eta_{t}.\end{split}$
(16)

Since $\widehat{\Phi}_{C}$ is independent of the second trajectory,
$\\{\widehat{\Phi}_{C}\eta_{2,t}\\}_{t=0}^{T_{2}}$ are i.i.d. noises
independent of other variables of the second trajectory. Moreover,
$\widehat{\mathcal{M}}$ satisfy Assumption 3.1. Therefore, we can apply Sys-
Oracle. Let $r_{c}=\text{rank}(C)$. With
$T_{2}\geq\kappa_{1}\text{poly}\left(r,r_{c},m\right)=\kappa_{1}\text{poly}\left(r,m\right)$,
outputs $\widehat{A},\widehat{B},\tilde{C}$ satisfy the following for some
invertible matrix $S$ with probability at least $1-\frac{\delta}{2}$

$\begin{split}{}&\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|\widehat{\Phi}_{C}{}^{\intercal}CS-\tilde{C}\right\|\right\\}\leq\kappa_{2}\sqrt{\frac{\text{poly}\left(r,r_{c},m\right)}{T_{2}}}\eqqcolon\Delta_{\widehat{\mathcal{M}}}.\end{split}$
(17)

Therefore, our final approximation for $C$ satisfies

$\begin{split}\left\|CS-\widehat{C}\right\|={}&\left\|CS-\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS+\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS-\widehat{\Phi}_{C}\tilde{C}\right\|\\\
\leq{}&\left\|CS-\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS\right\|+\left\|\widehat{\Phi}_{C}\widehat{\Phi}_{C}{}^{\intercal}CS-\widehat{\Phi}_{C}\tilde{C}\right\|\\\
\leq{}&\left\|\widehat{\Phi}_{C}^{\perp}\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}CS\right\|+\left\|\widehat{\Phi}_{C}{}^{\intercal}CS-\tilde{C}\right\|\\\
\leq{}&\left\|\widehat{\Phi}_{C}^{\perp}\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}CS\right\|+\Delta_{\widehat{\mathcal{M}}}.\end{split}$
(18)

For the first term, we know that

$\begin{split}\left\|\widehat{\Phi}_{C}^{\perp}\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}CS\right\|={}&\left\|\widehat{\Phi}_{C}^{\perp}\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}\Phi_{C}\Phi_{C}{}^{\intercal}CS\right\|=\left\|\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}\Phi_{C}\Phi_{C}{}^{\intercal}CS\right\|\\\
\leq{}&\left\|\left(\widehat{\Phi}_{C}^{\perp}\right){}^{\intercal}\Phi_{C}\right\|\left\|\Phi_{C}{}^{\intercal}\right\|\left\|CS\right\|\\\
\leq{}&\Delta_{\Phi}\left\|CS\right\|.\end{split}$ (19)

Therefore,

$\begin{split}\left\|CS-\widehat{C}\right\|\leq\Delta_{\Phi}\left\|CS\right\|+\Delta_{\widehat{\mathcal{M}}}.\end{split}$
(20)

To get the upper bound w.r.t. $\left\|\widehat{C}\right\|$, we notice that

$\begin{split}\left\|CS\right\|\leq\left\|\widehat{C}\right\|+\left\|CS-\widehat{C}\right\|\leq\left\|\widehat{C}\right\|+\Delta_{\widehat{\mathcal{M}}}+\Delta_{\Phi}\left\|CS\right\|.\end{split}$
(21)

Rearraging the terms gives

$\begin{split}\left\|CS\right\|\leq\frac{\left\|\widehat{C}\right\|+\Delta_{\widehat{\mathcal{M}}}}{1-\Delta_{\Phi}}.\end{split}$
(22)

Substituting back gives the following with probability at least $1-\delta$

$\begin{split}\left\|CS-\widehat{C}\right\|\leq{}&\frac{1}{1-\Delta_{\Phi}}\Delta_{\widehat{\mathcal{M}}}+\frac{\Delta_{\Phi}}{1-\Delta_{\Phi}}\left\|\widehat{C}\right\|\leq
2\Delta_{\widehat{\mathcal{M}}}+\frac{\Delta_{\Phi}}{1-\Delta_{\Phi}}\left\|\widehat{C}\right\|\\\
\overset{(i)}{\leq}{}&2\Delta_{\widehat{\mathcal{M}}}+2\Delta_{\Phi}\left\|\widehat{C}\right\|\\\
\lesssim{}&\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,r_{c},m\right)}{T_{2}}}+\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|\\\
\leq{}&\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,m\right)}{T_{2}}}+\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|\\\
\end{split}$ (23)

Here $(i)$ is because $\Delta_{\Phi}\leq 1/2$ due to LABEL:eq:3col_1. Finally,
we conclude that

$\max\left\\{\left\|S^{-1}AS-\widehat{A}\right\|,\left\|S^{-1}B-\widehat{B}\right\|,\left\|CS-\widehat{C}\right\|\right\\}\lesssim\kappa_{2}\cdot\sqrt{\frac{\text{poly}\left(r,m\right)}{T_{2}}}+\kappa_{4}\cdot\sqrt{\frac{n}{T_{1}}}\left\|\widehat{C}\right\|.$
(24)

∎

###  A.1 Upper Bounds for Col-Approx

The theoretical guarantee for Col-Approx is presented in the following lemma.

######  Lemma A.2 (Lemma 3.7 Restated).

Consider system $\mathcal{M}$, dataset
$\mathcal{D}_{1}=\mathcal{U}_{1}\cup\mathcal{Y}_{1}$ in HD-SYSID and constants
defined at the beginning of Appendix A. Suppose $\mathcal{M}$ satisfies
Assumption 3.1 with constants $\psi_{A}$ and $\rho_{A}$. If

$\begin{split}T_{1}\gtrsim\underbrace{\left(\frac{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{(1-\rho_{A}^{2})\phi_{u}^{2}\phi_{C}^{4}\phi_{R}^{4}}\right)^{2}\log^{2}\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)\log^{4}(\frac{r}{\delta})}_{\kappa_{5}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)}\cdot
n^{2}r^{3},\end{split}$ (25)

then $\widehat{\Phi}_{C}=\texttt{Col-Approx}(\mathcal{Y}_{1},\Sigma_{\eta})$
satisfies the following with probability at least $1-\delta$

$\begin{split}\widehat{r}_{c}=\text{rank}(C),\quad\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|\lesssim\underbrace{\frac{\psi_{\eta}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\log\frac{1}{\delta}}}_{\kappa_{4}\left(\mathcal{M},\mathcal{U}_{1},\delta\right)}\cdot\sqrt{\frac{n}{T_{1}}}.\end{split}$
(26)

###### Proof.

For simplicity, we omit all subscript $1$ for the rest of this section. From
the system dynamics, we know that

$\begin{split}\Sigma_{y}={}&\sum_{t=0}^{T}y_{t}y_{t}{}^{\intercal}=\sum_{t=0}^{T}Cx_{t}x_{t}{}^{\intercal}C{}^{\intercal}+\sum_{t=0}^{T}\eta_{t}\eta_{t}{}^{\intercal}+\sum_{t=0}^{T}\left(Cx_{t}\eta_{t}{}^{\intercal}+\eta_{t}x_{t}{}^{\intercal}C{}^{\intercal}\right)\\\
={}&\sum_{t=0}^{T}Cx_{t}x_{t}{}^{\intercal}C{}^{\intercal}+\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)+\sum_{t=0}^{T}\left(Cx_{t}\eta_{t}{}^{\intercal}+\eta_{t}x_{t}{}^{\intercal}C{}^{\intercal}\right)+(T+1)\sigma_{\eta}^{2}I\\\
\end{split}$ (27)

Here the first term is the information on $\text{col}(C)$, while the second
and third terms are noises. For the rest of the proof, we first upper bound
norms of the noise terms (step 1). With this, we show that
$\widehat{r}_{c}=\text{rank}(C)$ with high probability (step 2). We then apply
our subspace perturbation result to upper bound the influence of the noises on
the eigenspace of the first term (step 3).

Step 1: Noise Norm Upper Bounds. Define $r_{c}=\text{rank}(C)$. Notice that we
can write $C=\Phi_{C}\alpha$, where $\Phi_{C}\in\mathbb{R}^{n\times r_{c}}$
consists of orthonormal columns that form a basis of $\text{col}(C)$ and
$\alpha\in\mathbb{R}^{r_{c}\times r}$ is a full row rank matrix. It is then
clear that

$\begin{split}\sigma_{\min}(\alpha)=\sigma_{\min}(C)\geq\phi_{C},\quad\sigma_{1}(\alpha)=\sigma_{1}(C)\leq\psi_{C}.\end{split}$
(28)

Let $\Sigma_{C}=\sum_{t=0}^{T}Cx_{t}x_{t}{}^{\intercal}C{}^{\intercal}$,
$\bar{\Sigma}_{C}=\Sigma_{C}+TI$ and $\Sigma_{\alpha}=\sum_{t=0}^{T}\alpha
x_{t}x_{t}{}^{\intercal}\alpha{}^{\intercal}$. From the definitions, it is
clear that $\Sigma_{C}=\Phi_{C}\Sigma_{\alpha}\Phi_{C}{}^{\intercal}$ and
$\bar{\Sigma}_{C}=\Phi_{C}\Sigma_{\alpha}\Phi_{C}{}^{\intercal}+TI$. Then from
Lemma A.3, Lemma A.4, Lemma A.5, and Lemma A.6, we have the following for
$T\gtrsim\frac{\psi_{w}^{2}\psi_{C}^{4}\psi_{A}^{4}}{\phi_{u}^{2}\phi_{C}^{4}\phi_{R}^{4}}r^{3}\log\frac{r}{\delta}\cdot\tilde{U}^{2}$
(from Equation 25) with probability at least $1-\delta$

$\begin{split}{}&\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{8}TI\preceq\Sigma_{\alpha}\precsim\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}rT\log\frac{1}{\delta}I,\\\
{}&\left\|\left(\bar{\Sigma}_{C}\right)^{-\frac{1}{2}}\sum_{t=0}^{T}Cx_{t}\eta_{t}{}^{\intercal}\right\|\lesssim\tilde{U}\sqrt{\psi_{\eta}n\log\frac{1}{\delta}},\\\
{}&\left\|\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)\right\|\lesssim\psi_{\eta}\sqrt{nT\log\frac{1}{\delta}}.\end{split}$
(29)

Here
$\tilde{U}=\sqrt{\log\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)}$.

Step 2: Order Estimation Guarantee. Consider the following matrix

$\begin{split}\Sigma_{y}={}&\Sigma_{C}+\underbrace{\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)+\sum_{t=0}^{T}\left(Cx_{t}\eta_{t}{}^{\intercal}+\eta_{t}C{}^{\intercal}x_{t}{}^{\intercal}\right)}_{\Delta}+(T+1)\sigma_{\eta}^{2}I.\end{split}$
(30)

The inequalities of Step 1 imply

$\begin{split}\left\|\Delta\right\|\leq{}&\left\|\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)\right\|+2\left\|\sum_{t=0}^{T}Cx_{t}\eta_{t}{}^{\intercal}\right\|\\\
\lesssim{}&\psi_{\eta}\sqrt{nT\log\frac{1}{\delta}}+\tilde{U}\sqrt{\psi_{\eta}n\log\frac{1}{\delta}}\left\|(\bar{\Sigma}_{C})^{\frac{1}{2}}\right\|\\\
={}&\psi_{\eta}\sqrt{nT\log\frac{1}{\delta}}+\tilde{U}\sqrt{\psi_{\eta}n\log\frac{1}{\delta}}\left(\left\|(TI+\Sigma_{\alpha})^{\frac{1}{2}}\right\|\right)\\\
\lesssim{}&\psi_{\eta}\sqrt{nT\log\frac{1}{\delta}}+\tilde{U}\sqrt{\psi_{\eta}n\log\frac{1}{\delta}}\sqrt{\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}rT\log\frac{1}{\delta}+T}\\\
\lesssim{}&\tilde{U}\psi_{\eta}\sqrt{n\log\frac{1}{\delta}}\sqrt{\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}rT\log\frac{1}{\delta}}\\\
={}&\frac{\sqrt{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\sqrt{1-\rho_{A}^{2}}}\tilde{U}\log\left(\frac{1}{\delta}\right)\sqrt{nrT}\\\
\end{split}$ (31)

Therefore, for each $i\in[r_{c}]$, we have the following for some positive
constant $c_{1}$:

$\begin{split}\sigma_{i}\left(\Sigma_{y}\right)-(T+1)\sigma_{\eta}^{2}\geq{}&\sigma_{i}(\Sigma_{C})-\left\|\Delta\right\|\geq\sigma_{\min}(\Sigma_{\alpha})-\left\|\Delta\right\|\\\
\geq{}&\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{8}T-c_{1}\frac{\sqrt{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\sqrt{1-\rho_{A}^{2}}}\tilde{U}\log\left(\frac{1}{\delta}\right)\sqrt{nrT}\\\
\geq{}&\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{16}T.\end{split}$ (32)

Here the first inequality holds according to Theorem 1 in (Stewart, 1990) and
the last line is because
$T\gtrsim\frac{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{(1-\rho_{A}^{2})\phi_{u}^{2}\phi_{C}^{4}\phi_{R}^{4}}\log\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)\log^{2}(\frac{1}{\delta})\cdot
nr$ (from Equation 25). For $i\in[r_{c}+1,n]$,

$\begin{split}\sigma_{i}(\Sigma_{y})-(T+1)\sigma_{\eta}^{2}\leq\left\|\Delta\right\|\leq
c_{1}\frac{\sqrt{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\sqrt{1-\rho_{A}^{2}}}\tilde{U}\log\left(\frac{1}{\delta}\right)\sqrt{nrT}.\end{split}$
(33)

Based on the above three inequalities and Equation 25, we conclude that with
probability at least $1-\delta$, the following hold

$\begin{split}\sigma_{i}(\Sigma_{y})-\sigma_{j}(\Sigma_{y})\leq{}&2c_{1}\frac{\sqrt{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\sqrt{1-\rho_{A}^{2}}}\tilde{U}\log\left(\frac{1}{\delta}\right)\sqrt{nrT}\\\
<{}&T^{3/4},\quad\forall i<j\in[r_{c}+1,n]\\\
\sigma_{r_{c}}(\Sigma_{y})-\sigma_{r_{c}+1}(\Sigma_{y})\geq{}&\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{16}T-c_{1}\frac{\sqrt{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\sqrt{1-\rho_{A}^{2}}}\tilde{U}\log\left(\frac{1}{\delta}\right)\sqrt{nrT}\geq\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{32}T\\\
>{}&T^{3/4}.\end{split}$

The above is because
$T\gtrsim\left(\frac{\psi_{\eta}^{2}\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{(1-\rho_{A}^{2})\phi_{u}^{2}\phi_{C}^{4}\phi_{R}^{4}}\right)^{2}\log^{2}\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)\log^{4}(\frac{1}{\delta})\cdot
n^{2}r^{2}$. Therefore, from the definition of $\widehat{r}_{c}$, we know that
$\widehat{r}_{c}=r_{c}$ with probability at least $1-\delta$.

Step 3: Column Space Estimation Guarantee. With $\widehat{r}_{c}=r_{c}$, now
we try to apply our subspace perturbation result, i.e. Lemma A.7, on matrix
$\Sigma_{y}-\sigma_{\eta}^{2}(T+1)I+TI$. Notice that this matrix has exactly
the same eigenspace as $\Sigma_{y}$ and therefore the eigenspace of its first
$r_{c}$ eigenvectors is $\widehat{\Phi}_{C}$ (line 9 in Algorithm 2). This
matrix can be decomposed as

$\begin{split}\Sigma_{y}-\sigma_{\eta}^{2}(T+1)I+TI={}&\Sigma_{C}+\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)+\sum_{t=0}^{T}\left(Cx_{t}\eta_{t}{}^{\intercal}+\eta_{t}C{}^{\intercal}x_{t}{}^{\intercal}\right)+TI.\\\
={}&\underbrace{\bar{\Sigma}_{C}}_{M\text{ in Lemma
\ref{lem:perturbation}}}+\underbrace{\sum_{t=0}^{T}\left(\eta_{t}\eta_{t}{}^{\intercal}-\sigma_{\eta}^{2}I\right)}_{\Delta_{2}\text{
in Lemma
\ref{lem:perturbation}}}+\underbrace{\sum_{t=0}^{T}\left(Cx_{t}\eta_{t}{}^{\intercal}+\eta_{t}C{}^{\intercal}x_{t}{}^{\intercal}\right)}_{\Delta_{1}+\Delta_{1}{}^{\intercal}\text{
in Lemma \ref{lem:perturbation}}}.\end{split}$ (34)

For matrix
$\bar{\Sigma}_{C}=\Phi_{C}\Sigma_{\alpha}\Phi_{C}{}^{\intercal}+TI=\Phi_{C}\left(\Sigma_{\alpha}+TI\right)\Phi_{C}{}^{\intercal}+T\Phi_{C}^{\perp}\Phi_{C}^{\perp}{}^{\intercal}$,
it is clear that its SVD can be written as

$\begin{split}\bar{\Sigma}_{C}=\begin{bmatrix}\Phi_{C}&\Phi_{C}^{\perp}\end{bmatrix}\begin{bmatrix}\Lambda_{1}&\\\
&TI\end{bmatrix}\begin{bmatrix}\Phi_{C}{}^{\intercal}\\\
\Phi_{C}^{\perp}{}^{\intercal}\end{bmatrix},\quad\Lambda_{1}=\text{diag}\left(\sigma_{1}(\Sigma_{\alpha})+T,\dots,\sigma_{\min}(\Sigma_{\alpha})+T\right).\end{split}$
(35)

where $\Phi_{C}$ is an orthonormal basis of $\text{col}(C)$. Then we conclude
that the following holds for some large enough positive constant $c_{3}$

$\begin{split}{}&\sigma_{1}\left(\bar{\Sigma}_{C}\right)\leq
T+\sigma_{1}(\Sigma_{\alpha})\leq
c_{3}\left(1+\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)T\leq
2c_{3}\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)T,\\\
{}&\sigma_{r_{c}}\left(\bar{\Sigma}_{C}\right)-\sigma_{r_{c}+1}(\bar{\Sigma}_{C})=\sigma_{\min}(\Sigma_{\alpha})\geq\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{32}T,\quad\sigma_{\min}(\bar{\Sigma}_{C})=T.\end{split}$
(36)

Now we are ready to apply Lemma A.7 on $\bar{\Sigma}_{C}$ with the following
constants for $c_{4},c_{5}$ large enough

$\begin{split}{}&\alpha=c_{4}\tilde{U}\sqrt{\psi_{\eta}n\log\frac{1}{\delta}},\quad\beta=c_{5}\psi_{\eta}\sqrt{nT\log\frac{1}{\delta}},\\\
{}&\delta_{M}=\frac{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}{32}T,\quad\psi_{M}=2c_{3}\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}rT\log\frac{1}{\delta},\quad\phi_{M}=T,\quad\sigma_{1}(\Lambda_{2})=T.\end{split}$
(37)

Again,
$\tilde{U}=\sqrt{\log\left(\frac{\psi_{C}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)}$.
From Equation 25, it is clear that $\sqrt{\phi_{M}}\gtrsim\alpha$ and
$\delta_{M}\gtrsim\alpha\sqrt{\psi_{M}}+\beta$. Therefore,

$\begin{split}\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\Phi_{C}\right\|=\left\|\widehat{\Phi}_{C}^{\perp}{}^{\intercal}\tilde{\Phi}_{C}\right\|\lesssim{}&\frac{\psi_{\eta}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\frac{n}{T}\log\frac{1}{\delta}}+\tilde{U}\frac{\sqrt{\psi_{\eta}}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\frac{n}{T}\log\frac{1}{\delta}}+\tilde{U}^{2}\frac{\psi_{\eta}\sqrt{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}\sqrt{1-\rho_{A}^{2}}}\frac{n\sqrt{r}}{T}\log^{2}\frac{1}{\delta}\\\
\lesssim{}&\frac{\psi_{\eta}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\frac{n}{T}\log\frac{1}{\delta}}+\tilde{U}^{2}\frac{\psi_{\eta}\sqrt{\psi_{C}^{2}\psi_{w}\psi_{A}^{2}}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}\sqrt{1-\rho_{A}^{2}}}\frac{n\sqrt{r}}{T}\log^{2}\frac{1}{\delta}\\\
\lesssim{}&\frac{\psi_{\eta}}{\phi_{u}\phi_{C}^{2}\phi_{R}^{2}}\sqrt{\frac{n}{T}\log\frac{1}{\delta}}.\end{split}$
(38)

∎

####  A.1.1 Supporting Details

######  Lemma A.3.

Consider the same setting as Lemma A.2. Consider full row rank matrices
$\\{\alpha\in\mathbb{R}^{a\times r}\\}_{k\in[K]}$ for any positive integer
$a\leq r$. Let

$\begin{split}\Sigma_{\alpha}=\sum_{t\in[T]}\alpha
x_{t}x_{t}{}^{\intercal}\alpha{}^{\intercal},\quad\phi_{\alpha}=\min\\{\sigma_{\min}\left(\alpha\right),1\\},\quad\psi_{\alpha}=\max\\{\sigma_{1}\left(\alpha\right),1\\}.\end{split}$
(39)

If
$T\gtrsim\frac{\psi_{w}^{2}\psi_{\alpha}^{4}\psi_{A}^{4}}{\phi_{u}^{2}\phi_{\alpha}^{4}\phi_{R}^{4}}r^{3}\log\frac{r}{\delta}\cdot\log\left(\frac{\psi_{\alpha}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)$,
then the following holds with probability at least $1-\delta$

$\begin{split}\Sigma_{\alpha}\succeq\frac{\phi_{u}\phi_{\alpha}^{2}\phi_{R}^{2}}{8}TI.\end{split}$
(40)

###### Proof.

For further analysis, we let $\tilde{w}_{t}\coloneqq
Bu_{t}+w_{t}\sim\mathcal{N}\left(0,\Sigma_{\tilde{w}}\right)$ with
$\Sigma_{\tilde{w}}\coloneqq B\Sigma_{u}B{}^{\intercal}+\Sigma_{w}$ and we can
rewrite the system dynamics as follows

$\begin{split}x_{t+1}={}&Ax_{t}+Bu_{t}+w_{t}=Ax_{t}+\tilde{w}_{t},\\\
y_{t}={}&Cx_{t}+\eta_{t}.\end{split}$ (41)

For simplicity, we define $\Sigma_{\alpha,\tau}\coloneqq\alpha
A^{\tau}\left(\sum_{t=1}^{T-\tau}x_{\tau}x_{\tau}{}^{\intercal}\right)(\alpha
A^{\tau}){}^{\intercal}$. It is then clear that

$\begin{split}\Sigma_{\alpha}={}&\Sigma_{\alpha,0}=\sum_{t=1}^{T}\alpha
x_{t}x_{t}{}^{\intercal}\alpha{}^{\intercal}\\\
={}&\sum_{t=1}^{T}\alpha\left(Ax_{t-1}x_{t-1}{}^{\intercal}A{}^{\intercal}+\tilde{w}_{t-1}\tilde{w}_{t-1}{}^{\intercal}+Ax_{t-1}\tilde{w}_{t-1}{}^{\intercal}+\tilde{w}_{t-1}x_{t-1}{}^{\intercal}A{}^{\intercal}\right)\alpha{}^{\intercal}\\\
={}&\Sigma_{\alpha,1}+\sum_{t=0}^{T-1}\alpha\tilde{w}_{t}\tilde{w}_{t}{}^{\intercal}\alpha{}^{\intercal}+\sum_{t=0}^{T-1}\left(\alpha
Ax_{t}\left(\alpha\tilde{w}_{t}\right){}^{\intercal}+\left(\alpha\tilde{w}_{t}\right)x_{t}{}^{\intercal}A{}^{\intercal}\alpha{}^{\intercal}\right).\end{split}$

By Lemma A.4 and A.5, the following events hold with probability at least
$1-\delta/r$,

$\begin{split}\left\|\sum_{t=0}^{T-1}\alpha\tilde{w}_{t}\tilde{w}_{t}{}^{\intercal}\alpha{}^{\intercal}-T\cdot\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}\right\|\lesssim{}&\psi_{w}\psi_{\alpha}^{2}\sqrt{rT\log\frac{2r}{\delta}},\\\
\left\|\left(\Sigma_{\alpha,1}+TI\right)^{-\frac{1}{2}}\sum_{t=0}^{T-1}\alpha
Ax_{t}\left(\alpha\tilde{w}_{t}\right){}^{\intercal}\right\|\lesssim{}&\tilde{U}\sqrt{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}},\end{split}$
(42)

with
$\tilde{U}=\sqrt{\log\left(\frac{\psi_{\alpha}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{2r}{\delta}\right)}$.
From the first inequality, it is clear that the following holds for some large
enough positive constant $c_{1}$

$\begin{split}\sum_{t=0}^{T-1}\alpha\tilde{w}_{t}\tilde{w}_{t}{}^{\intercal}\alpha{}^{\intercal}\succeq
T\cdot\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}-c_{1}\psi_{w}\psi_{\alpha}^{2}\sqrt{rT\log\frac{2r}{\delta}}I.\end{split}$
(43)

From the second inequality, we have the following

$\begin{split}{}&\left\|\left(\Sigma_{\alpha,1}+TI\right)^{-\frac{1}{2}}\sum_{t=0}^{T-1}\alpha
Ax_{t}\left(\alpha\tilde{w}_{t}\right){}^{\intercal}\left(\Sigma_{\alpha,1}+TI\right)^{-\frac{1}{2}}\right\|\\\
\lesssim{}&\tilde{U}\sqrt{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}\left\|\left(\Sigma_{\alpha,1}+TI\right)^{-\frac{1}{2}}\right\|\\\
\lesssim{}&\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}.\end{split}$
(44)

This implies that the following holds for large enough positive constant
$c_{2}$

$\begin{split}{}&\sum_{t=0}^{T-1}\left(\alpha
Ax_{t}\left(\alpha\tilde{w}_{t}\right){}^{\intercal}+\left(\alpha\tilde{w}_{t}\right)x_{t}{}^{\intercal}A{}^{\intercal}\alpha{}^{\intercal}\right)\succeq-
c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\left(\Sigma_{\alpha,1}+TI\right).\end{split}$
(45)

Plugging back into Equation A.1.1 gives the following for some positive
constant $c_{3}$

$\begin{split}\Sigma_{\alpha,0}\succeq{}&\Sigma_{\alpha,1}+T\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}-c_{1}\psi_{w}\psi_{\alpha}^{2}\sqrt{rT\log\frac{2r}{\delta}}I-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\left(\Sigma_{\alpha,1}+TI\right)\\\
\succeq{}&\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\right)\Sigma_{\alpha,1}+T\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}\\\
{}&\hskip 20.00003pt-
c_{1}\psi_{w}\psi_{\alpha}^{2}\sqrt{rT\log\frac{2r}{\delta}}I-c_{2}\tilde{U}\sqrt{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}T}I\\\
\succeq{}&\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\right)\Sigma_{\alpha,1}+T\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\tilde{U}\sqrt{rT\log\frac{2r}{\delta}}I\\\
\end{split}$ (46)

Similarly, we expand
$\Sigma_{\alpha,1},\Sigma_{\alpha,2},\cdots,\Sigma_{\alpha,r-1}$ and have the
following with probability at least $1-\delta$

$\begin{split}\Sigma_{\alpha,0}\succeq{}&\left(T\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\tilde{U}\sqrt{rT\log\frac{2r}{\delta}}I\right)+\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\right)\Sigma_{\alpha,1}\\\
\succeq{}&\left(T\alpha\Sigma_{\tilde{w}}\alpha{}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\tilde{U}\sqrt{rT\log\frac{2r}{\delta}}I\right)\\\
{}&\hskip 20.00003pt+\left((T-1)(\alpha A)\Sigma_{\tilde{w}}(\alpha
A){}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\sqrt{r(T-1)\log\frac{2r}{\delta}}I\right)\cdot\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\right)\\\
{}&\hskip
20.00003pt+\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\log\frac{2r}{\delta}}{T}}\right)\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\log\frac{2r}{\delta}}{T-1}}\right)\Sigma_{\alpha,2}\\\
\succeq{}&\cdots\succeq\left((T-r+1)\sum_{i=0}^{r-1}\left(\alpha
A^{i}\right)\Sigma_{\tilde{w}}\left(\alpha
A^{i}\right){}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\cdot
r\sqrt{rT\log\frac{2r}{\delta}}I\right)\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\log\frac{2r}{\delta}}{T-r+1}}\right)^{r-1}\\\
{}&\hskip
40.00006pt+\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\log\frac{2r}{\delta}}{T-r+1}}\right)^{r}\Sigma_{\alpha,r}\\\
\end{split}$

The above inequality is further simplified as follows

$\begin{split}\Sigma_{\alpha,0}\overset{(i)}{\succeq}{}&\left(\phi_{u}(T-r+1)\sum_{i=0}^{r-1}\left(\alpha
A^{i}\right)BB{}^{\intercal}\left(\alpha
A^{i}\right){}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\cdot
r\sqrt{rT\log\frac{2r}{\delta}}I\right)\left(1-c_{2}\tilde{U}\sqrt{\frac{r\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\log\frac{2r}{\delta}}{T-r+1}}\right)^{r}\\\
\overset{(ii)}{\succeq}{}&\left(\phi_{u}(T-r+1)\alpha\left(\sum_{i=0}^{r-1}A^{i}B(A^{i}B){}^{\intercal}\right)\alpha{}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\cdot
r\sqrt{rT\log\frac{2r}{\delta}}I\right)\left(1-\frac{1}{2r}\right)^{r}\\\
\overset{(iii)}{\succeq}{}&\frac{1}{2}\left(\phi_{u}(T-r+1)\alpha\left(\sum_{i=0}^{r-1}A^{i}B(A^{i}B){}^{\intercal}\right)\alpha{}^{\intercal}-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\cdot
r\sqrt{rT\log\frac{2r}{\delta}}I\right)\\\
\succeq{}&\frac{1}{2}\left(\phi_{u}\phi_{\alpha}^{2}\phi_{R}^{2}(T-r+1)-c_{3}\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}\tilde{U}\cdot
r\sqrt{rT\log\frac{2r}{\delta}}\right)I\\\
\overset{(iv)}{\succeq}{}&\frac{\phi_{u}\phi_{\alpha}^{2}\phi_{R}^{2}}{8}TI\\\
\end{split}$

Here $(i)$ is because $\Sigma_{\alpha,r}\succeq 0$ and
$\Sigma_{\tilde{w}}\succeq
B\Sigma_{u}B{}^{\intercal}\succeq\phi_{u}BB{}^{\intercal}$, $(ii)$ is because
$T\gtrsim\psi_{w}\psi_{\alpha}^{2}\psi_{A}^{2}r^{3}\log\frac{r}{\delta}\cdot\log\left(\frac{\psi_{\alpha}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)$,
$(iii)$ is because $(1-\frac{1}{2r})^{r}\geq\frac{1}{2}$ for all positive
integers and $(iv)$ is because
$T\gtrsim\frac{\psi_{w}^{2}\psi_{\alpha}^{4}\psi_{A}^{4}}{\phi_{u}^{2}\phi_{\alpha}^{4}\phi_{R}^{4}}r^{3}\log\frac{r}{\delta}\cdot\log\left(\frac{\psi_{\alpha}^{2}\psi_{w}\psi_{A}^{4}}{1-\rho_{A}^{2}}r\log\frac{r}{\delta}\right)$.
∎

######  Lemma A.4.

Consider the same setting as Lemma A.2. For any positive integer $a$, let
$\\{\zeta_{t}\in\mathbb{R}^{a}\\}_{t=0}^{T}$ be a sequence of i.i.d Gaussian
vectors from $\mathcal{N}(0,\Sigma_{\zeta})$ such that $\zeta_{t}$ is
independent of $x_{t}$. Define $\psi_{\zeta}=\sigma_{1}(\Sigma_{\zeta})$. We
make the following defition for any matrix $P\in\mathbb{R}^{b\times r}$ for
positive integer $b$

$\begin{split}\bar{\Sigma}_{P}=\sum_{t\in[T]}Px_{t}x_{t}{}^{\intercal}P{}^{\intercal}+TI,\quad\psi_{P}=\sigma_{1}(P).\end{split}$
(47)

Then for any $\delta$, the following holds with probability at least
$1-\delta$,

$\begin{split}\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}{}^{\intercal}\right\|\lesssim\tilde{U}\sqrt{\max\\{r,a\\}\psi_{\zeta}\log\frac{1}{\delta}}.\end{split}$
(48)

Here
$\tilde{U}=\sqrt{\log\left(\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)}$.

###### Proof.

From Lemma A.9, we know that the following holds for some vector
$v\in\mathcal{S}^{a-1}$

$\begin{split}\mathbb{P}\left(\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}{}^{\intercal}\right\|>z\right)\leq
5^{a}\mathbb{P}\left(\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}{}^{\intercal}v\right\|>\frac{z}{2}\right).\end{split}$
(49)

Notice that $\zeta_{t}{}^{\intercal}v$ are independent Gaussian variables from
distribution $\mathcal{N}(0,v{}^{\intercal}\Sigma_{\zeta}v)$ for all
$t\in[T]$, which is $c_{1}\sqrt{\psi_{\zeta}}$-subGaussian for some positive
constant $c_{1}$. Then applying Theorem 1 in (Abbasi-yadkori et al., 2011) on
sequence $\\{Px_{t}\\}_{t\in[T]}$ and sequence
$\\{\zeta_{t}{}^{\intercal}v\\}_{t\in[T]}$, gives the following inequality

$\begin{split}\mathbb{P}\left(\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}{}^{\intercal}v\right\|>\sqrt{2c_{1}^{2}\psi_{\zeta}\log\left(\frac{\det(\bar{\Sigma}_{P})^{\frac{1}{2}}\det(TI)^{-\frac{1}{2}}}{\delta}\right)}\right)\leq\delta.\end{split}$
(50)

Substituting the above result back gives the following inequality

$\begin{split}\mathbb{P}\left(\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}\right\|>\sqrt{8c_{1}^{2}\psi_{\zeta}\log\left(\frac{\det(\bar{\Sigma}_{P})^{\frac{1}{2}}\det(TI)^{-\frac{1}{2}}}{\delta}\right)}\right)\leq
5^{a}\delta,\end{split}$ (51)

which implies the following inequality holds with probability at least
$1-\delta/2$

$\begin{split}\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}\right\|\leq\sqrt{8c_{1}^{2}\psi_{\zeta}\log\left(\frac{2\det(\bar{\Sigma}_{P})^{\frac{1}{2}}\det(TI)^{-\frac{1}{2}}}{\delta}\right)+8c_{1}^{2}\psi_{\zeta}a\log
5}.\end{split}$ (52)

Now consider
$\bar{\Sigma}_{P}=\sum_{t\in[T]}Px_{t}x_{t}{}^{\intercal}P{}^{\intercal}+TI$.
Then

$\begin{split}\det(TI)={}&T^{b},\quad\det(\bar{\Sigma}_{P})=T^{b-r}\prod_{i=1}^{r}\left(T+\lambda_{i}\left(\sum_{t\in[T]}Px_{t}x_{t}{}^{\intercal}P{}^{\intercal}\right)\right).\end{split}$
(53)

From Lemma A.6, we have the following with probability at least $1-\delta/2$

$\begin{split}\left\|\sum_{t\in[T]}x_{t}x_{t}{}^{\intercal}\right\|\lesssim\frac{\psi_{w}\psi_{A}^{2}r}{1-\rho_{A}^{2}}T\log\frac{2}{\delta}\lesssim\frac{\psi_{w}\psi_{A}^{2}r}{1-\rho_{A}^{2}}T\log\frac{1}{\delta}.\end{split}$
(54)

Therefore,

$\begin{split}\lambda_{i}\left(\sum_{t\in[T]}Px_{t}x_{t}{}^{\intercal}P{}^{\intercal}\right)\leq{}&\left\|\sum_{t\in[T]}x_{t}x_{t}{}^{\intercal}\right\|\left\|PP{}^{\intercal}\right\|\\\
\lesssim{}&\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}rT\log\frac{1}{\delta}\\\
\end{split}$ (55)

Substituting back gives the following for some positive constant $c_{2}$

$\begin{split}\det(\bar{\Sigma}_{P})={}&T^{b-r}\prod_{i=1}^{r}\left(T+\lambda_{i}\left(\sum_{t\in[T]}Px_{t}x_{t}{}^{\intercal}P{}^{\intercal}\right)\right)\\\
\leq{}&T^{b}\left(1+c_{2}\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)^{r},\end{split}$
(56)

which gives

$\begin{split}\det(\bar{\Sigma}_{P})^{\frac{1}{2}}\det(TI)^{-\frac{1}{2}}\leq\left(1+c_{2}\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)^{\frac{r}{2}}.\end{split}$
(57)

Finally, with a union bound over above events, we get the following with
probability at least $1-\delta$ from Equation 52 and 56

$\begin{split}{}&\left\|\left(\bar{\Sigma}_{P}\right)^{-\frac{1}{2}}\sum_{t\in[T]}Px_{t}\zeta_{t}\right\|\\\
\leq{}&\sqrt{8c_{1}^{2}\psi_{\zeta}\log\left(\frac{2\det(\bar{\Sigma}_{P})^{\frac{1}{2}}\det(TI)^{-\frac{1}{2}}}{\delta}\right)+8c_{1}^{2}\psi_{\zeta}a\log
5}\\\
\leq{}&\sqrt{4c_{1}^{2}r\psi_{\zeta}\log\left(2+2c_{2}\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)+8c_{1}^{2}\psi_{\zeta}\log\frac{1}{\delta}+8c_{1}^{2}\psi_{\zeta}a\log
5}\\\
\lesssim{}&\sqrt{\max\\{r,a\\}\psi_{\zeta}\log\frac{1}{\delta}}\cdot\sqrt{\log\left(\frac{\psi_{P}^{2}\psi_{w}\psi_{A}^{2}}{1-\rho_{A}^{2}}r\log\frac{1}{\delta}\right)}.\end{split}$
(58)

This completes the proof. ∎

######  Lemma A.5.

For any positive integer $a$, let $\\{w_{t}\in\mathbb{R}^{a}\\}_{t\in[T]}$ be
independent Gaussian vectors from $w_{t}\sim\mathcal{N}(0,\Sigma_{w,t})$. Let

$\begin{split}\psi_{w}=\max_{t\in[T]}\sigma_{1}(\Sigma_{w,t}),\quad\phi_{w}=\min_{t\in[T]}\sigma_{\min}(\Sigma_{w,t}).\end{split}$
(59)

Then with probability at least $1-\delta$,

$\begin{split}\left\|\sum_{t=1}^{T}w_{t}w_{t}{}^{\intercal}-\sum_{t=1}^{T}\Sigma_{w,t}\right\|\lesssim\psi_{w}\sqrt{aT\log\frac{1}{\delta}}.\end{split}$
(60)

###### Proof.

From Lemma A.9, we know that there exist $x,y\in\mathcal{S}^{a-1}$ s.t.

$\begin{split}\mathbb{P}\left(\left\|\sum_{t=1}^{T}w_{t}w_{t}{}^{\intercal}-\sum_{t=1}^{T}\Sigma_{w,t}\right\|>z\right)\leq
5^{2a}\mathbb{P}\left(\sum_{t=1}^{T}\Big{[}(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\Big{]}>\frac{1}{4}z\right)\end{split}$
(61)

For any $x,y$, following the definition of $w_{t}$, we have

$\begin{split}\|(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\|_{\psi_{1}}\lesssim\left\|(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})\right\|_{\psi_{1}}\leq\left\|y{}^{\intercal}w_{t}\right\|_{\psi_{2}}\left\|x{}^{\intercal}w_{t}\right\|_{\psi_{2}}\lesssim\psi_{w}.\end{split}$
(62)

Here the first and second inequalities are from Lemma 2.6.8 and Lemma 2.7.7 in
(Vershynin, 2018). Directly applying Bernstein’s inequality on
$\\{(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\\}_{t=1}^{T}$
(subexponential random variables) gives the following for some positive
constant $c_{1}$ with probability at least $1-\delta$

$\begin{split}\left\|\sum_{t=1}^{T}\Big{[}(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\Big{]}\right\|\leq
c_{1}\psi_{w}\sqrt{T\log(\frac{2}{\delta})},\end{split}$ (63)

which is equivalent to

$\begin{split}5^{2a}\mathbb{P}\left(\left\|\sum_{t=1}^{T}\Big{[}(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\Big{]}\right\|>c_{1}\psi_{w}\sqrt{T\log(\frac{2}{\delta})}\right)\leq
5^{2a}\delta.\end{split}$ (64)

Letting $\delta^{\prime}=5^{2a}\delta$ gives

$\begin{split}5^{2a}\mathbb{P}\left(\left\|\sum_{t=1}^{T}\Big{[}(y{}^{\intercal}w_{t})(x{}^{\intercal}w_{t})-y{}^{\intercal}\Sigma_{w,t}x\Big{]}\right\|>c_{1}\psi_{w}\sqrt{2aT\log(\frac{10}{\delta^{\prime}})}\right)\leq\delta^{\prime}.\end{split}$
(65)

Substituting back into Equation 61 gives

$\begin{split}\mathbb{P}\left(\left\|\sum_{t=1}^{T}w_{t}w_{t}{}^{\intercal}-\sum_{t=1}^{T}\Sigma_{w,t}\right\|>4c_{1}\psi_{w}\sqrt{2aT\log(\frac{10}{\delta^{\prime}})}\right)\leq\delta^{\prime}.\end{split}$
(66)

This finishes the first part of the proof.

∎

######  Lemma A.6.

Consider the same setting as Lemma A.2. For any $\delta$, with probability at
least $1-\delta$,

$\begin{split}\left\|\sum_{t=0}^{T}x_{t}x_{t}{}^{\intercal}\right\|\lesssim\frac{\psi_{w}\psi_{A}^{2}r}{1-\rho_{A}^{2}}T\log\frac{1}{\delta}.\end{split}$
(67)

###### Proof.

Directly applying Proposition 8.4 in (Sarkar & Rakhlin, 2019) gives the
following with probability at least $1-\delta$

$\begin{split}\left\|\sum_{t=0}^{T}x_{t}x_{t}{}^{\intercal}\right\|\lesssim{}&\psi_{w}\text{tr}\left(\sum_{t=0}^{T-1}\Gamma_{t}(A)\right)\log\frac{1}{\delta}\leq\psi_{w}r\left\|\sum_{t=0}^{T-1}\Gamma_{t}(A)\right\|\log\frac{1}{\delta}.\end{split}$
(68)

From Assumption 3.1, we know that

$\begin{split}\left\|\Gamma_{t}(A)\right\|=\left\|\sum_{\tau=0}^{t}A^{\tau}A^{\tau}{}^{\intercal}\right\|\preceq\sum_{\tau=0}^{t}\left\|A^{\tau}\right\|^{2}\leq
1+\frac{\psi_{A}^{2}}{1-\rho_{A}^{2}}\lesssim\frac{\psi_{A}^{2}}{1-\rho_{A}^{2}}.\end{split}$
(69)

Substituting back gives

$\begin{split}\left\|\sum_{t=0}^{T}x_{t}x_{t}{}^{\intercal}\right\|\lesssim\frac{\psi_{w}\psi_{A}^{2}r}{1-\rho_{A}^{2}}T\log\frac{1}{\delta}.\end{split}$
(70)

∎

####  A.1.2 Other Lemmas

######  Lemma A.7.

Given perturbation matrices $\Delta_{1},\Delta_{2}$ and p.d square matrix $M$
of the same dimension. Suppose they satisfy the following for some positive
constant $\alpha\in[0,1)$ and $\beta\geq 0$

$\begin{split}\left\|M^{-\frac{1}{2}}\Delta_{1}\right\|\leq\alpha<\sqrt{\phi_{M}},\quad\left\|\Delta_{2}\right\|\leq\beta.\end{split}$
(71)

Denote the SVD of the matrices as follows

$\begin{split}M=\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}\begin{bmatrix}\Lambda_{1}&\\\
&\Lambda_{2}\end{bmatrix}\begin{bmatrix}U_{1}{}^{\intercal}\\\
U_{2}{}^{\intercal}\end{bmatrix},\quad\tilde{M}\coloneqq
M+\left(\Delta_{1}+\Delta_{1}{}^{\intercal}\right)+\Delta_{2}=\begin{bmatrix}\tilde{U}_{1}&\tilde{U}_{2}\end{bmatrix}\begin{bmatrix}\tilde{\Lambda}_{1}&\\\
&\tilde{\Lambda}_{2}\end{bmatrix}\begin{bmatrix}\tilde{U}_{1}{}^{\intercal}\\\
\tilde{U}_{2}{}^{\intercal}\end{bmatrix},\end{split}$ (72)

and define the following constants

$\begin{split}\psi_{M}=\sigma_{1}(M),\quad\phi_{M}=\sigma_{\min}(M),\quad
0<\delta_{M}\leq\sigma_{\min}(\Lambda_{1})-\sigma_{1}(\Lambda_{2}).\end{split}$
(73)

Then we have

$\begin{split}\left\|\tilde{U}_{2}{}^{\intercal}U_{1}\right\|\leq\frac{\alpha^{2}+\beta}{\delta_{M}-4\alpha\sqrt{\psi_{M}}-3\alpha^{2}-\beta}+\frac{\alpha}{\sqrt{\phi_{M}}}\left(1+\frac{2-\alpha/\sqrt{\phi_{M}}}{1-\alpha/\sqrt{\phi_{M}}}\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\delta_{M}-2\alpha\sqrt{\psi_{M}}-\alpha^{2}}\right)\end{split}$
(74)

In the case where $\sqrt{\phi_{M}}\geq C\alpha$, $\delta_{M}\geq
C\left(\alpha\sqrt{\psi_{M}}+\beta\right)$ for large enough $C$, the above
result is further simplified to

$\begin{split}\left\|\tilde{U}_{2}{}^{\intercal}U_{1}\right\|\lesssim\frac{\beta}{\delta_{M}}+\frac{\alpha\sigma_{1}(\Lambda_{2})/\sqrt{\phi_{M}}}{\delta_{M}}+\frac{\alpha^{2}\sqrt{\psi_{M}/\phi_{M}}}{\delta_{M}}.\end{split}$
(75)

######  Remark A.8.

This result is adapted from previous subspace relative perturbation results
(Theorem 3.2 in (Li, 1998)). It provides tighter bound for matrices with
$\phi_{M}$ is close to $\psi_{M}$, as compared to the standard Davis-Kahan
$\sin\Theta$ theorem. To be more specific, Davis-Kahan gives the bound of
order $\frac{\alpha\sqrt{\psi_{M}}}{\delta_{M}}$. Here our result features the
bound $\frac{\alpha\sigma_{1}(\Lambda_{2})/\sqrt{\phi_{M}}}{\delta_{M}}$.

###### Proof.

Define

$\begin{split}\widehat{M}=M+\left(\Delta_{1}+\Delta_{1}{}^{\intercal}\right)+\Delta_{1}{}^{\intercal}M^{-1}\Delta_{1}=\begin{bmatrix}\widehat{U}_{1}&\widehat{U}_{2}\end{bmatrix}\begin{bmatrix}\widehat{\Lambda}_{1}&\\\
&\widehat{\Lambda}_{2}\end{bmatrix}\begin{bmatrix}\widehat{U}_{1}{}^{\intercal}\\\
\widehat{U}_{2}{}^{\intercal}\end{bmatrix}.\end{split}$ (76)

From its definition, we immediately have

$\begin{split}\left\|\widehat{M}-M\right\|\leq{}&2\left\|\Delta_{1}\right\|+\left\|\Delta_{1}{}^{\intercal}M^{-1}\Delta_{1}\right\|\\\
\overset{(i)}{\leq}{}&2\alpha\sqrt{\psi_{M}}+\left\|M^{-\frac{1}{2}}\Delta_{1}\right\|^{2}\\\
\leq{}&2\alpha\sqrt{\psi_{M}}+\alpha^{2}.\end{split}$ (77)

Here the first term in $(i)$ is because

$\begin{split}\alpha\geq\left\|M^{-\frac{1}{2}}\Delta_{1}\right\|\geq\sigma_{\min}\left(M^{-\frac{1}{2}}\right)\left\|\Delta_{1}\right\|=\frac{1}{\sqrt{\psi_{M}}}\left\|\Delta_{1}\right\|.\end{split}$
(78)

And the second term is because
$\left\|A{}^{\intercal}A\right\|=\left\|A\right\|^{2}$. Then from Theorem 1 in
(Stewart, 1990), we know that

$\begin{split}\sigma_{1}\left(\widehat{\Lambda}_{2}\right)\leq\sigma_{1}\left(\Lambda_{2}\right)+\left\|\widehat{M}-M\right\|\leq\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}.\end{split}$
(79)

which in turn gives

$\begin{split}\widehat{\delta}\coloneqq\sigma_{\min}(\Lambda_{1})-\sigma_{1}(\widehat{\Lambda}_{2})\geq\delta_{M}-2\alpha\sqrt{\psi_{M}}-\alpha^{2}.\end{split}$
(80)

For the rest of the proof, we upper bound
$\left\|\widehat{U}_{2}{}^{\intercal}U_{1}\right\|$,
$\left\|\tilde{U}_{2}{}^{\intercal}\widehat{U}_{1}\right\|$ and then derives
the desired result. To upper bound
$\left\|\widehat{U}_{2}{}^{\intercal}U_{1}\right\|$, we let
$D=M^{-1}\Delta_{1}$ with

$\begin{split}{}&\left\|D\right\|\leq\left\|M^{-\frac{1}{2}}\right\|\left\|M^{-\frac{1}{2}}\Delta_{1}\right\|\leq\frac{\alpha}{\sqrt{\phi_{M}}}<1,\\\
{}&\left\|\left(I+D\right)^{-1}\right\|=\frac{1}{\sigma_{\min}\left(I+D\right)}\overset{(i)}{\leq}\frac{1}{1-\sigma_{1}(D)}\leq\frac{1}{1-\alpha/\sqrt{\phi_{M}}}.\end{split}$
(81)

Here $(i)$ is because $|\sigma_{\min}(I+D)-\sigma_{\min}(I)|\leq\sigma_{1}(D)$
from Thereom 1 in (Stewart, 1990). And the definition of $D$ gives

$\begin{split}\widehat{M}=\left(I+\Delta_{1}{}^{\intercal}M^{-1}\right)M\left(I+M^{-1}\Delta_{1}\right)=\left(I+D{}^{\intercal}\right)M\left(I+D\right).\end{split}$
(82)

Now from Equation 80 we apply Theorem 3.2 in (Li, 1998) gives

$\begin{split}\left\|\widehat{U}_{2}{}^{\intercal}U_{1}\right\|\leq{}&\left\|I-(I+D{}^{\intercal})\right\|+\left(\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}\right)\frac{\left\|(I+D){}^{\intercal}-(I+D)^{-1}\right\|}{\widehat{\delta}}\\\
\leq{}&\frac{\alpha}{\sqrt{\phi_{M}}}+\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\widehat{\delta}}\left(\left\|(I+D){}^{\intercal}-I\right\|+\left\|I-\left(I+D\right)^{-1}\right\|\right)\\\
\leq{}&\frac{\alpha}{\sqrt{\phi_{M}}}+\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\widehat{\delta}}\left(\frac{\alpha}{\sqrt{\phi_{M}}}+\left\|\left(I+D\right)^{-1}\right\|\left\|I+D-I\right\|\right)\\\
\leq{}&\frac{\alpha}{\sqrt{\phi_{M}}}+\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\widehat{\delta}}\left(\frac{\alpha}{\sqrt{\phi_{M}}}+\frac{1}{1-\alpha/\sqrt{\phi_{M}}}\frac{\alpha}{\sqrt{\phi_{M}}}\right)\\\
={}&\frac{\alpha}{\sqrt{\phi_{M}}}\left(1+\frac{2-\alpha/\sqrt{\phi_{M}}}{1-\alpha/\sqrt{\phi_{M}}}\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\widehat{\delta}}\right).\end{split}$
(83)

We now upper bound
$\left\|\tilde{U}_{2}{}^{\intercal}\widehat{U}_{1}\right\|$. We start from
upper bounding $\left\|\widehat{M}-\tilde{M}\right\|$

$\begin{split}\left\|\widehat{M}-\tilde{M}\right\|=\left\|\Delta_{1}{}^{\intercal}M^{-1}\Delta_{1}-\Delta_{2}\right\|\leq\left\|\Delta_{1}{}^{\intercal}M^{-1}\Delta_{1}\right\|+\left\|\Delta_{2}\right\|\leq\alpha^{2}+\beta.\end{split}$
(84)

Moreover, we know that

$\begin{split}{}&\sigma_{\min}(\widehat{\Lambda}_{1})-\sigma_{1}(\tilde{\Lambda}_{2})\\\
={}&\sigma_{\min}(\widehat{\Lambda}_{1})-\sigma_{1}(\widehat{\Lambda}_{2})+\sigma_{1}(\widehat{\Lambda}_{2})-\sigma_{1}(\tilde{\Lambda}_{2})\\\
\geq{}&\sigma_{\min}(\widehat{\Lambda}_{1})-\sigma_{1}(\widehat{\Lambda}_{2})-\alpha^{2}-\beta\\\
={}&\sigma_{\min}(\widehat{\Lambda}_{1})-\sigma_{\min}(\Lambda_{1})+\sigma_{\min}(\Lambda_{1})-\sigma_{1}(\Lambda_{2})+\sigma_{1}(\Lambda_{2})-\sigma_{1}(\widehat{\Lambda}_{2})-\alpha^{2}-\beta\\\
\geq{}&-2\left\|\widehat{M}-M\right\|-\alpha^{2}-\beta+\delta_{M}\\\
\geq{}&-4\alpha\sqrt{\psi_{M}}-2\alpha^{2}-\alpha^{2}-\beta+\delta_{M}.\end{split}$
(85)

Here the first and second inequalities are due to Stewart (1990, Theorem 1) .
Directly applying the Davis-Kahan $\sin\theta$ Theorem in (Davis & Kahan,
1970) gives

$\begin{split}\left\|\tilde{U}_{2}{}^{\intercal}\widehat{U}_{1}\right\|\leq\frac{\left\|\widehat{M}-\tilde{M}\right\|}{\sigma_{\min}(\widehat{\Lambda}_{1})-\sigma_{1}(\tilde{\Lambda}_{2})}\leq\frac{\alpha^{2}+\beta}{-4\alpha\sqrt{\psi_{M}}-3\alpha^{2}-\beta+\delta_{M}}.\end{split}$
(86)

Finally, we conclude that

$\begin{split}\left\|\tilde{U}_{2}{}^{\intercal}U_{1}\right\|={}&\left\|\tilde{U}_{2}{}^{\intercal}\left(\widehat{U}_{1}\widehat{U}_{1}{}^{\intercal}+\widehat{U}_{2}\widehat{U}_{2}{}^{\intercal}\right)U_{1}\right\|\leq\left\|\tilde{U}_{2}{}^{\intercal}\widehat{U}_{1}\right\|+\left\|\widehat{U}_{2}{}^{\intercal}U_{1}\right\|\\\
\leq{}&\frac{\alpha^{2}+\beta}{\delta_{M}-4\alpha\sqrt{\psi_{M}}-3\alpha^{2}-\beta}+\frac{\alpha}{\sqrt{\phi_{M}}}\left(1+\frac{2-\alpha/\sqrt{\phi_{M}}}{1-\alpha/\sqrt{\phi_{M}}}\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\widehat{\delta}}\right)\\\
\leq{}&\frac{\alpha^{2}+\beta}{\delta_{M}-4\alpha\sqrt{\psi_{M}}-3\alpha^{2}-\beta}+\frac{\alpha}{\sqrt{\phi_{M}}}\left(1+\frac{2-\alpha/\sqrt{\phi_{M}}}{1-\alpha/\sqrt{\phi_{M}}}\frac{\sigma_{1}\left(\Lambda_{2}\right)+2\alpha\sqrt{\psi_{M}}+\alpha^{2}}{\delta_{M}-2\alpha\sqrt{\psi_{M}}-\alpha^{2}}\right).\end{split}$
(87)

∎

######  Lemma A.9.

For any positive integer $a,b$, let $M\in\mathbb{R}^{a\times b}$ be any random
matrix. There exist $x\in\mathcal{S}^{b-1}$ and $y\in\mathcal{S}^{a-1}$ such
that for all $\varepsilon<1$

$\begin{split}\mathbb{P}(\left\|M\right\|>z)\leq\left(1+\frac{2}{\varepsilon}\right)^{a+b}\mathbb{P}\left(\left\|y{}^{\intercal}Mx\right\|>(1-\varepsilon)^{2}z\right).\end{split}$
(88)

###### Proof.

From Proposition 8.1 of (Sarkar & Rakhlin, 2019), we know that there exist an
$x\in\mathcal{S}^{b-1}$ s.t. the following holds for any $\varepsilon<1$

$\begin{split}\mathbb{P}(\left\|M\right\|>z)\leq\left(1+\frac{2}{\varepsilon}\right)^{b}\mathbb{P}\left(\left\|Mx\right\|>(1-\varepsilon)z\right)=\left(1+\frac{2}{\varepsilon}\right)^{b}\mathbb{P}\left(\left\|x{}^{\intercal}M{}^{\intercal}\right\|>(1-\varepsilon)z\right).\end{split}$
(89)

Applying the proposition again on $x{}^{\intercal}M{}^{\intercal}$ gives the
following for some vector $y\in\mathcal{S}^{a-1}$

$\begin{split}\mathbb{P}\left(\left\|x{}^{\intercal}M{}^{\intercal}\right\|>(1-\varepsilon)z\right)\leq\left(1+\frac{2}{\varepsilon}\right)^{a}\mathbb{P}\left(\left\|x{}^{\intercal}My\right\|>(1-\varepsilon)^{2}z\right).\end{split}$
(90)

Therefore, we have

$\begin{split}\mathbb{P}(\left\|M\right\|>z)\leq\left(1+\frac{2}{\varepsilon}\right)^{a+b}\mathbb{P}\left(\left\|x{}^{\intercal}M{}^{\intercal}y\right\|>(1-\varepsilon)^{2}z\right)\end{split}$
(91)

∎

##  Appendix B Lower Bounds for HD-SYSID — Proof of Theorem B.2

An essential part for the proof is the Birge’s Inequality stated below.

######  Theorem B.1.

Let $\\{\mathbb{P}_{i}\\}_{i\in[N]}$ be probability distributions on
$(\Omega,\mathcal{F})$. Let $\\{F_{i}\\}_{i\in[N]}\in\mathcal{F}$ be pairwise
disjoint events. If
$\delta=\min_{i\in[N]}\mathbb{P}_{i}(F_{i})\geq\frac{1}{2}$ (the success
rate), then

$\begin{split}(2\delta-1)\log\left(\frac{\delta}{1-\delta}(N-1)\right)\leq\frac{1}{N-1}\sum_{i\in[N]}\text{KL}(\mathbb{P}_{i}||\mathbb{P}_{1}).\end{split}$
(92)

###### Proof.

From the Birge’s inequality (Boucheron et al., 2013), we directly have

$\begin{split}\delta\log\left(\delta/\frac{1-\delta}{N-1}\right)+(1-\delta)\log\left((1-\delta)/(1-\frac{1-\delta}{N-1})\right)\leq\frac{1}{N-1}\sum_{i\in[N]}\text{KL}(\mathbb{P}_{i}||\mathbb{P}_{1}),\end{split}$
(93)

On the other hand, the following holds for any $N\geq 3$

$\begin{split}{}&\delta\log\left(\delta/\frac{1-\delta}{N-1}\right)+(1-\delta)\log\left((1-\delta)/(1-\frac{1-\delta}{N-1})\right)\\\
\overset{(i)}{\geq}{}&\delta\log\left(\delta/\frac{1-\delta}{N-1}\right)+(1-\delta)\log\left((1-\delta)/((N-1)\delta)\right)\\\
={}&(2\delta-1)\log\left(\frac{\delta}{1-\delta}(N-1)\right)\end{split}$ (94)

Here $(i)$ holds because $1-\frac{1-\delta}{N-1}\leq 1\leq
2\delta\leq(N-1)\delta$. For $N=2$, the above inequality naturally holds.
Combining the two cases finishes the proof. ∎

Now we prove Theorem 4.1 by proving the following more general version that
holds for multiple trajectories.

######  Theorem B.2.

Suppose $n\geq 2$, $K\geq 1$, and choose positive scalars
$\delta\leq\frac{1}{2}$ and $\epsilon\leq 0.6$. Consider the class of minimal
systems $\mathcal{M}=(r,n,m,A,B,C,\Sigma_{w},\Sigma_{\eta})$ with different
$A$,$B$,$C$ matrices. All other parameters are fixed and known. Moreover,
$r<n$ and $\Sigma_{\eta}$ is positive definite. For every $k\in[K]$, let
$\mathcal{D}_{k}=\\{y_{k,t}\\}_{t=0}^{T_{k}}\cup\\{u_{k,t}\\}_{t=0}^{T_{k}-1}\cup\text{\\{all
known parameters\\}}$ denote the associated single trajectory dataset. Here
the inputs $u_{k,t}$ satisfy: 1). $u_{k,t}$ is sampled independently;2).
$\mathbb{E}(u_{k,t})=0$. Consider any estimator $\widehat{f}$ mapping the
datasets $\mathcal{D}\coloneqq\cup_{k\in[K]}\mathcal{D}_{k}$ to
$(\widehat{A}(\mathcal{D}),\widehat{B}(\mathcal{D}),\widehat{C}(\mathcal{D}))\in\mathbb{R}^{r\times
r}\times\mathbb{R}^{r\times m}\times\mathbb{R}^{n\times r}$. If

$\mathcal{T}\coloneqq\sum_{k\in[K]}T_{k}<\frac{\phi_{\eta}(1-2\delta)\log
1.4}{50\left(\psi_{w}+\psi_{u}\right)}\cdot\frac{n+\log\frac{1}{2\delta}}{\epsilon^{2}},$
(95)

there exists a system
$\mathcal{M}_{0}=(r,n,m,A_{0},B_{0},C_{0},\Sigma_{w},\Sigma_{\eta})$ with
dataset $\mathcal{D}$ such that

$\begin{split}\mathbb{P}{}&\left\\{\left\|\widehat{C}\left(\mathcal{D}\right)\widehat{B}\left(\mathcal{D}\right)-C_{0}B_{0}\right\|\geq\epsilon\right\\}\geq\delta.\end{split}$
(96)

Here $\mathbb{P}$ denotes the distribution of $\mathcal{D}$ generated by
system $\mathcal{M}_{0}$. Related constants are defined as follows

$\begin{split}\phi_{\eta}&\coloneqq\sigma_{\min}(\Sigma_{\eta}),\quad\psi_{w}\coloneqq\sigma_{1}(\Sigma_{w}),\\\
\quad\psi_{u}&\coloneqq\max_{k,t}\sigma_{\max}\left(\mathbb{E}(u_{k,t}u_{k,t}{}^{\intercal})\right).\end{split}$

###### Proof.

Fix a constant $\epsilon<0.1$ and a $\frac{1}{2}$-packing of the intersection
of the unit ball in $\mathbb{R}^{n}$, denoted by $\tilde{\mathcal{P}}$. Now
let
$\mathcal{P}=\tilde{\mathcal{P}}\backslash\\{(-1,0,\dots,0){}^{\intercal}\\}$.
We know $|\mathcal{P}|\geq 2^{n}-1$ and we let
$\mathcal{P}=\\{p_{i}\\}_{i=1}^{|\mathcal{P}|}$. Now consider the following
set of matrices $\\{A_{i},B_{i},C_{i}\\}_{i=1}^{|\mathcal{P}|}$ where

$\begin{split}{}&A_{i}=\epsilon\left(E_{r,1}+\sum_{j\in[r-1]}E_{j,j+1}\right)\in\mathbb{R}^{r\times
r},\quad B_{i}=E_{1,1}\in\mathbb{R}^{r\times m},\\\
{}&C_{i}=4\epsilon\left(p_{i}e_{1}{}^{\intercal}+\sum_{j\in[r]}E_{j,j}\right)\in\mathbb{R}^{n\times
r}.\end{split}$

Step 1: We first prove that all systems are both observable and controllable,
and therefore a minimal realizations. Firstly, notice that
$\text{rank}(C_{i})=r$. Therefore the observability matrix is full column
rank, implying that the system is observable. For controllability, we first
show by induction that for any interger $a$ satisfying $1\leq a<r$, the
following holds

$\begin{split}A_{i}^{a}=\epsilon^{a}\left(\sum_{j=1}^{a}E_{r-a+j,j}+\sum_{j=1}^{r-a}E_{j,j+a}\right).\end{split}$
(97)

The base case where $a=1$ is easily verified from the definition of $A_{i}$.
Now suppose the induction hypothesis holds for all positive integers $a\leq
a_{0}$ for integer $a_{0}\in[1,r-2]$. Then for $a=a_{0}+1<r$, we anow that

$\begin{split}A_{i}^{a_{0}+1}={}&A_{i}^{a_{0}}A_{i}=\epsilon^{a+1}\left(\sum_{j_{0}=1}^{a_{0}}E_{r-a_{0}+j_{0},j_{0}}+\sum_{j_{0}=1}^{r-a_{0}}E_{j_{0},j_{0}+a_{0}}\right)\left(E_{r,1}+\sum_{j\in[r-1]}E_{j,j+1}\right)\\\
={}&\epsilon^{a+1}\left(\sum_{j_{0}=1}^{a_{0}}E_{r-a_{0}+j_{0},j_{0}}E_{r,1}+\sum_{j_{0}=1}^{a_{0}}E_{r-a_{0}+j_{0},j_{0}}\sum_{j\in[r-1]}E_{j,j+1}\right.\\\
{}&\hskip
20.00003pt+\left.\sum_{j_{0}=1}^{r-a_{0}}E_{j_{0},j_{0}+a_{0}}E_{r,1}+\sum_{j_{0}=1}^{r-a_{0}}E_{j_{0},j_{0}+a_{0}}\sum_{j\in[r-1]}E_{j,j+1}\right)\\\
\overset{(i)}{=}{}&\epsilon^{a+1}\left(0+\sum_{j_{0}=1}^{a_{0}}E_{r-a_{0}+j_{0},j_{0}+1}+E_{r-a_{0},1}+\sum_{j_{0}=1}^{r-a_{0}}E_{j_{0},j_{0}+a_{0}+1}\right)\\\
={}&\epsilon^{a+1}\left(\sum_{j_{0}=1}^{a_{0}}E_{r-(a_{0}+1)+(j_{0}+1),j_{0}+1}+E_{r-(a_{0}+1)+1,1}+\sum_{j_{0}=1}^{r-(a_{0}+1)}E_{j_{0},j_{0}+a_{0}+1}\right)\\\
={}&\epsilon^{a+1}\left(\sum_{j_{0}=1}^{a_{0}+1}E_{r-(a_{0}+1)+j_{0},j_{0}}+\sum_{j_{0}=1}^{r-(a_{0}+1)}E_{j_{0},j_{0}+a_{0}+1}\right).\end{split}$
(98)

Here, again, $(i)$ is because $E_{a_{1},b_{1}}E_{a_{2},b_{2}}=E_{a_{1},b_{2}}$
if $b_{1}=a_{2}$. Otherwise the product is zero. This finishes the induction.
Now with $A_{i}^{a}$ for positive integer $a<r$, we directly have

$\begin{split}(A_{i}^{a}B_{i})_{1}={}&\epsilon^{a}\left(\left(\sum_{j=1}^{a}E_{r-a+j,j}+\sum_{j=1}^{r-a}E_{j,j+a}\right)\sum_{j\in\min\\{r,m\\}}E_{j,j}\right)_{1}\\\
={}&\epsilon^{a}\left(\left(\sum_{j=1}^{a}E_{r-a+j,j}+\sum_{j=1}^{r-a}E_{j,j+a}\right)E_{1,1}\right)_{1}\\\
={}&\epsilon^{a}E_{r-a+1,1}.\end{split}$ (99)

Here for any matrix $M$, we use $(M)_{1}$ to denote its first column.
Therefore, we have

$\begin{split}{}&\text{rank}\left(\begin{bmatrix}B&AB&\cdots&A^{r-1}B\end{bmatrix}\right)\\\
\geq{}&\text{rank}\left(\begin{bmatrix}(B)_{1}&(AB)_{1}&\cdots&(A^{r-1}B)_{1}\end{bmatrix}\right)\\\
={}&\text{rank}\left(\begin{bmatrix}e_{1}&\epsilon
e_{r}&\epsilon^{2}e_{r-1}\cdots&\epsilon^{r-1}e_{2}\end{bmatrix}\right)\\\
={}&r.\end{split}$ (100)

Therefore, the controllability matrix is full row rank, implying that the
system is controllable. Therefore, for all $i$, the matrices
$\\{A_{i},B_{i},C_{i}\\}$ represents a minimal system.

Step 2: Now we show learning at least one of those minimal realizations is
hard. Consider any estimator $\widehat{f}$ mapping $\mathcal{D}$ to
$\left(\widehat{A}(\mathcal{D}),\widehat{B}(\mathcal{D}),\widehat{C}(\mathcal{D})\right)\in\mathbb{R}^{r\times
r}\times\mathbb{R}^{r\times m}\times\mathbb{R}^{n\times r}$.

We now define events $\\{F_{i}\\}_{i\in[|\mathcal{P}|]}$ as follows: