


frame \(\Phi =(\phi _1,\dots ,\phi _{|\mathfrak {S}|})\)  and design parameters \(q,\epsilon >0\) 
set of vectors \(\Phi _s\)  and weight function \(w_s\)

\(\Phi _s=\varnothing \) , \(w_s(.) = 0\) , \(\mathbf {S}_s= {\bf 0}\) 
\(k=1\)  to \(q\) 
sample an element from \(\Phi \)  with probability distribution \(\pi \)  \(\rightarrow \phi \) 
update weight function
: \(w_s(\phi ) \leftarrow w_s(\phi ) + (q\pi (\phi ))^{-1}\)

\(\phi \notin \Phi _s\) ,
add \( \phi \)  to \(\Phi _s\) 
update the frame matrix:
\( \mathbf {S}_s \leftarrow \mathbf {S}_s+\phi \phi ^T\)


Theorem 5 
For a given frame \(\Phi \)  in \(\mathbb {R}^n\) , let us fix parameter \(\epsilon \in (1/\sqrt{n},1]\)  and the number of samples \(q=O(n \log n/\epsilon ^2)\) . Then, the resulting set of elements \(\Phi _s\)  from Algorithm REF  is also a frame for \(\mathbb {R}^n\)  with probability at least \(1/2\) . Furthermore, with probability at least \(1/4\) , the estimation-quality losses satisfyUsing monotonicity property of the estimation measures, one can show that upper bounds in (REF ) and () are nonnegative.
\(\dfrac{\rho _d(\Phi _s)-\rho _d(\Phi )}{\rho _d(\Phi )} &\le & -1+\sqrt{\dfrac{4\hspace{1.42271pt}\bar{\chi }}{1-\epsilon }}  \\{\rho _e(\Phi _s)-\rho _e(\Phi )} & \le & n \log \left( \dfrac{4\hspace{1.42271pt}\bar{\chi }}{1-\epsilon } \right), \) 
where \(\bar{\chi }:=\mathbb {E}\left\lbrace  \chi \right\rbrace \)  and \(\chi \)  is a random variable given by
\(\chi :=\inf \Bigg \lbrace {\gamma >0}~\Bigg |~ \sum _{\phi \in \Phi } w_s(\phi ) \big (\gamma - w_s(\phi ) \big )\hspace{1.42271pt}\phi \phi ^T \succeq {\bf 0} \Bigg \rbrace .\) 

The backbone of this result is based on Theorem 1 of [1]} and asserts that Algorithm REF  trims off a given (highly redundant) frame and returns, with probability more than \(0.5\) , a new frame whose size is almost linear in network size. Moreover, it is shown that, with probability at least \(0.25\) , the estimation measures of the new (sparsified) frame stays within constant multiples/difference of the estimation measure of the original (redundant) frame.
Corollary 4 
The random variable \(\chi \) , which is defined by (REF ), satisfies
\(\chi \hspace{1.42271pt}\le \hspace{1.42271pt}\max _{\phi \in \Phi }\hspace{1.42271pt}w_s(\phi )\) 
almost surely.
Moreover, under the settings of Theorem REF , inequalities
\(\dfrac{\rho _d(\Phi _s)-\rho _d(\Phi )}{\rho _d(\Phi )} &\le & -1+ \sqrt{\dfrac{4\hspace{1.42271pt}w_{\max }}{1-\epsilon }}  \\{\rho _e(\Phi _s)-\rho _e(\Phi )} & \le & n \log \left( \dfrac{4\hspace{1.42271pt}w_{\max }}{1-\epsilon } \right) \) 
holds with probability at least \(1/4\) , where
\(w_{\textrm {max}}:=\mathbb {E}\left\lbrace  \hspace{1.42271pt}\displaystyle \max _{\phi \in \Phi }\hspace{1.42271pt}w_s(\phi ) \hspace{1.42271pt}\right\rbrace .\) 

This corollary shows that there exists a clear relationship (even if it is not tight) between the performance loss and the magnitude of parameter \(w_{\max }\) .
The leverage scores disclose the importance of every component with respect to the entire frame for the sake of estimation. For instance, if the network is asymptotically stable, a component with a relatively large time label is expected to have a relatively small leverage score. Thus, such insignificant components are less likely to be sampled by Algorithm REF  and can be trimmed off to achieve a comparable estimation quality.

Running Time Analysis: Computing the inverse of \(\mathbf {S}\)  can be done in \(O(n^3)\)  operations, while computing the leverage scores using this matrix requires \(O(|\Phi |n^2)\) . We need to check for repeated samples, which does not increase the running time of the algorithm. Computing the frame matrix \(\mathbf {S}_s\)  can be done in \(O(n\log n /\epsilon ^2 \times n^2)=O(n^3 \log n/\epsilon ^2)\)  operations. Hence, the total running time of Algorithm REF  is \(O(n^3 \log n /\epsilon ^2+|\Phi |n^2)\) .
 Random Partitioning and Kadison-Singer Paving Solution
If the leverage scores (REF ) are uniformly bounded by a small enough number, then components of a frame can be partitioned in a balanced manner in order to obtain two separate subframes with explicit bounds on their spectra. Such spectral bounds are useful to find bounds on the estimation quality of the resulting subframes. The next theorem is based on Corollary 1.3 of the recent seminal paper [2]} that gives us a paving solution to the famous Kadison-Singer problem.
Proposition 6  Suppose that the leverage scores in (REF ) satisfy
\(r_{\phi }(\mathbf {S}) \le r^* \) 
for all \(\phi \in \Phi \)  and some positive number \(r^* < 1.5 - \sqrt{2}\) . Let us randomly partition \(\Phi \)  into two subfamilies \(\Phi _1\)  and \(\Phi _2\)  such that every element of \(\Phi \) , independent of others, belongs to either of the partitions with probability \(1/2\) .
Then, with a positive probability, the resulting partition will satisfy
\(\begin{adjustbox}{max width=300pt}\left(1-\dfrac{(1+\sqrt{2r^*})^2}{2} \right)\mathbf {S}\preceq \sum \limits _{\phi \in \Phi _j} \phi \phi ^T\preceq \left(1+\dfrac{(1+\sqrt{2r^*})^2}{2} \right)\mathbf {S}\end{adjustbox}\) 

for \(j=1\)  and 2.

For a given (highly redundant) observability frame, the result of Proposition REF  allows us to calculate the relative/absolute estimation-quality degradation of the resulting partitions.
Theorem 6  Suppose that the leverage scores (REF ) satisfy (REF ). Then, the randomly partitioned subfamilies \(\Phi _1\)  and \(\Phi _2\)  from Proposition REF  are both frames for \(\mathbb {R}^n\)  with a positive probability and the estimation-quality can be bounded as follows
\(\dfrac{\rho _d(\Phi _j)-\rho _d(\Phi )}{\rho _d(\Phi )} & \le & { \kappa }(r^*)  \\{\rho _e(\Phi _j)-\rho _e(\Phi )} & \le & -n \log \left(1 - \frac{(\sqrt{2r^*} + 1)^2}{2}\right) \) 
for \(j=1,2\) , where
\({\kappa (r)}:=\left(1 - \frac{(\sqrt{2r} + 1)^2}{2}\right)^{-1/2} - 1.\) 

The quantity \(\kappa (r^*)\)  is a worst-case bound on the relative performance degradation of the randomly partitioned subframes \(\Phi _1\)  and \(\Phi _2\) . This function has been illustrated in Fig. REF . In simulations, we observe that comparably better bounds are achievable.
<FIGURE>
Applicability of Random Partitioning: When we deal with massive incoming samples (data) from the sensors, we expect the leverage scores to be small. In fact, we observe that the leverage scores satisfy
\(\sum _{\phi \in \Phi } r_{\phi }(\mathbf {S})=n ~\Rightarrow ~ \bar{r}(\mathbf {S}) = \dfrac{n}{|\Phi |},\) 
where \(\bar{r}(.)\)  stands for the average of the leverage scores. Hence, as a rule of thumb, one should expect that for fairly balanced observability frames with size
\(|\Phi | > \big \lceil (6+4\sqrt{2})n \big \rceil ,\) 
the leverage score of a typical component in \(\Phi \) , in average, is less than \((6+4\sqrt{2})^{-1}=1.5-\sqrt{2}\approx 0.0858\) .

Running Time Analysis: Computing the frame matrix for each partition can be done in \(O(|\Phi | n^2)\) . Hence, the random partitioning can be done in \(O(|\Phi | n^2)\)  operations.

 Greedy Sparsification
In order to maintain a predetermined level of estimation quality and sparsity, one may consider using greedy algorithms that have been demonstrated to be useful in practice with satisfactory performance for a broad range of combinatorial problems [3]}, [4]}. In greedy frame sparsification, the core idea at every iteration is to eliminate that component of the frame which will increase value of a given estimation measure less than the others.
At iteration \(k\) , let us denote the remaining frame and its frame matrix by \(\Phi _k\)  and \(\mathbf {S}_k\) , respectively. Eliminating a component \(\phi \)  from \(\Phi _k\)  corresponds to the following rank-one update
\(\mathbf {S}_{k+1}=\mathbf {S}_k-\phi \phi ^T\) 

