1. Introduction

Data acquisition is a costly and time consuming operation in any form of tomography. During the data acquisition phase, a larger number of measurements are made to collect enough data to achieve the target precision in reconstructed images. Data collection can reduce the speed of the scanner, or the temporal resolution of the platform, to the point that the system cannot catch up with the speed of changes in the imaging of dynamically evolving objects. Increasing the number of measurements can also increase the operation cost and in some tomography platforms (e.g., X-ray based imaging systems) cause damage to the medium under test. Therefore, it is desirable to reduce the number of measurements as long as the target precision is achieved.

The information content of different measurements are not equal. Some measurements are more informative than others. If, instead of conducting all possible measurements, we compute and only conduct a handful of the most informative measurements, we can speed up the data acquisition process without noticeably sacrificing the precision. The main objective of this article is to introduce a systematic procedure to compute the optimal measurement in each step of the data acquisition process in real-time using the prior information and the outcome of all previous measurements when scanning the sample.

It is possible to improve the information content of measurements by increasing the number or the positions of sources/detectors installed on the scanner or the geometry of the scanner in general. However, without loss of generality, in this article, the focus is on a scanner with a given geometry and a constant number of sources/detectors. Here, we optimize measurements by computing the radiation intensity of each source, a.k.a. the illumination pattern, to make the measurements more informative. Same algorithms can be used to optimize source or detector positioning or other changes in the scanner’s geometry.

In previous studies, two general approaches were tested to find optimal illumination patterns. The first approach was based on the idea that stimulating a system with complex input patterns can provide more information about the system than measuring a set of impulse responses similar to raster scanning. Several illumination patterns were tested, including checkerboard [1], sinusoidal [2], wavelet [3], or diffractive optics patterns [4]. The improvement achieved by these algorithms were not significant. The second approach was designed based on the theory of optimal experiments mostly formulated via the maximum likelihood estimation theory. Here, illumination patterns were designed to minimize the trace of the covariance matrix (A-optimal) as a measure of uncertainty in reconstruction, minimize the determinant (D-optimal) or the norm (E-optimal) of the error covariance matrix [5]. An example of this approach is the study in which they computed illumination patterns to optimize the condition of the Fisher information matrix, which is an example of E-optimal solution, [6]. While this second approach was more successful in improving the quality of reconstructed images, it still has some major drawbacks. In this approach, all illumination patterns are computed in advance at the beginning of scanning mainly based on the geometry of the scanner. Therefore, the prior information about the medium or the outcome of ongoing measurements are not taken into consideration. This is a major drawback specifically when the scanner is imaging dynamically evolving objects. In the current manuscript, the optimization is formulated for a dynamical system and optimal experiments are computed at each iteration based on all available information collect to that point. The theory is formulated based on optimizing experiments in a Kalman estimation filter and it is under the category of A-optimal experiment design.

Kalman filter was introduced in 1960 [7] as a mathematical framework to compute the minimum variance estimation of the state by combining the prior estimation with the outcome of a noise contaminated measurements [8]. Sequential measurements followed by state estimation updates form the main structure of Kalman filter iterations. The theory has been quite successful from early days and has found a large number of applications in different industries [9]. Several more advanced versions of the theory was developed including the Extended Kalman Filter (EKF) algorithm which was designed to make the filter applicable to nonlinear systems [10]. In his original work, Kalman did not offer a recipe for the design of measurements to further improve the precision of computed estimations. In this article, a systematic procedure is introduced to design optimal measurements for Kalman iterations in general. Through this analysis, we also demonstrate the connection between the Kalman estimation filter and the theory of principal components (PCs). This approach can also be used for system identification [11] or sensor position optimization applications in tomography [12]. In the article, the developed theory was applied to a typical optical tomography problem to evaluate the improvement we can be achieved via this optimization process.

2. Theoretical background

In the published literature, Kalman filter, together with the well-known Bayesian learning and the theory of maximum likelihood, cast the main framework of the estimation theory [8]. In the current article, optimal measurements are computed based on the formulation of Kalman filter. Nonetheless, since these theories are closely related, it is possible to formulate the same results, for example based on Bayesian learning. In the following section, we first review the main formulation of the Kalman filter. Then, we show the connection between Kalman filter, Bayesian learning, and the theory of maximum likelihood.

2.1 Basics of the Kalman filter

Consider a linear dynamical system formulated by the following transition model:

