Joint Learning of Linear Time-Invariant Dynamical Systems

The problem of identifying the transition matrices in linear time-invariant (LTI) systems has been extensively studied in the literature [8, 26, 29]. Recent papers establish finite-time rates for accurately learning the dynamics in various online and offline settings [16, 36, 39]. Notably, existing results are established when the goal is to identify the transition matrix of a single system.

However, in many application areas of LTI systems, one observes state trajectories of multiple dynamical systems. So, in order to be able to efficiently use the full data of all state trajectories and utilize the possible commonalities the systems share, we need to estimate the transition matrices of all systems jointly. The range of applications is remarkably extensive, including dynamics of economic indicators in US states [35, 40, 42], flight dynamics of airplanes at different altitudes [6], drivers of gene expressions across related species [5, 19], time series data of multiple subjects that suffer from the same disease [38, 41], and commonalities among multiple subsystems in control engineering [43].

In all these settings, there are strong similarities in the dynamics of the systems, which are unknown and need to be learned from the data. Hence, it becomes of interest to develop a joint learning strategy for the system parameters, by pooling the data of the underlying systems together and learn the unknown similarities in their dynamics. In particular, this strategy is of extra importance in settings wherein the available data is limited, for example when the state trajectories are short or the dimensions are not small.

In general, joint learning (also referred to as multitask learning) approaches aim to study estimation methods subject to unknown similarities across the data generation mechanisms. Joint learning methods are studied in supervised learning and online settings [10, 4, 32, 33, 3]. Their theoretical analyses are obtained rely on a number of technical assumptions regarding the data, including independence, identical distributions, boundedness, richness, and isotropy.

However, for the problem of joint learning of dynamical systems, additional technical challenges are present. First, the observations are temporally dependent. Second, the number of unknown parameters is the square of the dimension of the system, which impacts the learning accuracy. Third, since in many applications the dynamics matrices of the underlying LTI systems might possess eigenvalues of (almost) unit magnitude, conventional approaches for dependent data (e.g., mixing) inapplicable [16, 36, 39]. Fourth, the spectral properties of the transition matrices play a critical role on the magnitude of the estimation errors. Technically, the state vectors of the systems can scale exponentially with the multiplicities of the eigenvalues of the transition matrices (which can be as large as the dimension). Accordingly, novel techniques are required for considering all important factors and new analytical tools are needed for establishing useful rates for estimation error. Further details and technical discussions are provided in Section 3.

We focus on a commonly used setting for joint learning that involves two layers of uncertainties. It lets all systems share a common basis, while coefficients of the linear combinations are idiosyncratic for each system. Such settings are adopted in multitask regression, linear bandits, and Markov decision processes [14, 22, 31, 44]. From another point of view, this assumption that the system transition matrices are unknown linear combinations of unknown basis matrices can be considered as a first-order approximation for unknown non-linear dynamical systems [27, 30]. Further, these compound layers of uncertainties subsume a recently studied case for mixtures of LTI systems where under additional assumptions such as exponential stability and distinguishable transition matrices, joint learning from unlabeled state trajectories outperforms individual system identification [11].

The main contributions of this work can be summarized as follows. We provide novel finite-time estimation error bounds for jointly learning multiple systems, and establish that pooling the data of state trajectories can drastically decrease the estimation error. Our analysis also presents effects of different parameters on estimation accuracy, including dimension, spectral radius, eigenvalues multiplicity, tail properties of the noise processes, and heterogeneity among the systems. Further, we study learning accuracy in the presence of model misspecifications and show that the developed joint estimator can robustly handle moderate violations of the shared structure in the dynamics matrices.

In order to obtain the results, we employ advanced techniques from random matrix theory and prove sharp concentration results for sums of multiple dependent random matrices. Then, we establish tight and simultaneous high-probability confidence bounds for the sample covariance matrices of the systems under study. The analyses precisely characterize the dependence of the presented bounds on the spectral properties of the transition matrices, condition numbers, and block-sizes in the Jordan decomposition. Further, to address the issue of temporal dependence, we extend self-normalized martingale bounds to multiple matrix-valued martingales, subject to shared structures across the systems. We also present a robustness result by showing that the error due to misspecifications can be effectively controlled.