with \(\mathbf {S}_0=\mathbf {S}\) . According to the Sherman-Morrison formula [5]}, one gets update rule
\(\mathbf {S}_{k+1}^{-1}=\mathbf {S}_k^{-1}+\dfrac{\mathbf {S}_k^{-1} \phi \phi ^T \mathbf {S}_k^{-1}}{1-\phi ^T\mathbf {S}_k^{-1}\phi }.\) 
[t]
Greedy Frame Sparsification


initialize:




frame \(\Phi \)  and its frame matrix \(\mathbf {S}\succ 0\) , \(\mathfrak {s} \in (0,1)\) , \(\mathfrak {e}>0\) 
frame \(\Phi _s\) 
\(\mathbf {S}_s=\mathbf {S}\) , \(\Phi _s=\Phi \)

\(\dfrac{\rho _{\Box }(\Phi _s)-\rho _{\Box }(\Phi )}{\rho _{\Box }(\Phi )} \hspace{1.42271pt}\le \hspace{1.42271pt}\mathfrak {e}\)  and \(\dfrac{|\Phi _s|}{|\Phi |} \hspace{1.42271pt}\ge \hspace{1.42271pt}\mathfrak {s}\)

find minimizer \(\phi _{\Box }^*\)  via solving (REF ) or (REF )
update the frame matrix
\(\mathbf {S}_{s}^{-1}\leftarrow \mathbf {S}_s^{-1}+\dfrac{\mathbf {S}_s^{-1} \phi _{\Box }^*\phi _{\Box }^{*T} \mathbf {S}_s^{-1}}{1-\phi _{\Box }^{*T}\mathbf {S}_s^{-1}\phi _{\Box }^*}\) 
update \(\Phi _s \leftarrow \Phi _s \backslash \phi _{\Box }^*\)

Proposition 7  Upon eliminating a component \(\phi \)  from an observability frame \(\Phi _k\) , the estimation measures are updated according to
\(\rho _d(\Phi _{k+1}) & = & \sqrt{\rho _d^2(\Phi _{k}) \hspace{1.42271pt}+ \hspace{1.42271pt}\dfrac{\sigma ^2 \hspace{1.42271pt}r_{\phi }(\mathbf {S}_k^2)}{1-r_{\phi }(\mathbf {S}_k)}}  \\\rho _e(\Phi _{k+1}) & = & \rho _e(\Phi _k)-\log \big (1-r_{\phi }(\mathbf {S}_k) \big ). \) 

At every iteration, the optimizer of \(\rho _d(\Phi _{k+1})\)  can be determined by solving the optimization problem
\(\phi _d^* = \arg \min _{\phi \in \Phi _k} ~ \frac{ \left\Vert  \mathbf {S}_k^{-1} \phi \right\Vert ^2}{1-\phi ^T\mathbf {S}_k^{-1}\phi } \) 
and for \( \rho _e(\Phi _{k+1})\)  by solving
\(\phi _e^* = \arg \min _{\phi \in \Phi _k} ~\phi ^T\mathbf {S}_k^{-1}\phi .\) 
Algorithm REF  details all necessary steps to compute a sparsfication of a (redundant) observability frame, where
we use notation \(\Box \in \lbrace d,e\rbrace \) . The algorithm stops whenever either a desired sparsity level \( \mathfrak {s} \in (0,1)\)  or a maximum allowable relative estimation error \(\mathfrak {e} >0\)  has been achieved. This algorithm resembles the procedure of updating a performance measure of a linear consensus network when a new coupling link is added to the network [6]}. The performance guarantees of the greedy methods in this context is a well-studied subject. In general, derivation of performance bounds heavily depends on the curvature conditions of the specific class of objective functions, e.g., sub-modularity, super-modularity or weak forms of these properties; please see [7]} and references in there. 

Running Time Analysis: Computing \(\mathbf {S}^{-1}\)  at the beginning requires \(O(n^3)\) . Then, at each iteration, one needs to compute and update the value of estimation measure for every vector in \(\Phi \) , which takes \(O(|\Phi | n^2)\) . Thus, in order to achieve sparsity level \(\mathfrak {s}\) , one needs \(O \big ((1-\mathfrak {s})|\Phi |\times |\Phi | n^2 \big )=O \big ((1-\mathfrak {s})|\Phi |^2n^2 \big )\)  operations. Since \(|\Phi | \ge n\) , Algorithm REF  can be implemented in \(O \big ((1-\mathfrak {s})|\Phi |^2 n^2 \big )\) . Compared to running time of the randomized sparsification \(O(n^3 \log n /\epsilon ^2+|\Phi |n^2)\) , the running time of the greedy method can be higher by an order of \(|\Phi |\) .
Remark 4 Our proposed algorithms in this section employ some results from [1]}, [2]}. The idea of sparsification via effective resistances are recently applied in controls community to obtain network abstraction as well as actuator scheduling in large-scale networked control systems [10]} and [11]}. The Kadison-Singer paving solution of [2]} is also utilized in [10]} as a method of randomized actuator scheduling.

Fundamental Limits and Tradeoffs
For a given linear network (REF ) whose state vector is sampled based on an arbitrary sampling strategy, we show that there are fundamental limits and tradeoffs on the best achievable values for the estimation measures and space-time sparsity.
Theorem 7 
Suppose that state of the \(n\) -dimensional linear network (REF ) is sampled under a sampling strategy \(\mathfrak {S}\)  with total number of samples \(|\mathfrak {S}|\) . Then, the best achievable estimation measures are bounded from below by constants that are quantified by
\(\rho _d(\Phi ) \ge {\dfrac{\sigma n}{\nu \hspace{1.42271pt}\sqrt{|\mathfrak {S}|} }}\) 

and
\(\rho _e(\Phi ) \ge n \log \left( \dfrac{n}{\nu ^2 \hspace{1.42271pt}|\mathfrak {S}|} \right),\) 

where \(\Phi \)  is the observability frame corresponding to \(\mathfrak {S}\)  and \(\nu :=\nu (\mathbf {A},\mathfrak {S})\)  is defined by
\(\nu (\mathbf {A},\mathfrak {S})=\max _{ t \in \Theta _1 \cup \dots \cup \Theta _{|\Omega |} } \big \Vert \mathrm {e}^{\mathbf {A}t} \big \Vert . \) 

The inequalities (REF ) and (REF ) give us some convenient rules of thumb about the scaling properties of the estimation measures. For instance, if \(|\mathfrak {S}| =O(n^2)\) , then it can be deduced from (REF ) thatIn the inequality (REF ), the value of \(\sigma /\nu \)  can be interpreted as noise-to-signal ratio because the value of \(\nu \)  relates to the norm of samples and \(\sigma \)  is the standard deviation of measurement noise. If \(\mathbf {A}\)  is Hurwitz, \(\nu \)  is always a finite number. On the other hand, for unstable networks, the noise-to-signal ratio loses its significance as sampling process is prolonged. As a consequence, the value of lower bounds in the inequalities (REF ) and (REF ) are usually small(er) for unstable networks.
\(\rho _d(\Phi ) \ge \frac{c\sigma }{\nu }\) 
for some constant \(c\) . This implies that the estimation quality cannot be enhanced beyond a hard limit. Such limitations are important in network design as they are independent of the network size and sampling strategy. For more discussions on significant role of fundamental limits in control, we refer to [14]}, [15]}.
Theorem 8 
For a given linear network (REF ), let us assume that there exists \(\delta >0\)  such that all distinct eigenvalues of its state matrix satisfy
\(\big (\lambda _i(\mathbf {A})-\lambda _k(\mathbf {A}) \big ) \hspace{1.42271pt}\delta \notin 2 \pi \mathrm {j}\ \mathbb {Z}\) 
and it is sampled according to a sampling strategy with propertyThere is a practical implication for this assumption: sensors usually take samples with some certain frequency based on a synchronized digital clock. \(\Theta _i \subset \delta \mathbb {Z}_+\)  for all \(i \in \Omega \) . If \(\mathbf {A}\)  is Hurwitz, then universal (i.e., independent of number of samples) fundamental limits on the best achievable estimation measures emerge as follows
\(\rho _d(\Phi ) & \ge & \sigma \sqrt{\mathrm {Tr}(\mathbf {Q}^{-1})}, \\\rho _e(\Phi ) & \ge & -\mathrm {Tr}\big (\log (\mathbf {Q}) \big ) ,\) 
where \(\mathbf {Q} \succ 0\)  is the observability Gramian, i.e., the unique solution of Lyapunov equation
\(\mathrm {e}^{\mathbf {A}^T \delta } \hspace{1.42271pt}\mathbf {Q} \hspace{1.42271pt}\mathrm {e}^{\mathbf {A}\delta } \hspace{1.42271pt}- \hspace{1.42271pt}\mathbf {Q} \hspace{1.42271pt}+ \hspace{1.42271pt}\mathbf {C}_{\Omega }^T \hspace{1.42271pt}\mathbf {C}_{\Omega } \hspace{1.42271pt}= \hspace{1.42271pt}0.\) 

