Multi-Task System Identification of Similar Linear Time-Invariant Dynamical Systems

In this section, we consider the case where the system matrices $A_{1},\ldots A_{N}$ have zeros in a common set of entries; i.e., $(A_{1})_{ij}=(A_{2})_{ij}=\ldots=(A_{N})_{ij}=0$ for a subset of indices $i\in[n]$ and $j\in[n]$. In the statistical learning literature, this structural similarity leads to a setup where the unknowns $\{A_{1},\ldots,A_{N}\}$ feature a group sparsity [28, 29]. Leveraging the technical approach of [28, 29], we then formulate the multi-task system identification problem for LTI systems with similar connectivity as follows:

where $\lambda\geq 0$ and the function $\mathcal{R}(A_{1},\ldots,A_{N})=\sum_{i=1}^{N}\sum_{j=1}^{N}\|[(A_{1})_{ij},(A_{2})_{ij},\ldots,(A_{N})_{ij}]\|_{2}$ is utilized to enforce group sparsity; that is, the solution of (4) is such that either the whole vector $[(A_{1})_{ij},(A_{2})_{ij},\ldots,(A_{N})_{ij}]$ is zero or not. The parameter $\lambda$ in (3) strikes a balance between the LS fit (in our case, the LS fit for individual LTI systems) and a number of vectors $[(A_{1})_{ij},(A_{2})_{ij},\ldots,(A_{N})_{ij}]$ that are set to zero (as shown shortly).

Problem (4) is an unconstrained convex program; given the composite cost, we consider a proximal-gradient method (with line search) for solving (4) (see, e.g., [38, 39]). The proximal-gradient method is tabulated as Algorithm 1.

We note that the step [S2.3] is in fact a closed-form expression for the proximal map:

and it involves $n^{2}$ parallel computations as in (5). The convergence behavior of Algorithm 1 can be readily analyzed by leveraging the results in [38, Chapter 2]. Moreover, Algorithm 1 can be converted into a “standard” proximal-gradient method if the line search is not performed [40].

From the thresholding operation (5), it is clear that increasing $\lambda$ has the effect of forcing a higher number of entries of the system matrices to be zero. Unfortunately, tuning $\lambda$ is not an easy task, and cross-validation procedures are typically utilized to find the value of $\lambda$ such that the estimated matrices yield the lowest error on test data; see, for example [21, 20].

In this section, we analyze the performance of the multi-task system identification method. The performance of regression problems with sparsity-enforcing regularization terms is oftentimes compared against an “oracle” that represents the best achievable estimation when considering a (group) sparse model, under a given model compatibility condition [30]. We will provide error bounds with respect to both the oracle and the true system matrices.

To simplify the notation, we outline the results for the case where $P_{i}=P$ for all $i\in[N]$ (though, similar results hold when the $\{P_{i}\}$ are different). Let $A^{S}\in\mathbb{R}^{n\times n}$ be a matrix formed by selecting the entries of $A\in\mathbb{R}^{n\times n}$ indexed by $S\subset[n]\times[n]$ and setting zero to the other entries. For $N$ matrices $A_{i}\in\mathbb{R}^{n\times n}$, $i\in[N]$, define for brevity $\|\{A_{i}\}_{i\in[N]}\|_{2,1}=\sum_{i=1}^{n}\sum_{j=1}^{n}\|[(A_{1})_{ij},(A_{2})_{ij},\ldots,(A_{N})_{ij}]\|_{2}$. With this notation in place, we first state the main assumptions and outline the compatibility condition associated with the group sparse model [30, Chapter 8].

We note that, when $S$ coincides with the support of the matrices $\{A_{i}\}_{i\in[n]}$, this technical condition provides bounds on the values of the non-zero entries of the matrices.

Hereafter, we let $\{A_{i}^{MT}\}_{i\in[N]}$ be an optimal solution of the system identification problem (4), and we denote as $\{A_{i}^{\star}\}_{i\in[N]}$ the true matrices of the LTI systems in (1). The following theorem provides an error bound for the proposed multi-task system identification (4) relative to the (true) system matrices $\{A_{i}^{\star}\}_{i\in[N]}$.

From Theorem 1, it can be seen that the average error $E(N,P):=\frac{1}{PN}\sum_{i=1}^{N}\sum_{i=1}^{N}\sum_{\tau=1}^{P}\|(A_{i}^{\star}-A_{i}^{MT})x_{i}(\tau)\|_{2}^{2}$ is of the order

when taking $\lambda=O(\sqrt{P}\sqrt{N+\sqrt{N}+1})$. Interestingly, by increasing the number of LTI systems in the multi-task system identification (i.e. $N$ increases), the average error $E(N,P)$ decreases; however, when $N\rightarrow\infty$, the error does not tend to $0$. This can be understood as a plateau in the ability to “transfer information” between systems.