The remainder of the paper is organized as follows. The problem is formulated in Section 2. In Section 3, we describe the joint-learning procedure, study the per-system estimation error, and provide the roles of various key quantities. Then, investigation of robustness to model misspecification and the impact of violating the shared structure are discussed in Section 4. We provide numerical illustrations for joint learning in Section 5 and present the proofs of our results in the subsequent sections. Finally, the paper is concluded in Section 10.

Notation. For a matrix $A$, $A^{\prime}$ denotes the transpose of $A$. For square matrices, we use the following order of eigenvalues in terms of their magnitudes: $\left|\lambda_{\max}(A)\right|=\left|\lambda_{1}(A)\right|\geq\left|\lambda_{2}(A)\right|\geq\cdots\geq\left|\lambda_{d}(A)\right|=\left|\lambda_{\min}(A)\right|$. For singular values, we employ $\sigma_{\min}(A)$ and $\sigma_{\max}(A)$. For any vector $v\in\mathbb{C}^{d}$, let ${\left|\kern-1.29167pt\left|v\right|\kern-1.29167pt\right|}_{p}$ denote its $\ell_{p}$ norm. We use ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\cdot\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}_{\gamma\rightarrow\beta}$ to denote the matrix operator-norm for $\beta,\gamma\in[1,\infty]$ and $A\in\mathbb{C}^{d_{1}\times d_{2}}$: $\left|\!\left|\!\left|A\right|\!\right|\!\right|_{{\gamma\rightarrow\beta}}=\sup\limits_{v\neq{\bf 0}}{{\left|\kern-1.29167pt\left|Av\right|\kern-1.29167pt\right|}_{\beta}}/{{\left|\kern-1.29167pt\left|v\right|\kern-1.29167pt\right|}_{\gamma}}$. When $\gamma=\beta$, we simply write ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|A\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}_{\beta}$. For functions $f,g:\mathcal{X}\rightarrow\mathbb{R}$, we write $f\lesssim g$, if $f(x)\leq cg(x)$ for a universal constant $c>0$. Similarly, we use $f=O(g)$ and $f=\Omega(h)$, if $0\leq f(n)\leq c_{1}g(n)$ for all $n\geq n_{1}$, and $0\leq c_{2}h(n)\leq f(n)$ for all $n\geq n_{2}$, respectively, where $c_{1},c_{2},n_{1},n_{2}$ are large enough constants. For any two matrices of the same dimensions, we define the inner product $\left\langle A,B\right\rangle=\mathop{\mathrm{tr}}\left(A^{\prime}B\right)$. Then, the Frobenius norm becomes ${\left|\kern-1.29167pt\left|A\right|\kern-1.29167pt\right|}_{F}=\sqrt{\left\langle A,A\right\rangle}$. The sigma-field generated by $X_{1},X_{2},\ldots,X_{n}$ is denoted by $\sigma(X_{1},X_{2},\ldots X_{n})$. We denote the $i$-th component of the vector $x\in\mathbb{R}^{d}$ by $x[i]$. Finally, for $n\in\mathbb{N}$, the shorthand $[n]$ is the set $\left\{1,2,\ldots,n\right\}$.

Our main goal is to study the rates of jointly learning dynamics of multiple LTI systems. Data consists of state trajectories of length $T$ from $M$ different systems. Specifically, for $m\in[M]$ and $t=0,1,\ldots,T$, let $x_{m}(t)\in\mathbb{R}^{d}$ denote the state of the $m$-th system, that evolves according to the Vector Auto-Regressive (VAR) process

Above, $A_{m}\in\mathbb{R}^{d\times d}$ denotes the true unknown transition matrix of the $m$-th system and $\eta_{m}(t+1)$ is a mean zero noise. For succinctness, we use $\Theta^{*}$ to denote the set of all $M$ transition matrices $\{A_{m}\}_{m=1}^{M}$. The transition matrices are related as will be specified in Assumption 3.

Note that the above setting includes systems with longer memories. Indeed, if the states $\tilde{x}_{m}(t)\in\mathbb{R}^{\tilde{d}}$ obey