Also, the lower bounds can be achieved if and only if
\(\mathfrak {S} \hspace{1.42271pt}= \hspace{1.42271pt}\Omega \times \delta \mathbb {Z}_+.\)

In an exponentially stable linear network, as time goes by, the magnitude of state dwindles compared to measurement noise. For such systems, Theorem REF  predicts that estimation quality cannot be improved beyond some certain threshold even if the number of samples tends to infinity.
For a given sampling strategy \(\mathfrak {S}\) , the average number of samples per subsystem is quantified by
\(\bar{\Theta }(\mathfrak {S}):=\dfrac{1}{|\Omega |}\sum _{i \in \Omega } |\Theta _i(\mathfrak {S})|.\) 
Theorem 9 
For given desired levels of estimation qualities \(\rho _d^*, \rho _e^* >0 \)  and parameter \(\nu ^* >0\) , let us consider
all \(n\) -dimensional linear networks (REF ) whose pair of state matrices and sampling strategies belong to
\( \mathfrak {N}_{\Box }=\Big \lbrace  (\mathbf {A},\mathfrak {S})~\Big |~ \nu (\mathbf {A},\mathfrak {S}) =\nu ^*~~\textrm {and}~~\rho _{\Box }\big (\Phi (\mathbf {A},\mathfrak {S}) \big ) = \rho ^*_{\Box } \Big \rbrace , \) 
where \( \Box \in \lbrace d,e\rbrace \) . Intrinsic tradeoffs between the number of sampling locations and the average number of samples per subsystem transpire over \(\mathfrak {N}_{\Box }\)  that are characterized by
\(\bar{\Theta }(\mathfrak {S}) \cdot |\Omega | ~ \ge ~ \left({\dfrac{\sigma n}{\nu ^* \hspace{1.42271pt}\rho _d^* }}\right)^2,\) 
\(\bar{\Theta }(\mathfrak {S}) \cdot |\Omega | ~\ge ~ \dfrac{n}{\nu ^{*2}} \exp \left( -\dfrac{\rho _e^*}{n} \right).\) 

The result of Theorem REF  asserts that intrinsic tradeoffs emerges among all linear networks with similar estimation quality and parameter (REF ): reducing number of sampling locations must be compensated by increasing the average number of samples per subsystem and vice versa.

Numerical Simulations
Let us consider a linear dynamical network (REF ) that consists of \(n\)  subsystems, which are randomly and uniformly distributed over a square-shape spatial domain \([0,1] \times [0,1]\) . The Euclidean (spatial) distance between the subsystems \(i\)  and \(j\)  is denoted by \(\mathrm {dis}(i,j)\) . If two subsystems lie in each others connectivity range, then there will be a coupling between the two subsystems and the corresponding entries in the state space will be nonzero numbers. More precisely, the state matrix \(\mathbf {A}=[a_{ij}]\)  is defined by
\(a_{ij}=\left\lbrace \begin{array}{lr} \zeta _{ij} \hspace{1.42271pt}\mathrm {e}^{-\mathfrak {a} \hspace{1.42271pt}\mathrm {dis}(i,j)^{\mathfrak {b}}}& \mathrm {dis}(i,j) \le d\\ 0 & \mathrm {dis}(i,j)> d\end{array} \right..\) 

for some \(d >0\) . The parameters \(\mathfrak {a} >0\)  and \(0 < \mathfrak {b} <1\)  determine decay rate of the couplings and spatial localization properties of the network. For instance, larger values of \(\mathfrak {a}\)  and \(\mathfrak {b}\)  result in more localized networks with short range couplings. To make our study generic, the coefficients \(\zeta _{ij}\)  are independently and randomly chosen from \(\mathcal {N}(0,1)\) . In our simulations, we set \(n=40\) , \(\mathfrak {a}=1\) , \(\mathfrak {b}=0.5\) , and \(d=0.3\) . We generate and save one state matrix \(\mathbf {A}\)  that has both stable and unstable modes and use it in the following simulation studies. The spatial locations of subsystems and their coupling topology are illustrated in Fig. REF .
<FIGURE>
Constructing Observability Frame: We set \(\Omega =\lbrace 1,\dots ,n\rbrace \) , \(\tau =0.12\) , and \(M_i=M=44\)  and utilize Corollary REF  to construct an observability frame \(\Phi \)  with format (REF ). The resulting frame contains \(nM=1760\)  vectors with a space-time representation illustrated by blue dots in Fig. REF .
<FIGURE>For a noise intensity of \(\sigma =0.1\) , the value of the (least-squares) estimation measure is
\(\rho _d(\Phi ) \approx 0.0967.\) 
Randomized Frame Sparsification: For design parameters \(\epsilon =0.5\)  and \(q=590\) , we apply Theorem REF  and find a sparsified frame with \(|\Phi _s|=503\) .
The spatial locations and time stamps of the sampling strategy corresponding to the sparsified frame are illustrated by black circles in Fig. REF . The value of the estimation measure is
\(\rho _d(\Phi _s)\approx 0.1883.\) 
The number of space-time samples has been reduced by \(71\%\)  for the price of \(95\%\)  relative estimation-quality loss.

Sparsification via Random Partitioning: Let us reconsider the frame \(\Phi \)  shown in Fig. REF . By applying Proposition REF , random partitioning of \(\Phi \)  leads into two subsets \(\Phi _1\)  and \(\Phi _2\) . In our simulations, we repeat the random partitioning procedure \(5\times 10^5\)  times. In each case, the minimum and maximum relative estimation-quality losses for \(\Phi _1\)  and \(\Phi _2\) , i.e.,
\(\max _{j=1,2}\dfrac{\rho _d(\Phi _j)-\rho _d(\Phi )}{\rho _d(\Phi )}~~\textrm {and}~~\min _{j=1,2}\dfrac{\rho _d(\Phi _j)-\rho _d(\Phi )}{\rho _d(\Phi )}\) 
is computed and saved. The histogram of this data is depicted in Fig. REF . In this simulation, the minimum and maximum degradations are less than \(0.55\)  with a high probability. The theoretical estimate from Theorem REF  is
\(r^*=\max _{\phi \in \Phi }~r_{\phi }(\mathbf {S}) \hspace{1.42271pt}\approx \hspace{1.42271pt}0.0315\Rightarrow \kappa (r^*)\approx 1.1441.\) 
Our extensive simulations reveal that, in practice, one typically achieves comparably better estimation quality than our theoretical bounds (REF ). In these simulations, the number of space-time samples are reduced by almost \(50\%\)  for the price of \(55\%\)  estimation-quality loss (in most outcomes of the simulations).
<FIGURE><FIGURE>
Performance of Randomized vs. Greedy Algorithms: In this simulation, we compare the estimation quality of the resulting sparsified frames from Algorithm REF  and REF . Using Algorithm REF , we construct 25 different sparsified frames by selecting 25 different values for \(\epsilon \)  in \((1/\sqrt{n},1]\) . We treat \(q\)  as a control parameter and vary its value between 5902 and 153. For a fixed \(q\) , we compute the value of the estimation measure for all 25 frames and save the one with the minimum value. When applying Algorithm REF , we change the desired sparsity level \(\mathfrak {s}\)  to get a sequence of sparsified frames. The outcome of our simulations is depicted in Fig. REF , where one can observe that both methods result in almost similar estimation qualities. The only difference we can report is their running time (on a personal computer with an Intel processor using MATLAB): the randomized method (including 25 experiments per \(q\) ) took about \(1.68\)  seconds, while the greedy method took \(26.71\)  seconds. This is consistent with our running time analysis for both algorithms.
<FIGURE>
Sequential Frame Construction:
We compute the value of \(\delta ^*\) , which is defined in Theorem REF , using the saved state matrix \(\mathbf {A}\)  and get \(\delta ^*\approx 0.0434\) . We merge \(N=12\)  subframes in the time horizon, where \(c_j=0.25\)  is for every subframe \(\Phi _j\)  for \(j=1,\dots ,N\) . Sampling times \(t_{ij}\)  for each frame is chosen randomly and uniformly from time interval \([t_j,t_{j+1}]\) . The space-time representation of these frames is illustrated in Fig. REF . The resulting concatenated frame has \(|\Phi |=Nn=480\)  components with estimation quality
\(\rho _d(\Phi )\approx 0.1862.\) 
Estimation Quality Deterioration with Time Shifts: Let us consider the first observability frame \(\Phi _1\)  in Fig. REF  with \(n\)  components. First, we construct a family of observability frame \(\Phi _\delta \) , which is defined in (REF ), by increasing \(\delta \)  from 0 to \(0.5\) . We compute exact value of the (least-squares) estimation measure for every \(\Phi _\delta \) . The result of our simulations is depicted in Fig. REF  along with our theoretical upper bound (REF ). One observes that our proposed upper bound is rather tight for all values of \(\delta \)  in \([0,0.5]\) .
<FIGURE>
Discussion and Conclusion
In this paper, we assume that measurement noises in (REF ) are Gaussian and independent of each other. When the measurement noises are dependent, the covariance of the estimation error () will be
\(\mathbf {\Sigma }_{\eta }= \left(T \mathbf {\Sigma }_{\xi }^{-1} )^{-1}\right.where \mathbf {\Sigma }_{\xi } is the covariance of the measurement noise. In general, we may not be able to expand (\ref {covar-eta}) as a sum of rank-one matrices made of frame components. This was a useful property for our developments in Section \ref {sec:sparse}. This case needs a thorough analysis which is beyond the scope of this paper. {\color {black}However, in the case that \mathbf {\Sigma }_{\xi } is diagonal and with entries \mathbf {\Sigma }_{\xi }=\mathrm {diag}(\sigma _{i,t}^2)_{(i,t) \in \mathfrak {S}}, i.e., the observations are spatially and temporally independent, one can conduct a similar analysis bydefining the frame components as\Phi (\mathbf {A},\mathfrak {S}) \hspace{1.42271pt}=\hspace{1.42271pt}\left( \sigma _{i,t}^{-1}\, \mathrm {e}^{\mathbf {A}^Tt} e_i ~\big |\ ~ (i,t) \in \mathfrak {S} \right),where \sigma _{i,t}^2 is the variance of the observation error at (i,t).}\) 
Our results can be extended to include linear dynamical networks with arbitrary output matrices. In such networks, the sampling locations will be different from subsystem locations and the components of the observability frame (REF ) will take form \(\mathrm {e}^{\mathbf {A}^T t} c^T\) , where \(c\)  is a row of output matrix. The main reason for working with output matrix (REF ) is to highlight inherent tradeoffs between minimum required number of samples in space and time in order to achieve a certain estimation quality in linear dynamical networks; see the results of Section .
In order to quantify estimation quality, we consider two estimation measures: standard deviation and differential entropy of the estimation error. Depending on the specific design criteria, one may choose another type of estimation measure. For instance, if the reliability of the estimators is of significant importance for the network designer, one can utilize risk measures to evaluate the estimation quality. Two useful risk measures are: risk of large aggregate deviations, i.e.,
\(\rho _{a} (\Phi ):=\inf \Big \lbrace \Delta \in \mathbb {R}_{+} ~\Big | ~\mathbb {P} \big \lbrace  \Vert \eta \Vert  \ge \Delta \big \rbrace  \hspace{1.42271pt}\le \hspace{1.42271pt}\epsilon \Big \rbrace \) 
or risk of large element-wise deviations, i.e.,
\(\rho _{r}(\Phi ) :=\inf \Big \lbrace \Delta \in \mathbb {R}_{+} ~\Big | ~\mathbb {P} \big \lbrace  |\eta _i| \ge \Delta \big \rbrace  \hspace{1.42271pt}\le \hspace{1.42271pt}\epsilon \Big \rbrace .\) 