One possible shortcoming of Theorem 1 is that the set $S^{\star}$ describing the (common) support of the matrices $\{A_{i}^{\star}\}_{i\in[N]}$ is assumed to satisfy the compatibility condition. In the following, we offer an additional error bound for cases where $\{A_{i}^{\star}\}_{i\in[N]}$ is not guaranteed to satisfy the compatibility condition; the error bound leverages the notion of oracle [30].

Theorem 2 asserts that the error incurred by the multi-task system identification (4) is bounded by the estimation error associated with the oracle plus an additional term modeling the error between the estimated matrices and the oracle. Here, the oracle represents the best achievable estimation when considering matrices with support index set that satisfies the compatibility condition.

Consider the case where, for any pair $A_{i},A_{j}$, $i,j\in[N]$, there exists $\epsilon>0$ such $\|A_{i}-A_{j}\|_{F}^{2}\leq\epsilon$. This case is referred to as “small heterogeneity” in [23]. With this prior, the multi-task system identification problem can be formulated as

where we recall that $\mathcal{L}_{i}(A_{i})$ is the LS fit for the $i$th system and where the regularization term penalizes large deviations between the estimated matrices [20].

Problem (12) is convex and can be solved in closed form. However, to avoid computationally-heavy matrix inversions (especially when the dimension $b$ is large and several systems are considered in the system identification process), Algorithm 1 can be utilized to solve (12) by replacing [S2.2] with the following $n^{2}$ parallel computations:

where $s_{ij}=\sum_{\ell=1}^{N}(Z_{\ell})_{ij}$.

From the update (13), it can be seen that $y_{ij}\rightarrow(s_{ij}/N)[1,1,\cdots,1]$ as $\lambda\rightarrow\infty$, thus setting all the system matrices to be the same. On the other hand, by setting $\lambda=0$ one recovers the LS method. Even in this case, cross-validation procedures can be utilized to find the value of $\lambda$ such that the estimated matrices yield the lowest error on test data [21].

The estimation performance of (12) can be analyzed by deriving bounds between an optimal solution of  (12) and the one of the LS method. Since the bound is straightforward, and because of space limitations, we omit this result from the paper.

Lastly, we comment on an additional “similarity” where the system matrices are (approximately) linearly dependent. Precisely, we consider a scenario where for a subset of systems, there exists coefficients $\{\alpha_{i,j}\in\mathbb{R}\}$ with $\alpha_{i,j}\neq 0$ for some $j\in[N]$ such that $A_{i}\approx\sum_{j=1,j\neq i}^{N}\alpha_{ij}A_{j}$.

To formalize the setup, suppose that $q\ll N$ of the matrices $\{A_{i}\}_{i\in[N]}$ are such that the remaining $N-q$ can be represented as a linear combination of these $q$ matrices. Consider then building the $n^{2}\times N$ matrix $[\operatorname{vec}(A_{1}),\operatorname{vec}(A_{2}),\ldots,\operatorname{vec}(A_{N})]$; it follows that this matrix has rank $q\ll N$. Based on this observation, we propose to formulate a multi-task system identification problem for this case as

where the regularization function promotes sparsity in the singular values of the matrix $[\operatorname{vec}(A_{1}),\operatorname{vec}(A_{2}),\ldots,\operatorname{vec}(A_{N})]$ (see, e.g., [31, 32]).

Problem (IV-B) is convex and with a composite cost where the regularization function is not differentiable [31, 32]. Still, the proximal-gradient algorithm tabulated as Algorithm 1 can be modified to solve (IV-B). In particular, one can replace the step [S2.2] is Algorithm 1 with $\textrm{prox}_{\alpha\lambda\left\|[\operatorname{vec}(A_{1}),\operatorname{vec}(A_{2}),\ldots,\operatorname{vec}(A_{N})]\right\|_{*}}(\{Z_{i}\}_{i\in[N]})$; this proximal map affords a closed-form solution given by:

where the singular value decomposition of $[\operatorname{vec}(A_{1}),\operatorname{vec}(A_{2}),\ldots,\operatorname{vec}(A_{N})]$ is $U\textrm{diag}(\{\sigma_{i}\})V$. The matrices $\{Y_{\ell}\}_{\ell\in[N]}$ are then extracted from the columns of $\bar{Y}$. From the computation of $\bar{Y}$, it can be seen that higher values of $\lambda$ lead to a higher number of singular values of $\bar{Y}$ that are set to zero.

While the effectiveness of nuclear norm minimization approaches has been verified numerically in the literature, identifiability results and analytical bounds for the estimation error are available only under given assumptions on the regressors [31, 32] that may not be applicable to (IV-B); see also [33]. Deriving analytical error bounds for (IV-B) is subject of our ongoing investigations.