then, by concatenating $\tilde{x}_{m}(t-1),\cdots,\tilde{x}_{m}(t-q)$ in one larger vector $x_{m}(t-1)$, the new state dynamics is (1), for $d=q\tilde{d}$ and $A_{m}=\begin{bmatrix}B_{m,1}\cdots B_{m,q-1}&B_{m,q}\\ I_{(q-1)\tilde{d}}&0\end{bmatrix}.$

We assume that the system states do not explode in the sense that the spectral radius of the transition matrix $A_{m}$ can be slightly larger than one. This is required for the systems to be able to operate for a reasonable time length [25, 15]. Note that this assumption still lets the state vectors grow with time, as shown in Figure 1.

In addition to the magnitudes of the eigenvalues, further properties of the transition matrices heavily determine the temporal evolution of the systems. A very important one is the size of the largest block in the Jordan decomposition of $A_{m}$, which will be rigorously defined shortly. This quantity is denoted by $l$ in (4). The impact of $l$ on the state trajectories is illustrated in Figure 1, wherein we plot the logarithm of the magnitude of state vectors for linear systems of dimension $d=32$. The upper plot depicts state magnitude for stable systems and for blocks of the size $l=2,4,8,16$ in the Jordan decomposition of the transition matrices. It illustrates that the state vector scales exponentially with $l$. Note that $l$ can be as large as the system dimension $d$.

Moreover, the case of transition matrices with eigenvalues close to (or exactly on) the unit circle is provided in the lower panel in Figure 1. It illustrates that the state vectors grow polynomially with time, whereas the scaling with the block-size $l$ is exponential. Therefore, in design and analysis of joint learning methods, one needs to carefully consider the effects of $l$ and $\left|\lambda_{1}(A_{m})\right|$.

Next, we express the probabilistic properties of the stochastic processes driving the dynamical systems. Let $\mathcal{F}_{t}=\sigma(x_{1:M}(0),\eta_{1:M}(1),\dots,\eta_{1:M}(t))$ denote the filtration generated by the the initial state and the sequence of noise vectors. Based on this, we adopt the following ubiquitous setting that lets the noise process $\{\eta_{m}(t)\}_{t=1}^{\infty}$ be a sub-Gaussian martingale difference sequence. Note that by definition, $\eta_{m}(t)$ is $\mathcal{F}_{t}$-measurable.

The above assumption is widely-used in the finite-sample analysis of statistical learning methods [1, 17]. It includes normally distributed martingale difference sequences, for which Assumption 2 is satisfied with $\sigma^{2}=\lambda_{\max}(C)$. Moreover, if the coordinates of $\eta_{m}(t)$ are (conditionally) independent and have sub-Gaussian distributions with constant $\sigma_{i}$, it suffices to let $\sigma^{2}=\sum_{i=1}^{d}\sigma_{i}^{2}$. We let a common noise covariance matrix for the ease of expression. However, the results simply generalize to covariance matrices that vary with time and across the systems, by appropriately replacing upper- and lower-bounds of the matrices [16, 36, 39].

For a single system $m\in[M]$, its underlying transition matrices $A_{m}$ can be individually learned from its own state trajectory data by using the least squares estimator [16, 36]. We are interested in jointly learning the transition matrices of all $M$ systems under the assumption that they share the following common structure.

This assumption is commonly-used in the literature of jointly learning multiple parameters [14, 44]. Intuitively, it states that each system evolves by combining the effects of $k$ systems. These $k$ unknown systems behind the scene are shared by all systems $m\in[M]$, the weight of each of which is reflected by the idiosyncratic coefficients that are collected in $\beta^{*}_{m}$ for system $m$. Thereby, the model allows for a rich heterogeneity across systems.

The main goal is to estimate $\Theta^{*}=\{A_{m}\}_{m=1}^{M}$ by observing $x_{m}(t)$ for $1\leq m\leq M$ and $0\leq t\leq T$. To that end, we need a reliable joint estimator that can leverage the unknown shared structure to learn from the state trajectories more accurately than individual estimations of the dynamics. Importantly, to theoretically analyze effects of all quantities on the estimation error, we encounter some challenges for joint learning of multiple systems that do not appear in single-system identification.