One can show that these two risk measures are monotone according to Definition REF  and as a result, one can effectively employ these measures instead of \(\rho _d\)  and \(\rho _e\)  in our proposed methodology.
The significance of fundamental limits and tradeoffs in Section  is that they reveal what is achievable and what is not. This is practically plausible as it prevents us from searching for sampling strategies with unachievable estimation qualities.
The results of Section  provide three methods to sparsify a given observability frame. Our theoretical error bounds are rather conservative. However, our extensive simulations assert that our proposed algorithms can achieve comparably better error bounds in practice.

Appendix I: Proof of The Remaining Results
 Proof of Theorem REF :
Let us assume that the family of vectors (REF ) is a frame for \(\mathbb {R}^n\)  with analysis operator \({\bf T}\) . From Lemma REF , initial condition can be recovered via
\(x_0 = {\bf T}^{\dagger } y,\) 
where \(y=\left[\left\langle x , \phi _i \right\rangle \right]_{\phi _i \in \Phi }\) . In the next step, suppose that \((\ref {eq:phidef_general_thm})\)  is not a frame for \(\mathbb {R}^n\) . Then, the frame matrix \({\bf S}\)  will be singular by Proposition REF . Thus, for every \(x_1\in \mathbb {R}^n\) , two initial states \(x_1\)  and \(x_1+x_2 \) , in which \(x_2\)  is a nonzero element of \(\mathrm {null}({\bf S})\) , produce the same vector of observation \(y\) . This is contrary to our assumption on unique determination of the initial state of the network.
Proof of Lemma REF :
Based on the Cayley-Hamilton theorem
\(\mathrm {e}^{\mathbf {A}t}=\sum _{k=0}^{n-1} g_k (t) \mathbf {A}^k,\) 

for some functions \(g_k(t)\)  for \(k=0,\dots ,n-1\) . Using this fact, the analysis matrix corresponding to the sampling strategy \(\mathfrak {S}\)  can be written as
\( \mathbf {T}=\left[\sum _{k=0}^{n-1} g_k (t) \mathbf {A}^k e_i \right]_{(i,t) \in \mathfrak {S}}.\) 

Assume that \(\Phi \)  is a frame but the pair \((\mathbf {A},\mathbf {C}_\Omega )\)  is not observable. Thus, the observability matrix
\(O(\mathbf {A},\mathbf {C}_\Omega )=\left[ \big (\mathbf {A}^k \mathbf {C}_\Omega ^T\big )^T\right]_{k=0,\dots ,{n-1}}\) 

has rank less than \(n\) . Because \(\Phi \)  is a frame, rank of \( is \) n\(. However, comparing (\ref {eq:T_CH}) with the form of the observability matrix, we observe that rank of matrix \)  in (REF ) is at most equal to rank of \(O(\mathbf {A},\mathbf {C}_\Omega )\) . This is a contradiction, proving that \((\mathbf {A},\mathbf {C}_\Omega )\)  must be observable.

Proof of Proposition REF :
The expression for the estimation measure holds because we can write
\(\rho _d(\Phi )^2 &=\mathbb {E} \lbrace  \Vert  \eta \Vert _2^2\rbrace  = \mathbb {E} \Big \lbrace \xi ^T {\bf T} {\bf S}^{-2} {\bf T}^T \xi \Big \rbrace =\mathbb {E} \Big \lbrace  {\rm Tr} \big ( {\bf S}^{-2} {\bf T}^T \xi \xi ^T {\bf T}\big )\Big \rbrace \\& = \sigma ^2 {\rm Tr} ({\bf S}^{-1})=\sigma ^2\sum \limits _{i=1}^n \dfrac{1}{\lambda _i(\mathbf {S})}.\) 

To see that \(\rho _d\)  is a monotone operator, consider the following chain of observations:
\(\Phi _1 \subseteq \Phi _2 \Rightarrow \mathbf {S}_1 \preceq \mathbf {S}_2 \Rightarrow \mathbf {S}_2^{-1} \preceq \mathbf {S}_1^{-1} \Rightarrow \mathrm {Tr}(\mathbf {S}_2^{-1} )\le \mathrm {Tr}( \mathbf {S}_1^{-1}).\) 
The last inequality implies that \(\rho _d(\Phi _2) \le \rho _d(\Phi _1)\) .
Proof of Proposition REF :
For the random variable \(\eta \) , using (REF ), we can show that (see Chapter 8 in [16]})
\({h}(\eta )&=~ \dfrac{1}{2}\log \left( \mathrm {det} (\sigma ^2 \mathbf {S}^{-1}) \right)+\dfrac{n}{2}\log (2\pi e) \nonumber \\ &= ~\dfrac{1}{2} \sum _{i=1}^n \log \left(\dfrac{\sigma ^2}{\lambda _i(\mathbf {S})}\right) +\dfrac{n}{2}\log (2\pi e) \\&= \frac{1}{2} \rho _e(\Phi ) + \frac{n}{2} \big (1+ \log (2 \pi \sigma ^2) \big ).\) 