We test the proposed method for the problem of estimating the dynamics of brain networks, using data corresponding to the resting state fMRI from the Human Connectome Project (HCP)³³3Data available at https://wiki.humanconnectome.org/ [35, 34, 41]. Here, $x_{i}(t)$ is a $116$-dimensional blood-oxygen-level-dependent (BOLD) time series for $116$ parcellations of the brain of the $i$-th subject. Our goal here is to estimate $N=5$ dynamical systems of the form $x_{i}(t+1)=A_{i}x_{i}(t)+w_{i}(t)$, that model the evolution of BOLD signal when the individual is in a resting state, with $w_{i}(t)$ capturing process noise (the model does not contain external inputs $u_{i}$ due to the resting state condition).

Since the matrices $\{A_{i}\}_{i\in[5]}$ are unknown, we consider the following error for each system:

where $n$ is the length of the testing vector, $p$ is the number of testing data and $\bar{x}(k):=\frac{1}{p}\sum_{i=1}^{p}x_{i}(k)$. Note that $1-\mathcal{E}(A)$ is precisely the average $R^{2}$ indicator of [35].

We consider three different methods: (i) the LS estimator (2), which is utilized per individual; (ii) the Least Absolute Shrinkage and Selection Operator (LASSO), which is again utilized per individual as proposed in [35]; and, (iii) the proposed method (3) with the group-sparsity regularization function (4). The rationale behind the group-sparsity is that the brain dynamics should exhibit the same effective connectivity between parcellations, though the remaining entries acknowledge the diversity in intensities of the interactions across individuals. We note that the effectiveness of the LS and LASSO has been experimentally validated in [35], where their estimation accuracy has been compared with several identification methods. Moreover, we performed a cross-validation procedure to optimize the performance of the LASSO.

In Figure 2, we compare the LS, LASSO and our approach (which is labeled as “MT”) in two cases: (a) the same amount of training data is utilized for the five subjects; and, (b) for subject 1, we utilize only 25% of the training data points with respect to the other subjects 2-5. We use $100$ test points. In Figure 2(a) we plot the mean error across the subjects 1-5; in Figure 2(b) we plot the mean error across the subjects 2-5 and the error for subject 1, for which fewer fMRI readings are available. The proposed method outperforms the LS and the LASSO, on par with the number of fMRI scans in both cases. The merits of the proposed method are particularly evident in Figure 2(b), where the proposed method significantly outperforms the LASSO for the subject 1; on the other hand, the LS is ill-conditioned and does not return meaningful estimates. This shows the ability to leverage information and data (in this case, fMRI readings) from the dynamics of subjects 2-5 to assist the estimation of the dynamics in subject 1.

To provide additional comparisons other than the mean error, Figure 3 shows the box plots for the LS, the LASSO, and the proposed approach in three different scenarios. In particular, Figure 3(a) shows that proposed multi-task identification method can achieve a smaller or comparable error (on average) than LS or LASSO when trajectories of 900 time steps are used for each subject (and these training trajectories are sufficiently rich). Figure 3(b) considers the case where $75$ training data points are available for subject 1 and $300$ for subjects 2-5. Here, the LS does not perform well due to ill-conditioning. The performance of the LASSO is comparable with the one of the proposed method in terms of median; however, the proposed method shows smaller upper and lower quartiles. Moreover, Figure 3(c) considers the case where fewer fMRI readings are available for subjects 1 and 3; the proposed method performs better than the LASSO in terms of quartiles and has a significantly less error deviation across the parcellations.

Finally, a representative example of the estimated matrices $\hat{A}_{i}$ for the subjects 1-3 is provided in Figure 4. The estimated matrices are the ones obtained in the case considered in Figure 3(c), where subjects 1 and 3 have fewer training points. It is possible to notice that the three matrices have zeros in many common entries. Based on this result, we will explore additional regularization methods that will combine group sparsity with (entry-wise) sparsity.

We provide additional results on synthetic data. We consider 10 systems as in (1), where $\{A_{i}\}_{i\in[10]}\in\mathbb{R}^{50\times 50}$, $\{B_{i}\}_{i\in[10]}\in\mathbb{R}^{50\times 4}$, $u_{i}(t)$ is the vector of all ones in $\mathbb{R}^{4}$, i.e. $u_{i}(t)$ is constant vector and $w_{i}(t)\sim\mathcal{N}(0,0.1^{2})$. We consider two different cases: common sparsity and linear combinations. We compare the LS estimator (2) and the proposed method (3) with the group-sparsity regularization (4) and nuclear norm regularization (IV-B).

Figure 5 compares the LS and our approach in two cases: (a) all the 10 systems can be represented by a linear combination of 3 systems and only 25% of the training data points are accessible for the tenth system with respect to the other systems 1-9; (b) all the 10 systems have the same sparsity pattern and only 25% of the training data points are accessible for the tenth system with respect to the other systems 1-9. The testing is on $60$ data points. In Figure 5(a), we plot the mean error across systems 1-9 as well as the error for system 10. The proposed method outperforms the LS approach in both the mean error and the error for system 10, especially in the case of only a small number of data available. In Figure 5(b), we can observe similar results.