Technically, the least-squares estimate of the transition matrix of a single system admits a closed form that lets the main challenge of the analysis be concentration of the sample covariance matrix of the state vectors. However, since closed forms are not achievable for joint-estimators, learning accuracy cannot be directly analyzed. To address this, we first bound the prediction error and then use that for bounding the estimation error. To establish the former, after appropriately decomposing the joint prediction error, we study its scaling with the trajectory-length and dimension, as well as the trade-offs between the number of systems, number of basis matrices, and magnitudes of the state vectors. Then, we deconvolve the prediction error to the estimation error and the sample covariance matrices, and show useful bounds that can tightly relate the largest and smallest eigenvalues of the sample covariance matrices across all systems. Notably, this step that is not required in single-system identification is based on novel probabilistic analysis for dependent random matrices.

In the sequel, we introduce a joint estimator for utilizing the structure in Assumption 3 and analyze its accuracy. Then, in Section 4 we consider violations of the structure in (2) and establish robustness guarantees.

In this section, we propose an estimator for jointly learning the $M$ transition matrices. Then, we establish that the estimation error decays at a significantly faster rate than competing procedures that learn each transition matrix $A_{m}$ separately by using only the data trajectory of system $m$.

Based on the parameterization in (2), we solve for $\widehat{\mathbf{W}}=\{\widehat{W}_{i}\}_{i=1}^{k}$ and $\widehat{B}=\left[\widehat{\beta}_{1}|\widehat{\beta}_{2}|\cdots\widehat{\beta}_{M}\right]\in\mathbb{R}^{k\times M}$, as follows:

where $\mathcal{L}(\Theta^{*},\mathbf{W},B)$ is the averaged squared loss across all $M$ systems:

In the analysis, we assume that one can approximately find the minimizer in (3). Although the loss function in (3) is non-convex, thanks to its structure, computationally fast methods for accurately finding the minimizer are applicable. Specifically, the loss function in (3) is quadratic and the non-convexity is the bilinear dependence on $(\mathbf{W},B)$. The optimization in (3) is of the form of explicit rank-constrained representations  [9]. For such problems, it has been shown under mild conditions that gradient descent converges to a low-rank minimizer at a linear rate [48]. Moreover, it is known that methods such as stochastic gradient descent have global convergence, and these bilinear non-convexities do not lead to any spurious local minima [20]. In addition, since the loss function is biconvex in $\mathbf{W}$ and $B$, alternating minimization techniques converge to global optima, under standard assumptions [23]. Nonetheless, note that a near-optimal minimum for the objective function is sufficient, and we only need to estimate the product $WB$ accurately instead of recovering both $W$ and $B$. More specifically, the error of the joint estimator in (3) degrades gracefully in the presence of moderate optimization errors. For instance, suppose that the optimization problem is solved up to an error of $\epsilon$ from a global optimum. It can be shown that an additional term of magnitude $O\left(\epsilon/\lambda_{\min}(C)\right)$ arises in the estimation error, due to this optimization error. Numerical experiments in Section 5 illustrate the implementation of (3).

In the sequel, we provide key results for the joint estimator in (3) and establish the high probability decay rates of $\sum_{m=1}^{M}{\left|\kern-1.29167pt\left|A_{m}-\widehat{A}_{m}\right|\kern-1.29167pt\right|}_{F}^{2}$.

The analysis leverages high probability bounds on the sample covariance matrices of all systems, denoted by

For that purpose, we utilize the Jordan forms of matrices, as follows. For matrix $A_{m}$, its Jordan decomposition is $A_{m}=P^{-1}_{m}\Lambda_{m}P_{m}$, where $\Lambda_{m}$ is a block diagonal matrix; $\Lambda_{m}={\rm diag}(\Lambda_{m,1},\ldots\Lambda_{m,q_{m}})$, and for $i=1,\ldots q_{m}$, each block $\Lambda_{m,i}\in\mathbb{C}^{l_{m,i}\times l_{m,i}}$ is a Jordan matrix of the eigenvalue $\lambda_{m,i}$. A Jordan matrix of size $l$ for $\lambda\in\mathbb{C}$ is

Henceforth, we denote the size of each Jordan block by $l_{m,i}$, for $i=1,\cdots,q_{m}$, and the size of the largest Jordan block for system $m$ by $l^{*}_{m}$. Note that for diagonalizable matrices $A_{m}$, since $\Lambda_{m}$ is diagonal, we have $l^{*}_{m}=1$. Now, using this notation, we define