To see that \(\rho _e\)  is a monotone operator, consider the following chain of operations:
\(\Phi _1 \subseteq \Phi _2& \Rightarrow & \mathbf {S}_1 \preceq \mathbf {S}_2 \\& \Rightarrow & \mathbf {S}_2^{-1} \preceq \mathbf {S}_1^{-1} \\& \Rightarrow & \log (\mathbf {S}_2^{-1} )\preceq \log ( \mathbf {S}_1^{-1}) \\& \Rightarrow & \mathrm {Tr}( \log (\mathbf {S}_2^{-1} ) )\le \mathrm {Tr}( \log (\mathbf {S}_1^{-1} )).\) 
The third inequality holds because \(\log \) (.), as a map from the cone of positive definite matrices to the set of symmetric matrices, is analytic and increasing on the cone of positive definite matrices. The last inequality implies that \(\rho _e(\Phi _2) \le \rho _e(\Phi _1)\) .
Proof of Proposition REF :
Let us denote the analysis matrices corresponding to \(\Phi \)  and \(\Phi _\delta \)  by \( and \)\(, respectively. It follows that\)\(\delta =^{\mathbf {A}\delta }.\)\(Thus, the corresponding frame matrix for \)\( is\)\(\mathbf {S}_\delta =\delta ^T\delta =\mathrm {e}^{\mathbf {A}^T \delta }T^{\mathbf {A}\delta }=\mathrm {e}^{\mathbf {A}^T \delta }\mathbf {S}\mathrm {e}^{\mathbf {A}\delta },\)\(where \) S\( is the frame matrix of \)\(.As a result, we can see that\)\(\mathbf {S}_\delta ^{-1}= \mathrm {e}^{-\mathbf {A}\delta } \mathbf {S}^{-1} \mathrm {e}^{-\mathbf {A}^T\delta } .\)\(Therefore, estimation using the shifted frame \)\( results in a normal error (random) variable{\begin{@align}{1}{-1}\eta \sim \mathcal {N}\left(0,\sigma ^2 \mathbf {S}_{\delta }^{-1} \right)=\mathcal {N}\left(0,\sigma ^2 \mathrm {e}^{-\mathbf {A}\delta } \mathbf {S}^{-1} \mathrm {e}^{-\mathbf {A}^T\delta } \right).\end{@align}}In the next step, the estimation measures can be found as follows\begin{eqnarray*}\frac{1}{\sigma ^2} \hspace{1.42271pt}\rho _d(\Phi _\delta )^2 &=& \mathrm {Tr} \left( \mathrm {e}^{-\mathbf {A}\delta }\mathbf {S}^{-1} \mathrm {e}^{-\mathbf {A}^T \delta } \right)\\& = & \mathrm {Tr} \left( \mathrm {e}^{-\mathbf {A}^T \delta } \mathrm {e}^{-\mathbf {A}\delta }\mathbf {S}^{-1} \right) \\& = & \mathrm {Tr}\left( \big (\mathrm {e}^{\mathbf {A}\delta } \mathrm {e}^{\mathbf {A}^T \delta }\big )^{-1} \mathbf {S}^{-1}\right) \\& \le &\sum _{i=1}^n \sigma _i\left( \big (\mathrm {e}^{\mathbf {A}\delta } \mathrm {e}^{\mathbf {A}^T \delta }\big )^{-1} \right) \hspace{1.42271pt}\sigma _i(\mathbf {S}^{-1}),\end{eqnarray*}where in the last equality we have used Von Neumann^{\prime }s trace inequality; we refer to \cite {mirsky1975trace} for more details.We can further write\begin{eqnarray*}\rho _d(\Phi _\delta )^2 & \le & \sigma ^2 \sum _{i=1}^n \sigma _i^2 \big ( \mathrm {e}^{-\mathbf {A}\delta } \big ) \hspace{1.42271pt}\lambda _i(\mathbf {S}^{-1})\\& =& \sigma ^2 \sum _{i=1}^n \dfrac{1}{\sigma _i^2(\mathrm {e}^{\mathbf {A}\delta }) \cdot \lambda _i(\mathbf {S})}.\end{eqnarray*}For the last part of the proof, we have\begin{eqnarray*}\rho _e(\Phi _\delta ) & = & -\log (\det (\mathbf {S}_\delta ))=-\log \big (\det (\mathrm {e}^{\mathbf {A}^T \delta }\mathbf {S}\mathrm {e}^{\mathbf {A}\delta })\big )\\& = & -\log \big (\det (\mathrm {e}^{\mathbf {A}^T \delta })\big ) \det \big ((\mathrm {e}^{\mathbf {A}\delta })\big )-\log \det (\mathbf {S}) \\& = & \rho _e(\Phi )-\sum _{i=1}^n \log \big (\lambda _i(\mathrm {e}^{\mathbf {A}^T \delta }\mathrm {e}^{\mathbf {A}\delta })\big ) \\& = & \rho _e(\Phi )-\sum _{i=1}^n \log \big (\sigma _i^2(\mathrm {e}^{\mathbf {A}\delta })\big ).\end{eqnarray*}\)
Proof of Theorem REF : 
When \({\bf A}={\bf 0}\) , the proof becomes trivial as in this case \(\mathrm {e}^{{\bf A}t} =\bf I\)  for all \(t\in \mathbb {R}\)  and \(\mathbf {C}_\Omega \)  must be the identity matrix according to Assumption REF . When \({\bf A}\)  is nonzero, its minimal polynomial has degree \(\mathfrak {d}(\mathbf {A})\)  with \(1\le \mathfrak {d}(\mathbf {A})\le n\) . Then, it follows that
\(\mathrm {e}^{\mathbf {A}t} \hspace{1.42271pt}=\hspace{1.42271pt}\sum _{k=0}^\infty \hspace{1.42271pt}\frac{t^k}{k!} \hspace{1.42271pt}{\bf A}^k \hspace{1.42271pt}=\hspace{1.42271pt}\sum _{k=0}^{\mathfrak {d}(\mathbf {A})-1} \hspace{1.42271pt}g_k (t) \hspace{1.42271pt}{\bf A}^k\) 
for all \(t\in \mathbb {R}\) , where \(g_k\) 's are some functions of time.
Let \({\bf J}\)  be the Jordan canonical form of \({\bf A}\)  and
\({\bf A}= {\bf Q}^{-1} {\bf J} {\bf Q}.\) 
for some nonsingular matrix \({\bf Q}\) . By combining (REF ) and (REF ), one obtains
\(\mathrm {e}^{{{\bf J}t}}=\sum _{k=0}^{\mathfrak {d}(\mathbf {A})-1} \hspace{1.42271pt}g_k (t) \hspace{1.42271pt}{\bf J}^k\) 
for all \(t\in \mathbb {R}\) . Let us consider the minimal polynomial of \(\mathbf {A}\)  given by (REF ), where \(p_m\) 's are some integer numbers for \(1\le m\le q\) .
From the definition of a minimal polynomial, there are Jordan blocks
\({\bf J}_{\lambda _m, p_m}\) 
associated with every eigenvalue \(\lambda _m(\mathbf {A})\)  whose dimension is \(p_m\) .
This together with (REF ) implies that
\( \mathrm {e}^{{{\bf J}_{\lambda _m, p_m} t}} \hspace{1.42271pt}= \hspace{1.42271pt}\sum _{k=0}^{\mathfrak {d}(\mathbf {A})-1} \hspace{1.42271pt}g_k (t) \hspace{1.42271pt}\left({\bf J}_{\lambda _m, p_m} \right)^k,\) 
or equivalently,
\(\mathrm {e}^{\lambda _m(\mathbf {A}) t} \hspace{1.42271pt}\frac{ t^l}{l!} \hspace{1.42271pt}= \hspace{1.42271pt}\sum _{k=l}^{\mathfrak {d}(\mathbf {A})-1} \hspace{1.42271pt}g_k (t) \hspace{1.42271pt}{{k}\atopwithdelims (){l}} \hspace{1.42271pt}\lambda _m(\mathbf {A})^{k-l}\) 
for all \(0 \le l \le p_j-1\)  and \(t\in \mathbb {R}\) .
The last equivalence holds as \((j, j^{\prime })\) 'th entry of \(({\bf J}_{\lambda _m, p_m})^k\)  with property \(0\le j^{\prime }-j\le \min (k, p_m-1)\)  equals to \({{k}\atopwithdelims (){j^{\prime }-j}} \lambda _m(\mathbf {A})^{k-j+j^{\prime }}\)  and all other entries are equal to zero.
It is well known that functions \(\mathrm {e}^{\lambda _m(\mathbf {A}) t} t^k\)  for \(m=1,\dots ,q\)  and \(k=0,\dots , p_m-1\)  for \(t \in \mathbb {R}\)  are linearly independent.
This, together with (REF ), and the fact that
\(\mathfrak {d}(\mathbf {A})=\sum _{m=1}^q p_m\) 
implies the existence
of a nonsingular matrix \({\bf D}\)  that satisfies
\(G(t)= {\bf D} E(t),\) 
where \(E(t)\)  is defined in (REF )
and
\(G(t)=\big [g_k(t)\big ]_{0\le k\le \mathfrak {d}(\mathbf {A})-1}.\) 
Let us denote \(\Theta _i=\big \lbrace t_{i,j}~\big |~ 1\le j\le M_i\big \rbrace \) . According to our assumptions, matrix \({\bf E_i}\)  defined by (REF ) has full row rank \(\mathfrak {d}(\mathbf {A})\) .
This, together with (REF ) and (REF ), implies the existence of scalars
\(a_{i, j, k}\)  such that
\({\bf A}^k~=~ \sum _{j=1}^{M_i} ~a_{i, j, k}~ \mathrm {e}^{ {\bf A}^Tt_{i, j}}\) 
for \(0\le k\le \mathfrak {d}(\mathbf {A})-1\) . Let us denote
\( \mathbf {C}_\Omega = \left[\begin{array}{ccc}e_{i_1}~\big | & \dots &\big |~e_{i_{p}}\end{array}\right]^T\) 
and
\({\bf F}=\left[ ({\bf A}^T)^k \hspace{1.42271pt}e_i \right]_{i\in \Omega \atop 0\le k\le \mathfrak {d}(\mathbf {A})-1}.\) 
Based on Assumption REF ,
and the fact that \(\mathfrak {d}(\mathbf {A})\)  is greater than or equal to the observability index of \((\mathbf {A},\mathbf {C})\) , for \(n \times \mathfrak {d}(\mathbf {A})|\Omega |\)  matrix \({\bf F}\)  it holds that
\({\rm rank}({\bf F})=n,\) 
(for example see Section 6.3.1 in [17]}).
Let us define
\(\left[\mathrm {e}^{{\bf A}^Tt} \hspace{1.42271pt}e_i\right]^T_{ {i\in \Omega },~{t\in \Theta _i}}.\) 
From (REF ), it follows that the
rank of matrix \(, whose size is \) n|S|\(, islarger than or equal to rank of \) [CAk ]0k d(A)-1\(. This together with (\ref {appendixA.eq8})implies that \)  has rank \(n\) , which implies that the family of vectors (REF ) form a frame for \(\mathbb {R}^n\) .