In these equations, $n$ is the step in time, $x(n) \in \mathbb {R}^p$ is the state variable, $y(n) \in \mathbb {R}^q$ is the measurement outcome or observation, and $u(n) \in \mathbb {R}^r$ is the input applied to the system. Here, $\varepsilon _p(n) \sim \mathcal {N}(0,R_p)$ and $\varepsilon _m(n) \sim \mathcal {N}(0,R_m)$ are process and measurement noises modeled by a zero-mean normal distribution with covariance matrices $R_p$ and $R_m$, respectively. The process noise $\varepsilon _p(n)$ models our lack of confidence in the model we have for the dynamics of the system and the way input affects the dynamics. Based on Kalman theory, if $x^{n|n-1}$ and $x^{n|n}$ are estimations of the state variable before (prior) and after (posterior) the $n$th measurement, then [8]:

As the equation shows, the posterior estimation is a combination of the prior estimation and the difference between the actual and estimated value of the measurement. The matrix $K(n)$, a.k.a. Kalman gain, is the weight that decides whether we should relay more on the prior estimation or the measurement outcome to obtain a more accurate posterior. Kalman computed the optimal value for $K(n)$ by minimizing the trace of the posterior covariance via solving the equation $\nabla _{K(n)} Tr[P^{n|n}]=0$ to get:

Here, $P^{n|n-1}=Cov(x^{n|n-1})$ and $P^{n|n}=Cov(x^{n|n})$ are the prior and posterior covariance matrices. In Kalman filter, posterior covariance is computed iteratively using the following equation:

After updating the state estimation by Eq. (2), the dynamics is pushed one step forward to compute the prior for the next iteration:

Kalman iterations can continue with a finite or infinite horizon.

2.2 Correlation between Kalman filter, Bayesian learning, and maximum likelihood estimation

We assumed that the prior estimation has the following normal distribution:

By using this equation and the equation of measurement in (1), we can compute the probability distribution of the measurement:

From these two distributions, one can compute the covariance between the state and the corresponding measurement:

In this case, the joint probability distribution of the state and the measurement, $p(x,y)$, prior to the $n$th measurement, has the following normal distribution:

In Bayesian learning, the posterior is defined as the conditional probability $p(x|y)$ which is the probability of the state given the outcome of the measurement. We can prove (Appendix A) that the posterior estimation has the following normal distribution:

By comparing the equation for $x^{n|n}$ in (2) with the same equation in (11), we can see that if $K(n) = \Sigma _{12} ~\Sigma _{22}^{-1}$ the equation for the posterior estimation (and its covariance) in the Bayesian framework matches the Kalman estimation filter.

In Kalman formulation, the posterior covariance is usually calculated from the prior covariance via Eq. (4). It is possible to show that these two covariance matrices are also related by (Appendix B):

In the case of no prior information and before conducting the first measurement, we have infinite amount of uncertainty regarding the state of the system in the next time step or: $P^{1|0} \rightarrow \infty$, which translates to:

We can rewrite the Kalman gain Eq. (3) as:

Simply rearrange this equation to get:

If we use the equation of posterior covariance for the case of no prior, the Kalman gain becomes:

The posterior estimation equation in this case is:

Assume that the coordination is adjusted so that $x^{1|0}=0$, then:

It is possible to show that Kalman state estimation equation for the case of no prior information matches the estimation obtained by the theory of maximum likelihood. To see this, assume that a measurement is made following the Eq. (1), or $y=Cx + \varepsilon _m$. The probability distribution of the observed outcome is:

When the state variable is given, we have: $Cov(x) = 0$. Then, the likelihood probability is:

To find the state that maximizes the probability of this observation, we take the derivative of the logarithm of the likelihood probability w.r.t. the state variable and put it equal to zero:

From this, we compute the maximum likelihood estimation of the state variable, $\hat {x}_{ML}$:

We can compute the covariance of this estimation directly from Eq. (22):

Simplify this equation, to arrive at:

It is now clear that the state estimation and its covariance computed by the maximum likelihood theory, Eqs. (22) and (24), matches the ones computed by the Kalman filter in the case of no prior as stated in Eqs. (18) and (13).

3. Optimal experiment