where $\lambda_{m,1}=\lambda_{1}\left(A_{m}\right)$ and

The quantities in the definition of $\alpha\left(A_{m}\right)$ can be interpreted as follows. The term $\left|\!\left|\!\left|P_{m}^{-1}\right|\!\right|\!\right|_{{\infty\rightarrow 2}}\left|\!\left|\!\left|P_{m}\right|\!\right|\!\right|_{{\infty}}$ is similar to the condition number of the similarity matrix $P_{m}$ in the Jordan decomposition that is used to block-diagonalize the matrix. Moreover, $f\left(\Lambda_{m}\right)$ for stable matrices, and $e^{\rho+1}$ for transition matrices with (almost) unit eigenvalues, capture the long term influences of the eigenvalues. In other words, $f\left(\Lambda_{m}\right)$ indicates the amount that $\eta_{m}(t)$ contributes to the growth of ${\left|\kern-1.29167pt\left|x_{m}(s)\right|\kern-1.29167pt\right|}$, for $s\gg t$ and $\left|\lambda_{m,1}\right|<1-{\rho}/{T}$. When $|\lambda|\approx 1$, ${\left|\kern-1.29167pt\left|x_{m}(s)\right|\kern-1.29167pt\right|}$ scales polynomially with the trajectory length $T$, since influences of the noise vectors $\eta_{m}(t)$ do not decay as $s-t$ grows, because of the accumulations caused by the unit eigenvalues. The exact expressions are in Theorem 1 below. Note that while $f(\Lambda_{m})$ is used to obtain an analytical upper bound for the whole range $\left|\lambda_{m,1}\right|<1-{\rho}/{T}$, it is not tight for small values of $\lambda_{m,1}$ and tighter expressions can be obtained using the analysis in the proof of Theorem 1.

To introduce the following result, we define $\bar{b}_{m}$ next. First, for some $\delta_{C}>0$ that will be determined later, for system $m$, define $\bar{b}_{m}=b_{T}(\delta_{C}/3)+{\left|\kern-1.29167pt\left|x_{m}(0)\right|\kern-1.29167pt\right|}_{\infty}$, where $b_{T}(\delta)=\sqrt{2\sigma^{2}\log\left(2dMT\delta^{-1}\right)}$. Then, we establish high probability bounds on the sample covariance matrices $\Sigma_{m}$ with the detailed proof provided in Section 6.

The above two expressions for $\bar{\lambda}_{m}$ show that for $\left|\lambda_{m,1}\right|<1-{\rho}/{T}$, the largest eigenvalue of the covariance matrix grows linearly in $T$, whereas for $\left|\left|\lambda_{m,1}\right|-1\right|\leq{\rho}/{T}$, the bounds scale exponentially with the multiplicities of the eigenvalues. Note that the bounds in Theorem 1 and the estimation error results stated hereafter require the trajectories for each system to be longer than $T_{0}$. The precise definition for $T_{0}$ can be found in the statement of Lemma 2 in Section 6.

For establishing the above, we extend existing tools for learning linear systems [1, 16, 36, 45]. Specifically, we leverage truncation-based arguments and introduce the quantity $\alpha(A_{m})$ that captures the effect of the spectral properties of the transition matrices on the magnitudes of the state trajectories. Further, we develop strategies for finding high probability bounds for largest and smallest singular values of random matrices and for studying self-normalized matrix-valued martingales.

Importantly, Theorem 1 provides a tight characterization of the sample covariance matrix for each system, in terms of the magnitudes of eigenvalues of $A_{m}$, as well as the largest block-size in the Jordan decomposition of $A_{m}$. The upper bounds show that $\bar{\lambda}_{m}$ grows exponentially with the dimension $d$, whenever $l^{*}_{m}=\Omega(d)$. Further, if $A_{m}$ has eigenvalues with magnitudes close to $1$, then scaling with time $T$ can be as large as $T^{2d+1}$. The bounds in Theorem 1 are more general than $\mathop{\mathrm{tr}}\left(\sum_{t=0}^{T}A_{m}^{t}A^{\prime}_{m}{}^{t}\right)$ that appears in some analyses [36, 39], and can be used to calculate the latter term. Finally, Theorem 1 indicates that the classical framework of persistent excitation [7, 21, 24] is not applicable, since the lower and upper bounds of eigenvalues grow at drastically different rates.