Proof of Corollary REF :  First suppose that \(|\Theta _i|=\mathfrak {d}(\mathbf {A})\)  for each \(i \in \Omega \) . Observe that nonzero combinations of \(\mathrm {e}^{\lambda _i(\mathbf {A}) t} t^k\)  for \(i=1,\dots ,q\)  and \(k=0,\dots , p_i-1\)  have only finitely many zeros. This shows that for any \(\tau >0\)  and \(i \in \Omega \)
\(\det ({\mathbf {E}_i})=\mathrm {det}([E(t_j)]_{j=1,\dots ,\mathfrak {d}(\mathbf {A})}) \ne 0,\) 
for almost every choice of times \([t_1, \ldots , t_{\mathfrak {d}(\mathbf {A})}]\in [0, \tau ]^{\mathfrak {d}(\mathbf {A})}\) . Hence, if
we choose the sampling times \(t_j \in \Theta _i\)  with \(|\Theta _i|=\mathfrak {d}(\mathbf {A})\)  randomly in the range \([0, \tau ]\)  in an independent and uniform manner, then with probability one the requirement (REF ) is satisfied and \(\mathbf {E}_i\)  is full rank for each \(i \in \Omega \) . Thus, by Theorem REF , the resulting family of vectors is a frame with probability one. If the number of random samples per location \(i \in \Omega \)  increases beyond \(\mathfrak {d}(\mathbf {A})\) , the result is still a frame.

Proof of Theorem REF : First assume that for each location \(i\in \Omega \) , we choose \(|\Theta _i|=\mathfrak {d}(\mathbf {A})\) . In this case, for every location \(i \in \Omega \)
\({\bf E}_i= [E(k\delta )]_{k=0,\dots , \mathfrak {d}(\mathbf {A})-1}.\) 
one observes that
\(\det (\mathbf {E}_i)= C \prod _{1\le m<m^{\prime }\le q} \left(\mathrm {e}^{\lambda _m(\mathbf {A}) \delta }-\mathrm {e}^{\lambda _{m^{\prime }}(\mathbf {A})\delta } \right)^{p_mp_{m^{\prime }}}\) 
for some nonzero number \(C\)  depending only on \(p_m\)  with \(1\le m\le q\) . So the requirement (REF ) is satisfied if
\((\lambda _m(\mathbf {A})-\lambda _{m^{\prime }}(\mathbf {A}))\delta \notin 2\pi \mathrm {j}\, \mathbb {Z}\) 
for all distinct eigenvalues \(\lambda _m(\mathbf {A})\)  and \( \lambda _{m^{\prime }}(\mathbf {A})\) . Therefore by Theorem REF , for \(\mathbf {B}_\delta :=\mathrm {e}^{\mathbf {A}^T \delta }\) , we know that
\({\Phi }^{\prime }:= \left(\left.{{\bf B}_\delta ^{k} e_i} \right|\ i \in \Omega ,~ k=0,\dots ,\mathfrak {d}(\mathbf {A})-1 \right).\) 
is a frame. The family of vectors \(\Phi \)  given in the theorem satisfies \({\Phi }^{\prime } \subseteq \Phi \) . Thus, \(\Phi \)  is also a frame.
Proof of Theorem REF :
Since matrix \(\mathrm {e}^{\mathbf {A}^Tt^*}\)  is full rank for all \(t^*\ge 0\) ,
one can verify that \({\Phi }\)  in (REF ) forms a
frame for \(\mathbb {R}^n\)  if and only if \(\big ({\mathrm {e}^{\mathbf {A}^T(t-t^*)} e_i}~\big |~(i,t) \in \mathfrak {S} \big )\) 
forms a
frame. Due to this shift-invariance property, we may safely assume that \(t^*=0\)  by shifting every element in all \(\Theta _i\) 's by \(-t^*\) . Let us pick \(t_i\in \Theta _i\)  for \(1\le i\le n\) . Then, for every vector \(c=[c_1, \ldots , c_n]^T\) , we have
\(\left\Vert \sum _{i=1}^n c_i e^{\mathbf {A}t_i} e_i- \sum _{i=1}^n c_i e_i\right\Vert &= & \left\Vert  \sum _{m=1}^\infty \frac{\mathbf {A}^m}{m!} \sum _{i=1}^n c_i t_i^m e_i\right\Vert \\&\le & \sum _{m=1}^\infty \frac{\Vert \mathbf {A}\Vert ^m}{m!} ~\left\Vert \sum _{i=1}^n c_i t_i^m e_i\right\Vert \\&\le & \sum _{m=1}^\infty \frac{\Vert \mathbf {A}\Vert ^m (\delta ^*)^m}{m!} \Vert c\Vert  \\& = & \left(e^{\Vert \mathbf {A}\Vert  \delta ^*}-1 \right)\hspace{1.42271pt}\Vert c\Vert .\) 
This implies that
\(\left(2-e^{\Vert \mathbf {A}\Vert  \delta ^*} \right) \hspace{1.42271pt}\Vert c\Vert  \hspace{1.42271pt}\le \hspace{1.42271pt}\Big \Vert \sum _{i=1}^n c_i e^{\mathbf {A}t_i} e_i\Big \Vert  \hspace{1.42271pt}\le \hspace{1.42271pt}e^{\Vert \mathbf {A}\Vert  \delta ^*}\Vert c\Vert \) 
for all \(c\in \mathbb {R}^n\) . Hence, \(\Phi \) , with \(|\Theta _i|=1\)  for every sampling location \(i \in \Omega \) , consists of \(n\)  linearly independent vectors and forms a frame for \(\mathbb {R}^n\) . By increasing the number of samples per sampling location, \(\Phi \)  will remain to be a frame for \(\mathbb {R}^n\) .
Proof of Corollary REF : First, for each \(i \in \Omega \) , consider \(|\Theta _i|=1\)  and denote \(\Theta _i=\lbrace t_i\rbrace \) . Moreover, we define
\(\theta =[t_i]_{i=1,\dots ,n}.\) 

Now, we denote \(\mathbf {S}(\theta )\)  to be the frame matrix corresponding to family of vectors \((\ref {eq:nonuniform_complicated})\) , where the sampling times have been gathered in \(\theta \in \mathbb {R}^n\) . Based on Theorem REF , there exist a \(\theta \)  for which the real analytic function \(F(\theta ):~\mathbb {R}^n \rightarrow \mathbb {R}\)  defined by
\(F(\theta ):=\mathrm {det}(\mathbf {S}(\theta )),\) 
is nonzero. Thus, measure of points \(\theta \in [0,\tau ]^{n}\)  for which \(F(\theta )\)  vanishes is zero (e.g. see [18]}), and independent and uniform sampling of the time stamps in \([0,\tau ]^{n}\)  gives a frame with probability one. If we increase \(\Theta _i\)  beyond 1 probability of getting a frame is 1.
Proof of Theorem REF :
We build up our proof based on some of the steps taken in the proof of Theorem 5 in [19]}. Using the leverage scores (REF ), one can verify that
\({\bf 0} \hspace{1.42271pt}\prec \hspace{1.42271pt}(1-\epsilon ) \hspace{1.42271pt}{\bf S} \hspace{1.42271pt}\preceq \hspace{1.42271pt}\mathbf {S}_w\) 

