The system identification problem is to learn the parameters of a dynamical system, given the measurements of the inputs and outputs. While classical system identification techniques focused primarily on achieving asymptotic consistency [1, 2, 3], recent efforts have sought to characterize the number of samples needed to achieve a desired level of accuracy in the learned model. In particular, non-asymptotic analysis of system identification based on a single trajectory has been studied extensively over the past few years [4, 5, 6, 7, 8, 9].

In practice, performing system identification using a single trajectory could be problematic for many applications. For example, having the system run for a long time could incur risks when the system is unstable. Furthermore, only historical snippets of data about the system may be available, without the ability to easily observe long-run behavior. Additionally, in settings where one has the ability to restart the system or have several copies of the system running in parallel, one may obtain multiple trajectories generated by the system dynamics [10]. The paper [11] studies the sample complexity of identifying a system whose state is fully measured using only the final data points from multiple trajectories. Using a similar setup, the paper [12] explores the benefits of adding an $\ell_{1}$ regularizer. The paper [13] studies the sample complexity of partially-measured system identification by including nuclear norm regularization, again only using the final samples from each trajectory. For partially-measured systems, the paper [14] allows for more efficient use of data. As mentioned in [14], compared to the single trajectory setup, the multiple trajectories setup usually allows for more direct application of concentration inequalities due to the assumption of independence over multiple trajectories.

In addition to the potential lack of single long trajectories, in many settings we may not be able to actually apply inputs to the system in order to perform system identification; this could be due to the costs of applying inputs, or due to the fact that we are simply observing an autonomous system that we cannot control. The uncontrolled system may also be serving as a subsystem connected to the main system that one wants to control, and having a better model of the subsystem could be useful in controlling the main system. For partially-measured systems, the characterization of finite sample error of purely stochastic systems (systems that are entirely driven by unmeasurable noise) is more challenging as indicated in [15]. There, the goal is to estimate the system matrices as well as the steady state Kalman filter gain of the corresponding system. The paper [15] shows that classical stochastic system identification algorithm can achieve a learning rate of $\mathcal{O}(\frac{1}{\sqrt{\bar{N}}})$ (up to logarithmic factors) for both strictly stable and marginally stable systems, where $\bar{N}$ denotes the number of samples in a single trajectory.

In this paper, we are motivated by the challenge of system identification for partially-measured and autonomous stochastic linear systems (with no controlled inputs) as in [15], but for the case where a single long-run trajectory is unavailable. Existing results on consistency and learning rate of stochastic system identification algorithms (including [15]) typically convert the original system to a statistically equivalent form of the Kalman filter that is assumed to have reached steady state [1, 16, 15]. The analysis is then performed by assuming that the covariance matrix of the initial state of the system is the same as the steady state Kalman filter error covariance matrix, which simplifies the analysis. We note, however, that this assumption is invalid when one has no long run observation of the system trajectory, since it is in general unclear how long one should wait until the Kalman filter “converges” (even if it converges exponentially fast) for an unknown system. Furthermore, the available short trajectories may not be long enough to guarantee that the underlying filter has converged. Consequently, the single trajectory-based results cannot be directly applied to the multiple (short) trajectories case. Our goal in this paper is to estimate the system matrices (up to similarity transformations) using only multiple trajectories of transient responses of a partially-measured system that is entirely driven by noise.

Let $\mathbb{R}$ and $\mathbb{N}$ denote the sets of real numbers and natural numbers, respectively. Let $\sigma_{n}(\cdot)$ and $\sigma_{min}(\cdot)$ be the $n$-th largest and smallest singular value, respectively, of a symmetric matrix. Similarly, let $\lambda_{min}(\cdot)$ be the smallest eigenvalue of a symmetric matrix. We use $*$ to denote the conjugate transpose of a given matrix. The spectral radius of a given matrix is denoted as $\rho(\cdot)$. A square matrix $A$ is called strictly stable if $\rho(A)<1$, marginally stable if $\rho(A)\leq 1$, and unstable if $\rho(A)>1$. We use $\|\cdot\|$ to denote the spectral norm of a given matrix. Vectors are treated as column vectors. A Gaussian distributed random vector is denoted as $u\sim\mathcal{N}(\mu,\Sigma)$, where $\mu$ is the mean and $\Sigma$ is the covariance matrix. We use $I$ to denote the identity matrix. We use $\mathbb{E}$ to denote the expectation. The symbol $\prod_{t=i}^{j}A_{t}$ is used to denote the matrix product, $A_{i}A_{i+1}\cdots A_{j}$. The symbol $\dagger$ is used to denote the pseudoinverse.