In a dynamical system formulated by Eq. (1), $C(n)$ is the measurement matrix. Assume that every time we make a measurement, we have the option to choose our measurement matrix. By choosing the best possible measurement, we can further reduce the trace of the posterior covariance, which is our measure of uncertainty. Trace of the posterior is a function of Kalman gain and the measurement matrix. We can consider an iterative process in which we first optimize the trace over the Kalman gain, for a given measurement matrix, then we assume that the gain is constant and we optimize the posterior trace over the measurement matrix and we continue these iterations till convergence. Previously, we took the derivative of the trace of the posterior covariance matrix, w.r.t. Kalman gain, and we found the optimal gain for a given measurement matrix. Now, we take the derivative of the posterior trace w.r.t. the measurement matrix and solve $\nabla _{C(n)} Tr[P^{n|n}]=0$ to find the optimal measurement for the given gain. In fact, this leads to an interesting result. For $C^*(n)$ being the optimal measurement matrix, we get:

Here, $K^{\dagger }(n)=[K^T(n)~K(n)]^{-1}K(n)^T$. This equation suggests that, for a given gain, the optimal measurement matrix is the left pseudo-inverse of $K(n)$. We can compute the optimal gain and measurement matrix using Eqs. (3) and (25) iteratively and continue the iterations until the algorithm converges to an extremum of the trace function. Obviously, the algorithm converges to an extremum point of the trace function (possibly a saddle point) but not necessarily the global minimum, which is the ideal outcome. Nevertheless, Eq. (25) provides a new vision toward computing the optimal experiment.

One question here is: is it possible to find the global optimal measurement matrix without even running the iterations between $K(n)$ and $C(n)$? We focus on answering this question in the following.

It is known that the prior covariance is a symmetric positive definite matrix. Therefore, we can write the following eigenvalue decomposition for the prior covariance:

Here, $W=[w_1~\cdots ~w_p]$ is the matrix of eigenvectors of $P^{n|n-1}$ and $\Lambda =diag[\lambda _1~\cdots ~\lambda _p]$ the matrix of eigenvalues, assuming $\lambda _1 \geq \lambda _2\cdots \geq \lambda _p$. Eigenvectors of the covariance matrix are traditionally known as the principal components (PCs). Eigenvectors with largest eigenvalues are the main PCs of a covariance matrix. For the posterior covariance matrix, when the experiment is optimal, we have:

Assume that $Col[K(n)]$ and $Col^{\bot }[K(n)]$ are the $K(n)$’s column-space and the corresponding orthogonal space. The matrix $I - K(n) K^{\dagger }(n)$ in this equation is a projection operator that projects columns of $P^{n|n-1}$ to $Col^{\bot }[K(n)]$ [13]. This means that the posterior covariance will have no eigenvector in the $Col[K(n)]$.

Trace of the covariance matrix has extremums along it’s principal components. Iterations between $K(n)$ and $C(n)$, depending on the initial point of iterations, converge to one of these extremums. Since the measurement matrix has $q$ rows, $K(n)$ has $q$ columns. We can choose each column of $K(n)$ to be an eigenvector of the prior covariance. In this case, when $C^*(n)=K^{\dagger }(n)$, the drop in the trace of the posterior covariance is equal to the summation of all $q$ corresponding eigenvalues. However, if we choose $K(n)=[w_1~\cdots ~w_q]$ so that its columns are the first main principal components of the prior covariance with the largest eigenvalues, then we can expect to observe the maximum possible drop in the trace of the posterior covariance in an experiment with $q$ detectors. Therefore, $C^*(n)=[w_1~\cdots ~w_q]^{\dagger }$ is the equation for the global optimal experiment.

This analysis shows the intersection between the Kalman filter and the theory of principal components. The analysis suggests that best directions to measure components of the state vector are directions along which we have maximum variances which are the main principal components of the prior covariance.

Sequential iterations between computing optimal $K(n)$ and $C(n)$ leads to the development of a measurement matrix in which each row of $C(n)$ is a principal component of the prior covariance. However, the convergence in these iterations is only asymptotic. To show this, consider the following equation derived from Eq. (15):

In the case of optimal measurement, there is a paradox embedded in this equality which can be resolved only at the limit when the number of iterations goes to infinity. Since for optimal measurement $P^{n|n}$ belongs to the $Col^{\bot }[K(n)]$, it has no component in $Col[K(n)]$. Therefore, the right side of this equation can never match the left side, which is in $Col[K(n)]$, unless they are both zeros. Obviously, this can happen in the absence of noise when $R_m=0$, which means the noise power is zero. However, when $R_m \neq 0$, iterations resolve this paradox by continuously increasing the second norm of rows of the measurement matrix in each $K(n)$-$C(n)$ iteration, which translates to reduction of the second norm of columns of $K(n)$, since $C^*(n)K(n)=I$. In other words, the optimization process continuously increases the power of signal in measurements to minimize the effect of noise. When the number of iterations goes to infinity, the signal power is increased so much that noise is effectively canceled and both sides of the equation match, but only asymptotically.