with probability at least \(1/2\) , where
\( \mathbf {S}_w \hspace{1.42271pt}= \hspace{1.42271pt}\sum _{\phi \in \Phi }\hspace{1.42271pt}w_s(\phi ) \hspace{1.42271pt}\phi \phi ^T \hspace{1.42271pt}= \hspace{1.42271pt}\sum _{\phi \in \Phi } \hspace{1.42271pt}\frac{f_{\phi }}{q\pi (\phi )} \hspace{1.42271pt}\hspace{1.42271pt}\phi \phi ^T. \) 
and \(f_{\phi }\ge 0\)  is the frequency of the times that \(\phi \)  is sampled (due to sampling with replacement, some vectors can be sampled multiple times). Weights of those \(\phi \notin \Phi _s\)  are equal to 0. For an outcome of Algorithm REF , one has
\( {\bf 0} \hspace{1.42271pt}\prec \hspace{1.42271pt}{\bf S}_w \hspace{1.42271pt}\preceq \hspace{1.42271pt}\left( \max _{\phi \in \Phi _s} \hspace{1.42271pt}w_s(\phi ) \right) \hspace{1.42271pt}\sum _{\phi \in \Phi _s} \phi \phi ^T \) 
This implies that \({\bf S}_s = \sum _{\phi \in \Phi _s} \phi \phi ^T \succ {\bf 0}\)  with probability at least \(1/2\) . For the first part of our proof, \(\Phi _s\)  is a frame for \(\mathbb {R}^n\)  if and only if \(\mathbf {S}_s \succ {\bf 0}\) .
Denote the analysis operator of \(\Phi _s\)  by \(s\)  and set \({\bf W}_s=\mathrm {diag}\big (w_s(\phi )\big )\big |_{\phi \in \Phi _s}\) .
Let
\(y_s=s x_0+\xi _s\) 

be a noisy observation vector collected by \(\phi \in \Phi _s\) , in which
\(\xi _s \in \mathbb {R}^m\)  is a zero mean Gaussian measurement noise with independent components and covariance
\(\mathbb {E} \left\lbrace  \xi \xi ^T \right\rbrace =\sigma ^2 I\) . In the next step, let us consider an alternative estimator \(\tilde{x}_0\)  that is given by
\(\tilde{x}_0 & = & \big (s^T {\bf W}_s \hspace{1.42271pt}s\big )^{-1} \hspace{1.42271pt}s^T {\bf W}_s y_s.\) 
It is straightforward to verify that this is an unbiased estimator using only the observations corresponding to \(\Phi _s\) . Since covariance of noise is \(\mathbb {E} \left\lbrace  \xi \xi ^T \right\rbrace =\sigma ^2 I\) , the unweighted least-squares estimator gives the optimal estimator \(\hat{x}_0\) . Therefore,
\(\mathbb {E}\left\lbrace \tilde{x}_0 \tilde{x}_0^T\right\rbrace  \succeq \mathbb {E}\left\lbrace \hat{x}_0 \hat{x}_0^T \right\rbrace .\) 
This lets us write
\( \rho _d(\Phi _s)^2&=& \sigma ^2 \hspace{1.42271pt}\mathrm {Tr}\left(\mathbb {E}\lbrace \hat{x}_0 \hat{x}_0^T\rbrace \right) \\&\le &\sigma ^2 \hspace{1.42271pt}\mathrm {Tr}\left(\mathbb {E}\lbrace \tilde{x}_0 \tilde{x}_0^T\rbrace \right) \\ & = & \sigma ^2 \hspace{1.42271pt}\mathrm {Tr}\left(\mathbf {S}_w^{-1} \hspace{1.42271pt}T {\bf W}_s^2 \hspace{1.42271pt}\mathbf {S}_w^{-1} \right) \\& = &\sigma ^2 \hspace{1.42271pt}\mathrm {Tr}\left(\mathbf {S}_w^{-1} \hspace{1.42271pt}\hat{\mathbf {S}}_w \hspace{1.42271pt}\mathbf {S}_w^{-1} \right) \) 
where \(\mathbf {S}_w = s^T {\bf W}_s s\)  and \(\hat{\mathbf {S}}_w := s^T {\bf W}_s^2 \hspace{1.42271pt}s\) . From (REF ) and the definition of random variable \(\chi \)  in (REF ), it follows that
\(\rho _d(\Phi _s)^2 ~\le ~ \sigma ^2 \chi \hspace{1.42271pt}\mathrm {Tr}\big (\mathbf {S}_w^{-1} \mathbf {S}_w \mathbf {S}_w^{-1} \big )\) 

in which the middle matrix is replaced by its upper bound as one can show that for every three positive-definite matrices \(\mathbf {X}_1,\mathbf {X}_2,\mathbf {X}_3\)  with \(\mathbf {X}_1 \succeq \mathbf {X}_2\) , inequality
\(\mathbf {X}_3 \mathbf {X}_1 \mathbf {X_3} \succeq \mathbf {X}_3 \mathbf {X}_2 \mathbf {X_3}\) 
holds. According to Markov inequality, the next inequality holds
\(\chi \hspace{1.42271pt}\le \hspace{1.42271pt}\dfrac{1}{1-\frac{3}{4}} \hspace{1.42271pt}\mathbb {E}\lbrace \chi \rbrace  \hspace{1.42271pt}= \hspace{1.42271pt}4\bar{\chi }\) 

with probability at least \(3/4\) . From (REF ) and (REF ), we get
\(\rho _d(\Phi _s) ^2 ~\le ~4 \hspace{1.42271pt}\sigma ^2 \hspace{1.42271pt}\bar{\chi }\hspace{1.42271pt}\hspace{1.42271pt}\mathrm {Tr} \big (\mathbf {S}_w^{-1}\big ).\) 

On the other hand, by applying similar steps to the proof of Theorem 5 in [19]}, it follows that with probability at least \(1/2\)  we have
\((1-\epsilon ) \mathbf {S}\hspace{1.42271pt}\preceq \hspace{1.42271pt}\mathbf {S}_w \hspace{1.42271pt}\preceq \hspace{1.42271pt}(1+\epsilon ) \mathbf {S}.\) 

By taking inverse, we get
\((1+\epsilon )^{-1} \mathbf {S}^{-1} \hspace{1.42271pt}\preceq \hspace{1.42271pt}\mathbf {S}_w^{-1} \hspace{1.42271pt}\preceq \hspace{1.42271pt}(1-\epsilon )^{-1} \mathbf {S}^{-1}.\) 

If the event described in (REF ) is denoted by \(\mathfrak {A}\)  and the event described by (REF ) is denoted by \(\mathfrak {B}\) , then
\(\mathbb {P}(\mathfrak {A} \cap \mathfrak {B}) = \mathbb {P}(\mathfrak {A})+\mathbb {P}(\mathfrak {B})-\mathbb {P}(\mathfrak {A} \cup \mathfrak {B}) \ge \dfrac{3}{4}+\dfrac{1}{2} -1=\dfrac{1}{4}.\) 
Therefore both (REF ) and (REF ) hold with probability at least \(1/4\) .
Combining (REF ) and (REF ), one arrives at
\(\rho _d(\Phi _s)^2 \hspace{1.42271pt}\le \hspace{1.42271pt}\dfrac{4 \sigma ^2\bar{\chi }}{1-\epsilon } \hspace{1.42271pt}\hspace{1.42271pt}\mathrm {Tr}(\mathbf {S}^{-1})= \dfrac{4 \bar{\chi }}{1-\epsilon } \hspace{1.42271pt}\rho _d(\Phi ) ^2.\) 