Next, we express the joint estimation error rates.

The proof is provided in Section 7. By putting Theorems 1 and 2 together, the estimation error per-system¹¹1In order to obtain a guarantee for the maximum error over all systems, additional assumptions on the matrix $\left[\beta^{*}_{1}\ldots\beta^{*}_{M}\right]$ are required. This problem falls beyond the scope of this paper and we leave it to a future work. is

The above expression demonstrates the effects of learning the systems in a joint manner. The first term in (7) can be interpreted as the error in estimating the idiosyncratic components $\beta_{m}$ for each system. The convergence rate is $O\left({k}/{T}\right)$, as each $\beta_{m}$ is a $k$-dimensional parameter and for each system, we have a trajectory of length $T$. More importantly, the second term in (7) indicates that the joint estimator in (3) effectively increases the sample size for the shared components $\{W_{i}\}_{i=1}^{k}$, by pooling the data of all systems. So, the error decays as $O({d^{2}k}/{MT})$, showing that the effective sample size for $\{W_{i}\}_{i=1}^{k}$ is $MT$.

In contrast, for individual learning of LTI systems, the rate is known [16, 18, 36, 39] to be

Thus, the estimation error rate in (7) recovers the rate for a single system ($k=1$), and it significantly improves for joint learning, especially when

Note that the above conditions are as expected. First, when $k\approx d^{2}$, the structure in Assumption 3 does not provide any commonality among the systems. That is, for $k=d^{2}$, the LTI systems can be totally arbitrary and Assumption 3 is automatically satisfied. This prevents reductions in the effective dimension of the unknown transition matrices, and also prevents joint learning from being any different than individual learning. Similarly, $k\approx M$ precludes all commonalities and indicates that $\{A_{m}\}_{m=1}^{M}$ are too heterogeneous to allow for any improved learning via joint estimation.

Importantly, when the largest block-size $l^{*}_{m}$ varies significantly across the $M$ systems, a higher degree of shared structure is needed to improve the joint estimation error for all systems. Since $\bm{\kappa}$ and $\bm{\kappa}_{\infty}$ depend exponentially on $l^{*}_{m}$ (as shown in Figure 1 and Theorem 1) and $l^{*}_{m}$ can be as large as $d$, we can have $\log\bm{\kappa}_{\infty}=\log\bm{\kappa}=\Omega(d)$. Hence, in this situation we incur an additional dimension dependence in the error of the joint estimator. Note that such effects of $l^{*}_{m}$ are unavoidable (regardless of the employed estimator). Moreover, in this case, joint learning rates improve if $k\leq d\text{ and }kd\leq M$. Therefore, our analysis highlights the important effects of the large blocks in the Jordan form of the transition matrices.

The above is an inherent difference between estimating dynamics of LTI systems and learning from independent observations. In fact, the analysis established in this work includes stochastic matrix regressions that the data of system $m$ consists of

wherein the regressors $x_{m}(t)$ are drawn from some distribution $\mathcal{D}_{m}$, and $y_{m}(t)$ is the response. Assume that $\left(x_{m}(t),y_{m}(t)\right)$ are independent as $m,t$ vary. Now, the sample covariance matrix $\Sigma_{m}$ for each system does not depend on $A_{m}$. Hence, the error for the joint estimator is not affected by the block-sizes in the Jordan decomposition of $A_{m}$. Therefore, in this setting, joint learning always leads to improved per-system error rates, as long as the necessary conditions $k<d^{2}$ and $k<M$ hold.

In Theorem 2, we showed that Assumption 3 can be utilized for obtaining an improved estimation error, by jointly learning the $M$ systems. Next, we consider the impacts of misspecified models on the estimation error and study robustness of the proposed joint estimator against violations of the structure in Assumption 3.

Let us first consider the deviation of the dynamics of each system $m\in[M]$ from the shared structure. Specifically, by employing the matrix $D_{m}$ to denote the deviation of system $m$ from Assumption 3, suppose that