Now, we answer the question we asked before. Based on this analysis, it is possible to avoid $K(n)$-$C(n)$ iterations. We have two factors when synthesizing the optimal measurement matrix. First, we should choose the right directions (direction of rows of the measurement matrix), which are directions of the main principal components, and then we should increase the power of signal (second norm of rows of the measurement matrix) to the maximum level that is practically allowed.

Convergence of $K(n)$-$C(n)$ iterations for a simple two-dimensional state vector is shown in Fig. 1. Only one measurement is made at a time in this example, therefore $C(n)$ has only one row. The figure shows the result of the first $10$ iterations when iterations were started at different initial points. The graph shows how these iterations gradually align the measurement vector with the main principal component of the prior covariance and then increase the length of the vector to improve the signal power.

 

Fig. 1. Asymptotic convergence of $K(n)$-$C(n)$ iterations. A two-dimensional normal distribution of a state vector is shown in both panels. In this example only one measurement is made at a time. The figure shows the result of the first $10$ iterations for the measurement vector $C(n)$ while starting at two different initial points in the left and right panels. Even when the initial point is closer to the second principal component, the iterations move the measurement vector toward the main PC along which the variance is maximum. As shown, during iterations the second norm of the measurement vector is increasing. In other words, these iterations correct the direction of the measurements and increase the signal power to overcome the noise. If the second norm of the $C(n)$, representing the signal power, is limited to $1.0$, then the optimal measurement is a point along the direction of the main PC on the unit circle shown in both panels.

Download Full Size | PDF

While this approach seems logical, in most real-world applications, there are practical constrains that make it impossible to synthesize the optimal measurement matrix. In such cases, running $K(n)$-$C(n)$ iterations is a reasonable approach to solve the optimization problem and search for an optimal measurement matrix which is practically feasible. In the following, we demonstrate this scenario in an optical tomography problem.

4. Optimal scanning in diffuse optical tomography