Consider a discrete time linear time-invariant system with no user specified inputs:

where $x_{k}\in\mathbb{R}^{n}$, $y_{k}\in\mathbb{R}^{m}$, $w_{k}\in\mathbb{R}^{n}$, $v_{k}\in\mathbb{R}^{m}$, $A\in\mathbb{R}^{n\times n}$ and $C\in\mathbb{R}^{m\times n}$. The noise terms $w_{k}$ and $v_{k}$ are assumed to be i.i.d Gaussian, i.e., $w_{k}\sim\mathcal{N}(0,Q)$, $v_{k}\sim\mathcal{N}(0,R)$. The initial state is also assumed to be independent of $w_{k}$ and $v_{k}$, and is distributed as $x_{0}\sim\mathcal{N}(\mu,\Sigma_{0})$. In addition, whether $\mu$ is zero or non-zero is assumed to be known. If $\mu$ is non-zero, the system matrix $A$ is assumed to be strictly stable or marginally stable. The system order $n$ is also assumed to be known. We refer to the above system as an autonomous stochastic linear system. We will make the following assumption.

Under the above assumption, the Kalman Filter associated with system (1) is a system of the form

where $\hat{x}_{k}$ is an estimate of state $x_{k}$, with the initial estimate being the mean of the initial state in system (1), i.e., $\hat{x}_{0}=\mu$. The sequence of matrices $K_{k}\in\mathbb{R}^{n\times m}$ is called the Kalman gain, given by

where the estimation error covariance $P_{k}\in\mathbb{R}^{n\times n}$ is updated based on the Riccati equation

with $P_{0}=\Sigma_{0}$. Finally, $e_{k}=y_{k}-C\hat{x}_{k}$ are independent zero mean Gaussian innovations with covariance matrix given by

Since the outputs of system (1) and system (2) have identical statistical properties [17], we will perform analysis on system (2). The subspace identification problem for stochastic systems that we tackle in this paper is to identify the system matrices $(A,C)$ up to a similarity transformation, given multiple trajectories of outputs of the system (1). As a byproduct, we will also simultaneously learn the Kalman filter gain $K_{k}$ of the corresponding system, at some time step $k$. In particular, we are interested in the quality of the estimates of $(A,C)$ given a finite number of samples.

Here we describe a variant of classical subspace identification algorithm [17] to estimate $(A,C)$ matrices (up to a similarity transformation). We will first establish some definitions.

Let the batch matrix of initial states be $\hat{X}_{0}\triangleq\begin{bmatrix}\hat{x}_{0}^{1}&\hat{x}_{0}^{2}&\cdots&\hat{x}_{0}^{N}\\ \end{bmatrix}.$ Define the largest norm of innovation covariance matrices as $\mathcal{\bar{R}}_{T}\triangleq\max_{t\in 0,\dots,T-1}\|\bar{R}_{t}\|,$ where $\bar{R}_{t}$ is defined in (4). For any $l\geq 1$, the extended observability matrix $\mathcal{O}_{l}\in\mathbb{R}^{ml\times n}$ and the reversed extended controllability matrix $\mathcal{K}_{p}\in\mathbb{R}^{n\times mp}$ are defined as:

Define