Then, denote the total misspecification by $\bar{\zeta}^{2}=\sum_{m=1}^{M}{\left|\kern-1.29167pt\left|D_{m}\right|\kern-1.29167pt\right|}_{F}^{2}$. We study the consequences of the above deviations, assuming that the same joint learning method as before is used for estimating the transition matrices.

The proof of Theorem 3 is provided in Section 8. In (11), we observe that the total misspecification $\bar{\zeta}^{2}$ imposes an additional error of $(\bm{\kappa}_{\infty}+1)\bar{\zeta}^{2}$ for jointly learning all $M$ system. Hence, to obtain accurate estimates, we need the total misspecification $\bar{\zeta}^{2}$ to be smaller than the number of systems $M$, as one can expect. The discussion following Theorem 2 is still applicable in the misspecified setting and indicates that in order to have accurate estimates, the number of the shared bases $k$ must be smaller than $M$ as well. In addition, compared to individual learning, the joint estimation error improves despite the unknown model misspecifications, as long as

This shows that when the total misspecification is proportional to the number of systems; $\bar{\zeta}^{*}=\Omega(M)$, we pay a constant factor proportional to $\bm{\kappa}_{\infty}$ on the per-system estimation error. Note that in case all systems are stable, according to Theorem 1, the maximum condition number $\bm{\kappa}_{\infty}$ does not grow with $T$, but it scales exponentially with $l^{*}_{m}$. The latter again indicates an important consequence of the largest block-sizes in Jordan decomposition that this work introduces.

Moreover, when a transition matrix $A_{m}$ has eigenvalues close to or on the unit circle in the complex plane, by Theorem 1, the factor $\bm{\kappa}_{\infty}$ grows polynomially with $T$. Thus, for systems with infinite memories or accumulative behaviors, misspecifications can significantly deteriorate the benefits of joint learning. Intuitively, the reason is that effects of notably small misspecifications can accumulate over time and contaminate the whole data of state trajectories, because of the unit eigenvalues of the transition matrices $A_{m}$. Therefore, the above strong sensitivity to deviations from the shared model for systems with unit eigenvalues seems to be unavoidable.

For example, if for the total misspecification we have $\bar{\zeta}^{2}=O(M^{1-a})$, for some $a>0$, joint estimation improves over the individual estimators, as long as $T{\bm{\kappa}_{\infty}}\lesssim M^{a}{d^{2}}$. Hence, when all systems are stable, the joint estimation error rate improves when the number of systems satisfies $T^{1/a}\lesssim M$. Otherwise, idiosyncrasies in system dynamics dominate the commonalities. Note that larger values of $a$ correspond to smaller misspecifications. On the other hand, Theorem 3 implies that in systems with (almost) unit eigenvalues, the impact of $\bar{\zeta}^{2}$ is amplified. Indeed, by Theorem 1, for unit-root systems, joint learning improves over individual estimators when $d^{2}M^{a}\gg T^{2l^{*}_{m}+2}$. That is, for benefiting from the shared structure and utilizing pooled data, the number of systems $M$ needs to be as large as $T^{(2l^{*}_{m}+2)/a}/d^{2/a}$.

In contrast, if $\bar{\zeta}^{2}=O(M^{1-a})$ for some $a>0$, the joint estimation error for the regression problem in (9) incurs only an additive factor of $O(1/M^{a})$, regardless of the largest block-sizes in the Jordan decompositions and unit-root eigenvalues. Thus, Theorem 3 further highlights the stark difference between joint learning from independent, bounded, and stationary observations, and from state trajectories of LTI systems.

We complement our theoretical analyses with a set of numerical experiments which demonstrate the benefit of jointly learning the systems. We investigate two main aspects of our theoretical results: (i) benefits of joint learning when the $M$ systems share a common linear basis, for different values of $M$, and (ii) interplay of the spectral radii of the system matrices with the joint-estimation error. To that end, we compare the estimation error for the joint estimator in (3) against the ordinary least-squares (OLS) estimates of the transition matrices for each system individually. For solving (3), we use a minibatch gradient-descent-based implementation with Adam as the optimization algorithm [28]. Due to the bilinear form of the optimization objective, gradient descent methods can lead to convergence and computational issues for $\widehat{\mathbf{W}}$ and $\widehat{B}$. Although prior studies utilize norm regularizations to address this issue in some cases [44], we do not use any such regularization in our objective function in (3). Notably, our unregularized minimization exposes no convergence issue in the simulations we performed.