Consider the tomography scanner displayed in Fig. 2. Here, the turbid medium under test, modeled by absorption ($\mu _a$) and reduced scattering ($\mu '_s$) coefficients, is located between an array of point light sources and photodetectors [14]. Each source radiates light into the medium which is then detected by the detectors after interacting with the medium. Assume that the scanner has $s$ number of sources, $d$ detectors, and the medium is divided to $v$ pixels. The goal of tomography in this example is to identify the distribution of the absorption coefficient in each pixel, $\eta _k$ for $k=1,\ldots,v$, from the data collected by the detectors. The state vector in this example is $\eta =[\eta _1,\ldots,\eta _v]^T$. The energy radiated from the point source $S_i$, located at $r_{S_i}$, that reaches the pixel $V_k$ at $r_{V_k}$ is $I_i~G(r_{V_k};r_{S_i})$. Here, $I_i$ is the radiation intensity of the source $S_i$. The light that is scattered from this pixel propagates to detectors such that the light intensity sensed by the detector $D_j$, located at $r_{D_j}$, is: $\sum ^V_{k=1}I_i\eta _k G(r_{D_j};r_{V_k}) G(r_{V_k};r_{S_i})$, where $G(r_{V_k};r_{S_i})$ and $G(r_{D_j};r_{V_k})$ are the Green’s functions that model the propagation of the light from the $i$th point source to the $k$th pixel, and the propagation of the scattered light from the $k$th pixel to the $j$th detector. Using the matrix-vector notation, the equation of measurement for this scanner can be written as [11], [6]:

 

Fig. 2. Structure of a diffuse optical tomography scanner. Medium is surrounded by an array of sources and detectors. $G(r_{V_k};r_{S_i})$ and $G(r_{D_j};r_{V_k})$ are Green’s functions that model the propagation of the light from the $i$th source to $k$th pixel, and propagation of the scattered light from the $k$th pixel to the $j$th detector. $\eta _k$ is the concentration of absorbing molecules in pixel $V_k$. The scanner has $s$ number of point sources, $d$ point detectors, and the medium is divided to $v$ number of differential pixels.

Download Full Size | PDF

Here, $G_2$ and $G_1$ are matrices of dimensions $d \times v$ and $v \times s$ with elements $G^{jk}_2=G(r_{D_j};r_{V_k})$ and $G^{ki}_1=G(r_{V_k};r_{S_i})$. In the equation $diag(.)$ is the operator that takes a vector and returns a diagonal matrix. $I=[I_1,I_2,\ldots,I_s]^T$ is the source intensity vector or the illumination pattern. We assume that $0\le I_i \le I_{Max}$ where $I_{Max}$ is the maximum radiation intensity of each source. The measurement matrix in this example is $C(n)=G_2~diag[G_1~I(n)]$, which is determined by the geometry of the scanner and the illumination pattern. We can optimize the performance of the scanner by changing the geometry or illumination patterns. In this example, we assume that the geometry is given and $I(n)$ is the optimization variable. $y = [y_1,y_2,\ldots,y_d]^T$ is the output of the detectors and $\varepsilon _m$ is the measurement noise. Fixed geometry of the scanner and the upper limit for source intensity restrict our ability to synthesize desired measurement matrices.

Traditionally, for scanning a sample, sources are turned on sequentially such that only one source in turned on during each measurement with its maximum radiation intensity and the scanner records the output of all detectors. This protocol is known as raster scanning [14]. It is also possible to turn on all sources with random radiation intensities within the acceptable range. The random scanning protocol is used in some tomography platforms. Our goal is to compare the performance of these two commonly used scanning protocols with the optimal measurement algorithm where the illumination pattern, $I(n)$, is tuned in each measurement to maximize the information we obtain about the state variable $\eta$. By taking the trace of the covariance of both sides of the Eq. (2), we arrive at:

To find the optimal illumination pattern, we minimize $Tr(P^{n|n})$ over the variable $I(n)$. This optimization problem can be formulated as:

For a single row $C(n)$, this optimization is a quadratic programming which is convex [5]. But in general, when the measurement matrix $C(n)$ has more than one row, this optimization is non-convex. Fortunately, it is possible to compute the gradient of the cost function and use gradient-based optimization algorithms (such as the projected gradient descent method [15]) to search for local minimums. We can prove that (Appendix D):

Here, $diag^*(A)$ is the operator that takes the matrix $A$ and returns a vector of diagonal elements of $A$. Consider the two scanner configurations shown in Fig. 3. Scanner dimensions are in an arbitrary unit of length. In the scanner shown at the top, the medium is surrounded by $54$ sources and $23$ detectors and divided to $225$ pixels, $15$ pixels in each direction. The figure also shows the initial actual (left panel) and assumed (right panel) distributions of the absorbing molecules in pixels. The assumed distribution is our prior estimation of the distribution, which is clearly far from the reality in this example. The prior covariance matrix used had eigenvalues uniformly distributed between $0.001$ and $1$ with the condition number of $10^3$. We assumed that point sources and detectors are implemented as optical fiber tips that are assembled on the sides of the scanner. The scanner at the bottom of Fig. 3 is similar to the scanner shown on top, but the number of sources and detectors are reduced to $36$ and $13$, respectively.

 

Fig. 3. Configuration of the tomography scanner used in the simulation: For the top scanner, the turbid medium is surrounded by $54$ sources and $23$ detectors. For the scanner shown at the bottom, the number of sources and detectors are reduced to $36$ and $13$, respectively. The medium is divided to $225$ pixels, $15$ pixel in each direction. The actual and assumed initial distributions are shown in the left and right panels, for both scanners. Dimensions are in an arbitrary unit of length. The figure shows the actual and assumed distributions of the perturbation in the absorption coefficient, and not the background. The color bar shows the percentage of change in the absorption coefficient compared to the background.

Download Full Size | PDF

The following equation was used for the Green’s function:

Figure 4 shows the snapshots of the reconstructed images produced by the first scanner. It was assumed that the system has a dynamics and the distribution is changing over time by moving toward the right-top of the scanner and simultaneously diffusing into the medium (matrix of the dynamics is given in Appendix C). Snapshots of this dynamics are shown in the first row. We assumed that measurements are made much faster than changes in the distribution so that during the period of each measurement the change in the state variable is negligible. Optimal, random, and raster protocols were tested and results are shown in the figure. Borders of the actual and reconstructed distributions are marked by red and green contours, respectively. The simulation clearly demonstrates the superior performance of the optimal scanner regarding the precision of reconstructed images. Figure 5 shows snapshots of the reconstructed images generated by the second scanner with the reduced number of sources and detectors. From the comparison between Figs. 4 and 5 we can see the effect of number of sources/detectors in the quality of reconstructed images. Obviously, larger number of sources/detectors, specifically for small source/detector counts, improve the precision of reconstructed images. However, after some point, this improvement is almost negligible, because of the close correlation between the recordings when these numbers are large.

 

Fig. 4. Snapshots of the actual dynamics (top row), estimated distributions after optimal (second row), random (third row), and raster measurements (forth row) are shown separately in each row generated by the fist scanner. The perturbation in this example is changing as a function of time by slowly diffusing in the medium as moving toward the up right corner of the scanner. Borders of the actual and reconstructed distributions are shown in all snapshots. Snapshot are for state estimations at time steps $n = [0,72,143,215,286,358,429,500]$. The simulations show how optimal measurements outperform the traditional raster and random scanning in terms of revealing more accurate pictures of absorption coefficient distribution in the medium.

Download Full Size | PDF

 

Fig. 5. The same snapshots as Fig. 4 but generated by the second scanner with a smaller number of sources/detectors.

Download Full Size | PDF

Figure 6 shows changes in the normalized uncertainty and error recorded as a function of the measurement number for both scanners. Normalized uncertainty is the trace of the covariance matrix after each measurement divided by the same value before conducting any experiments. Normalized error is the sum of square of differences between the actual and reconstructed pixels divided by the same value before conducting any measurements. Superior performance of the optimal scanning protocol is clear in these plots. For instance, for the first scanner, the normalized error reaches $0.58$ after $188$ measurements. The random scanning protocol reaches the same error after $500$ iterations while raster scanning does not reach this level within the first $500$ measurements. A more interesting result comes from the caparison between the performance of the first and second scanners. The second scanner, with much smaller number of sources and detectors, when following the optimal scanning protocol, has superior performance compared to the first scanner when using either raster or random protocols. For instance, the second scanner reached the normalized error of $0.63$ via optimal scanning at iteration $258$. The first scanner reached this normalized error after about $400$ iterations using the random protocol and about $500$ iterations using the raster protocol. Therefore, following the optimal protocol make it possible to reduce the complexity of the hardware and the operation cost while improving the speed and temporal resolution of the scanner.

 

Fig. 6. Normalized uncertainty and normalized error plotted as a function of the measurement number for optimal, random, and raster scanning algorithms for the first (top) and second (bottom) scanners. Optimal protocol outperforms other algorithms based on these two quantitative measures. Data shows that the optimal algorithm reaches the normalized error $0.58$ after $188$ measurements. The random protocol, which performs better than the raster scan, reaches the same error rate after $500$ measurements.

Download Full Size | PDF

5. Conclusion

A mathematical approach for the design of optimal experiments in Kalman filter was introduced in the article. By optimizing experiments, it is possible to reduce the number of resources used for the measurement to reduce the cost of the operation while achieving the same precision level in the computed estimations. The compromise is in the computational power used to calculate the optimal measurement matrix in each step.

When parameters of the dynamics (matrices $A(n)$ or $B(n)$, the initial value of the state vector and its covariance) are known, we can compute all optimal measurements in advance. Nonetheless, in many applications, parameters of the dynamics are unknown and system identification algorithms, such as the expectation maximization [16], subspace analysis [17], or dynamic mode decomposition [18], [19], are run concurrently to simultaneously estimate these parameters from the collected data in addition to estimating the state variable. In those cases, we cannot compute the direction of optimal measurements in advance and optimal measurements should be computed online. Conducting optimal experiments can also significantly improve the performance of system identification algorithms and the precision of state estimations [11].

There are two aspects to further investigate in the future. The first aspect is around the concept of observability. In control theory, observability is a measure of how well the state of a system can be inferred from feasible partial measurement [13]. Observability is a function of the measurement matrix. By optimizing the measurement matrix, we can improve observability of the partial measurements in closed-loop control systems. The approach taken here can be applied to the optimization of the observability Gramian and synthesize experiments that offer optimal observability.

The second direction for future studies is in the field of signal sampling. It was shown in the theory of compressed sensing that a signal which has a sparse representation in an isometric transform domain can be reconstructed from far less number of samples compared to what is prescribed by the Nyquist theorem [20], [21]. The theory of optimal measurements in this article was formulated around the idea of using the covariance between different components of the hidden state vector to find the most informative experiments. In the signal sampling problem, it is possible to use the correlation between values of a signal at different points in time and only sample the signal at limited points and yet reconstruct the entire signal, with high confidence, based on these correlations. The sparsity that is used in compressed sensing is the result of some form of correlation between samples that is extracted by an isometric transform. But do we know which transform is optimal and generate the most sparse representation? The optimal experiment theory as discussed can perform the same task but it has the potential to find the optimal transform domain. Therefore, there is a potential to surpass the performance of the traditional compressed sensing algorithm.