Let $K\in\mathbb{R}^{n\times m}$ be the steady state Kalman gain $K=APC^{*}(CPC^{*}+R)^{-1}$, where $P\in\mathbb{R}^{n\times n}$ is the solution to the Riccati equation, $P=APA^{*}+Q-APC^{*}(CPC^{*}+R)^{-1}CPA^{*}.$ From Kalman filtering theory, the matrix $A-KC$ has spectral radius strictly less than 1 [18]. Denote the reversed extended controllability matrix formed by the steady state Kalman gain $K$ as

We further make the following assumption.

The rank condition on $\mathbf{K}_{p}$ is standard, e.g., [3, 15]. The rank condition on $\mathcal{K}_{p}$ is needed to ensure that $G$ has rank $n$, which can be satisfied in practice by choosing $p$ to be relatively large if the rank condition on $\mathbf{K}_{p}$ is satisfied (see Proposition 8 in the Appendix).

Finally, for any integer $a\geq 0$ and $b\geq 2$, define the block-Toeplitz matrix $\mathcal{T}_{b}^{a}\in\mathbb{R}^{mb\times mb}$ as:

In this section, we will present our main results on bounding the estimation error $\|\hat{G}-G\|$ from (10). We will show that the term $\|(Y_{-}Y_{-}^{*})^{-1}\|$ decreases with a rate of $\mathcal{O}(\frac{1}{N})$, and then upper bound the growth rate of other terms in (10) separately. Using recent results on the balanced realization algorithm with the adjustments to accommodate the non-steady state Kalman filter, we then show that the estimation error of the system matrices $A,C,K_{p-1}$ will also be bounded when the error $\|\hat{G}-G\|$ is small enough.

First, denote the covariance matrix of the weighted past innovation $\mathcal{T}_{p}^{0}E_{-}^{i}$, $1\leq i\leq N$, as:

Let $\sigma_{E}\triangleq\sigma_{min}(\Sigma_{E})$. We first show that the weighted innovation covariance matrix $\Sigma_{E}$ is strictly positive definite. The proof is similar to [15, Lemma 2].

Next we will show that the term $\|(Y_{-}Y_{-}^{*})^{-1}\|$ is decreasing with a rate of $\mathcal{O}(\frac{1}{N})$. Since $Y_{-}=\mathcal{O}_{p}\hat{X}_{0}+\mathcal{T}_{p}^{0}E_{-}$, the explicit form of $Y_{-}Y_{-}^{*}$ is

and we will bound these terms separately.

We will rely on the following lemma from [11, Lemma 2], which provides a non-asymptotic lower bound of a standard Wishart matrix.

To bound the cross terms due to the possibly non-zero batch matrix of initial states in (11), $\hat{X}_{0}E_{-}^{*}$, we will be using the following lemma from [5, Lemma A.1].

Now we are ready to show that $\|(Y_{-}Y_{-}^{*})^{-1}\|$ is decreasing with a rate of $\mathcal{O}({\frac{1}{N}})$.

To see that $N$ will eventually be greater than $N_{1}$ even if $\|\hat{X}_{0}\|$ is non-zero, note that when the system is strictly stable or marginally stable, $\|\mathcal{T}_{p}^{0}\|$ and $\|\mathcal{O}_{p}\|$ will grow no faster than $\mathcal{O}(p^{\mathbf{d}})$ for some constant $\mathbf{d}$ (see Proposition 5 for $\|\mathcal{T}_{p}^{0}\|$, and [15, Corollary E.1] for $\|\mathcal{O}_{p}\|$). Further, $\gamma_{p}$ is $\mathcal{O}(p^{\frac{1}{2}})$, and $\|\hat{X}_{0}\|=\sqrt{N}\|\mu\|=\mathcal{O}(\sqrt{N})$. Thus if $p=\mathcal{O}(\log{N})$, $N$ will eventually be greater than $N_{1}$ as $N$ increases.

Now we will show that the term $\|E_{+}E_{-}^{*}\|$ is $\mathcal{O}(\sqrt{N})$. We will leverage the following lemma from[11, Lemma 1] to bound the product of the innovation terms.

With Proposition 4, we are now in place to prove the bound on the estimation error of $G$ .