Taking the square root from both sides, we get the desired inequality (REF ).
For the entropy estimation measure, by following almost identical steps, we can show that with probability at least \(1/4\)  the following inequality holds
\(\rho _e(\Phi _s) &= &\det (\log (\mathbf {S}_s^{-1})) \\& \le & \log \left(\det \left(\dfrac{4\bar{\chi }}{1-\epsilon } \mathbf {S}^{-1}\right) \right) \\&=& n\log \left( \dfrac{4\bar{\chi }}{1-\epsilon } \right) +\rho _e(\Phi ).\) 
Proof of Corollary REF : We can write
\(\sum _{\phi \in \Phi } w_s(\phi )^2 \hspace{1.42271pt}\phi \phi ^T & \le & \sum _{\phi \in \Phi } w_s(\phi ) \Big ( \max _{j=1,\dots ,|\Phi |} w_s(\phi )\Big )\, \phi \phi ^T \\& = & \Big ( \max _{j=1,\dots ,|\Phi |} w_s(\phi )\Big ) \sum _{\phi \in \Phi } w_s(\phi ) \hspace{1.42271pt}\phi \phi ^T.\) 
Therefore, \(\max _{\phi \in \Phi }\hspace{1.42271pt}w_s(\phi )\)  is an upper-bound on \(\chi \) . The rest of the proof follows from the fact that we can replace \(\chi \)  with \(\max _{\phi \in \Phi }\hspace{1.42271pt}w_s(\phi )\)  in the proof of Theorem REF .
Proof of Theorem REF :
Let us denote \(\mathbf {S}_j\)  to be the frame matrix corresponding to family of vectors \(\Phi _j\)  that is resulted from the random partitioning according to Proposition REF . For the least-squares estimation measure, we have
\(\rho _d(\Phi _j) &=& \sigma \sqrt{\mathrm {Tr} \big (\mathbf {S}_j^{-1} \big )} \le \dfrac{\sigma }{\sqrt{1-\dfrac{(1+\sqrt{2r})^2}{2}}}\sqrt{\mathrm {Tr} \left(\mathbf {S}^{-1} \right)} \\& = &\Big (1-\dfrac{(1+\sqrt{2r})^2}{2}\Big )^{-1/2} {\rho _d(\Phi )}.\) 
This proves the first bound.
For the entropy estimation measure, it follows that
\(\rho _e(\Phi _j) & = & \det (\log (\mathbf {S}_j^{-1})) \\& \le & \log \left(\det \left(({{1-\dfrac{(1+\sqrt{2r})^2}{2}}})^{-1} \mathbf {S}^{-1}\right) \right) \\&= & - n\log \left( 1-\dfrac{(1+\sqrt{2r})^2}{2} \right) +\rho _e(\Phi ),\) 
which proves the second bound.
Proof of Proposition REF :
Taking trace from both sides of (REF ), we get
\(\rho _d^2(\Phi _{i+1})=\rho _d^2(\Phi _{i})+\dfrac{\sigma ^2}{1-\phi ^T\mathbf {S}_i^{-1}\phi } \mathrm {Tr}\left( \mathbf {S}_i^{-1} \phi \phi ^T \mathbf {S}_i^{-1}\right).\) 
Update rule (REF ) follows by utilizing the following equation
\(\mathrm {Tr}\left( \mathbf {S}_i^{-1} \phi \phi ^T \mathbf {S}_i^{-1}\right)=\left\Vert  \mathbf {S}_i^{-1} \phi \right\Vert ^2.\) 
The update rule () for the entropy estimation measure follows from
\(\rho _e(\Phi _{i+1}) &=& \log \left(\mathrm {det}\left(\mathbf {S}_{i+1}^{-1}\right)\right) \\& = &\log \left(\mathrm {det}\left(\mathbf {S}_{i}^{-1}+\dfrac{\mathbf {S}_i^{-1} \phi \phi ^T \mathbf {S}_i^{-1}}{1-\phi ^T\mathbf {S}_i^{-1}\phi }\right)\right) \\& =& \log \left(\left(1+\dfrac{(\mathbf {S}_i^{-1} \phi )^T (\mathbf {S}_i^{-1})^{-1} \mathbf {S}_i^{-1}\phi }{1-\phi ^T\mathbf {S}_i^{-1}\phi }\right)\mathrm {det}(\mathbf {S}_{i}^{-1})\right) \\& =& \log \left(\left(1+\dfrac{\phi ^T\mathbf {S}_i^{-1}\phi }{1-\phi ^T\mathbf {S}_i^{-1}\phi }\right)\mathrm {det}(\mathbf {S}_{i}^{-1})\right) \\& = & \log \left(\left(\dfrac{1}{1-\phi ^T\mathbf {S}_i^{-1}\phi }\right)\mathrm {det}(\mathbf {S}_{i}^{-1})\right) \\& = & \log (\mathrm {det}(\mathbf {S}_{i}^{-1}))+\log \left(\dfrac{ 1}{1-\phi ^T\mathbf {S}_i^{-1}\phi }\right) \\& = & \rho _e(\Phi _i)-\log (1-\phi ^T\mathbf {S}_i^{-1}\phi ).\) 
In the third line, the matrix determinant lemma is applied [21]}.
Proof of Theorem REF :
First, we prove a more general inequality. Let us consider the class of all estimation measures that have the following spectral representation
\(\rho (\Phi )=\sum _{i=1}^n \psi (\lambda _i(\mathbf {S})),\) 

for some convex and monotonically decreasing function \(\psi : \mathbb {R}_+ \rightarrow \mathbb {R}\) . Since \(\psi \)  is convex, we apply Jensen's inequality [22]} and write
\(\dfrac{1}{n}\rho (\Phi ) & = & \dfrac{1}{n}\sum _{i=1}^n \psi (\lambda _i(\mathbf {S})) \nonumber \\& \ge & \psi \left( \sum _{i=1}^n \dfrac{\lambda _i(\mathbf {S})}{n} \right)=\psi \left( \dfrac{\mathrm {Tr}(\mathbf {S})}{n} \right).\) 
On the other hand, we have
\(\mathrm {Tr}\left(\mathbf {S}\right) &=& \mathrm {Tr}\Bigg (~\sum _{i=1}^{|\mathfrak {S}|} \phi _i \phi _i^T ~\Bigg )=\sum _{i=1}^{|\mathfrak {S}|} \hspace{1.42271pt}\mathrm {Tr}\left(\phi _i \phi _i^T \right) \\& = & \sum _{i=1}^{|\mathfrak {S}|} \phi _i^T \phi _i= \sum _{i=1}^{|\mathfrak {S}|} \Vert \phi _i\Vert ^2. \) 
Each vector \(\phi _i\)  corresponds to some time stamp \(t\)  and some index \(j\in \Omega \) , i.e.,
\(\phi _i=\mathrm {e}^{\mathbf {A}^Tt} e_j.\) 
Thus, we use the bound on the matrix exponential to get
\(\Vert \phi _i\Vert  \hspace{1.42271pt}= \hspace{1.42271pt}\Vert \mathrm {e}^{\mathbf {A}^Tt} \hspace{1.42271pt}e_j\Vert  \hspace{1.42271pt}\le \hspace{1.42271pt}\Vert \mathrm {e}^{\mathbf {A}^Tt}\Vert  \hspace{1.42271pt}\Vert e_j\Vert  \hspace{1.42271pt}\le \hspace{1.42271pt}\nu (\mathbf {A},\mathfrak {S}).\) 
This inequality together with (REF ) gives us
\(\mathrm {Tr}\left(\mathbf {S}\right) \hspace{1.42271pt}\le \hspace{1.42271pt}\sum _{i=1}^{|\mathfrak {S}|} \hspace{1.42271pt}\Vert \phi _i\Vert ^2 \hspace{1.42271pt}\le \hspace{1.42271pt}\sum _{i=1}^{|\mathfrak {S}|} \hspace{1.42271pt}\nu (\mathbf {A},\mathfrak {S}) ^2=\nu (\mathbf {A},\mathfrak {S})^2 \hspace{1.42271pt}|\mathfrak {S}|.\) 
Since \(\psi \)  is monotonically decreasing, from (REF ), one may conclude that
\(\dfrac{1}{n}\rho (\Phi ) \hspace{1.42271pt}\ge \hspace{1.42271pt}\psi \left( \dfrac{\mathrm {Tr}(\mathbf {S})}{n} \right) \hspace{1.42271pt}\ge \hspace{1.42271pt}\psi \left( \dfrac{\nu (\mathbf {A},\mathfrak {S})^2 \hspace{1.42271pt}|\mathfrak {S}|}{n} \right).\) 
By applying inequality (REF ) to spectral functions \(\psi (\lambda )=\lambda ^{-1}\)  and \(\psi (\lambda )=-\log (\lambda )\) , we will get the desired inequalities (REF ) and (REF ), respectively.
Proof of Theorem REF :
Every sampling strategy \(\mathfrak {S}\)  that satisfies our assumptions also satisfies
\(\mathfrak {S} \subset \Omega \times \delta \mathbb {Z}.\) 
Therefore, the frame matrix corresponding to such sampling strategy satisfies
\(\mathbf {S}& \preceq &\sum _{i \in \Omega } \hspace{1.42271pt}\hspace{1.42271pt}\sum _{t \in \delta \mathbb {Z}_+} \mathrm {e}^{\mathbf {A}^T t} \hspace{1.42271pt}[\mathbf {C}_{\Omega }]_i^T[\mathbf {C}_{\Omega }]_i \hspace{1.42271pt}\mathrm {e}^{\mathbf {A}t} \\& = & \sum _{k=0}^\infty \Big ( \mathrm {e}^{\mathbf {A}^T \delta } \Big ) ^k \mathbf {C}_{\Omega }^T \mathbf {C}_{\Omega } \Big ( \mathrm {e}^{\mathbf {A}\delta } \Big ) ^k=\mathbf {Q}.\) 
The last equality holds for the following reason. Since \(\mathbf {A}\)  is Hurwitz, \(\mathrm {e}^{\mathbf {A}\delta }\)  is Schur and \((\mathrm {e}^{\mathbf {A}\delta },\mathbf {C}_\Omega )\)  is observable according to (REF ). Inequality (REF ) is equivalent to
\(\mathbf {Q}^{-1} \preceq \mathbf {S}^{-1}.\) 

Functions \(\sqrt{\mathrm {Tr}(\bf X)}\)  and \(\det (\log (\bf X))\)  are nondecreasing on the cone of positive semi-definite matrices. By apply these functions to both sides of (REF ), one obtains the desired inequalities (REF ) and ().