For generating the systems, we consider settings with the number of bases $k=10$, dimension $d=25$, trajectory length $T=200$, and the number of systems $M\in\{1,10,20,50,100,200\}$. We simulate two cases:
(i) the spectral radii are in the range $[0.7,0.9]$, and
(ii) all systems have an eigenvalue of magnitude $1$.

The matrices $\{W_{i}\}_{i=1}^{10}$ are generated randomly, such that each entry of $W_{i}$ is sampled independently from the standard normal distribution $N(0,1)$. Using these matrices, we generate $M$ systems by randomly generating the idiosyncratic components $\beta_{m}$ from a standard normal distribution. For generating the state trajectories, noise vectors are isotropic Gaussian with variance $4$. Additional numerical simulations using Bernoulli random matrices are provided in Appendix A.

We simulate the joint learning problem both with and without model misspecifications. For the latter, deviations from the shared structure are simulated by the components $D_{m}$, which are added randomly with probability $1/M^{a}$ for $a\in\{0,0.25,0.5\}$. The matrices $D_{m}$ are generated with independent Gaussian entries of variance $0.01$, leading to ${\left|\kern-1.29167pt\left|D_{m}\right|\kern-1.29167pt\right|}_{F}^{2}\approx 6.25$ and $\bar{\zeta}^{2}\approx 6.25~{}M^{1-a}$, according to the dimension $d=25$.

To report the results, for each value of $M$ in Figure 2 (resp. Figure 3), we average the errors from $10$ (resp. $20$) random replicates and plot the standard deviation as the error bar. Figure 2 depicts the estimation errors for both stable and unit-root transition matrices, versus $M$. It can be seen that the joint estimator exhibits the expected improvement against the individual one.

More interestingly, in Figure 3(a), we observe that for stable systems, the joint estimator performs worse than the individual one, when significant violations from the shared structure occurs in all systems (i.e., $a=0$). Note that it corroborates Theorem 3, since in this case the total misspecification $\bar{\zeta}^{2}$ scales linearly with $M$. However, if the proportion of systems which violate the shared structure in Assumption 3 decreases, the joint estimation error improves as expected ($a=0.25,0.5$).

Figure 3(b) depicts the estimation error for the joint estimator under misspecification for systems that have an eigenvalue on the unit circle in the complex plane. Our theoretical results suggest that the number of systems needs to be significantly larger in this case to circumvent the cost of misspecification in joint learning. The figure corroborates this result, wherein we observe that the joint estimation error is larger than the individual one, if all systems are misspecified (i.e., $a=0$). Decreases in the total misspecification (i.e., $a=0.25,0.5$) improves the error rate for joint learning, but requires larger number of systems than the stable case.

Finally, we discuss the choice of the number of bases $k$ for applying the joint estimator to real data. It can be handled by model selection methods such as elbow criterion and information criteria [2, 37], as well as robust estimation methods in panel data and factor models [12, 13]. In fact, for all $k^{\prime}\geq k$, the structural assumption is satisfied and leads to similar learning rates, while $k^{\prime}<k$ can lead to larger estimation errors. In Figure 4, we provide a simulation (with $T=250,M=50$) and report the per-system estimation error, as well as the prediction error on a validation data (which is a subset of size $50$). Across all $10$ runs in the experiment, we observed that if the hyperparameter $k^{\prime}$ is chosen according to the elbow criteria, the resulting number of basis models is either equal to the true value $k=10$, or slightly larger. For misspecified models, the optimal choice of $k^{\prime}$ can vary, in the sense that large misspecifications can be added to the shared basis (i.e., $k^{\prime}>k$).

In this and the following sections, we provide the detailed proofs for our results. We start by analysing the sample covariance matrix for each system which is then used to derive the estimation error rates in Theorem 2 and Theorem 3. For clarity, some details of the proofs are delegated to Appendix C. In Section 9, we provide the general probabilistic inequalities that are used throughout the proofs. Now, we prove high probability bounds for covariance matrices $\Sigma_{m}=\Sigma_{m}(T)=\sum_{t=0}^{T}x_{m}(t)x_{m}(t)^{\prime}$ in Theorem 1.