Greedy Sensor Selection for Weighted Linear Least Squares Estimation under Correlated Noise

A linear measurement equation for $p$ sensors $\mathbf{y}\in\mathbb{R}^{p}$ and a state vector of $r$ components $\mathbf{z}\in\mathbb{R}^{r}$ is corrupted by Gaussian noise $\mathbf{w}\sim\mathcal{N}(\mathbf{0}|\mathbf{R})\in\mathbb{R}^{p}$, which is independent of $\mathbf{z}$,

We assume that the sensor characteristic $\mathbf{C}\in\mathbb{R}^{p{\times}r}$ is known in advance and nonsingular, and the covariance of the measurement noise $\mathbf{R}\in\mathbb{R}^{p{\times}p}$ is positive definite and symmetric. An parameter vector $\tilde{\mathbf{z}}$ is estimated from Eq. 1:

Note that Eq. (2a) corresponds to the minimal norm solution in which the measurement noise is not considered, though Eq. (2a) is derived from the formulation including measurement noise. On the other hand, Eq. (2b) is a minimum variance unbiased estimation considering measurement noise as the generalized least squares estimation [56, Section 4.5]. The present study focuses on the formulation above, excluding any prior distribution of the state variables.

We also assume a large number of possible measurement points, e.g., $\mathbf{x}\in\mathbb{R}^{n}\,(n\gg r)$. The linear coefficients and noise covariance for all of the measurement points are expressed as $\mathbf{U}\in\mathbb{R}^{n\times r}$ and $\mathcal{R}\in\mathbb{R}^{n\times n}$, respectively. Actual calculations for these terms are introduced later herein at subsection 2.2. With these notations, the measurement is expressed by substituting $\left(\mathbf{y},\,\mathbf{C},\,\mathbf{R}\right)\leftarrow\left(\mathbf{x},\,\mathbf{U},\,\mathcal{R}\right)$. Here, the estimation Eq. 2a is redundant if $r$ is small, and thus the reduced measurement $p\ll n$ is sufficient in terms of both estimation quality and calculation efficiency. A sensor indication matrix $\mathbf{H}{(\mathcal{S}_{p})}\in\mathbb{R}^{p{\times}n}$ is defined for a set $\mathcal{S}_{p}$ of $p$ sensor indices selected from $n$ candidates. The position of unity in the $i$-th row of $\mathbf{H}{(\mathcal{S}_{p})}$ is associated with the $i$-th component of $\mathcal{S}_{p}$, whereas the rest of the row is zero. Measurements and linear coefficients for the selected sensors are denoted as $\mathbf{y}{(\mathcal{S}_{p})}=\mathbf{H}{(\mathcal{S}_{p})}\mathbf{x}$ and $\mathbf{C}{(\mathcal{S}_{p})}=\mathbf{H}{(\mathcal{S}_{p})}\mathbf{U}$, respectively. In addition, the covariance matrix for measurement noise is expressed by $\mathbf{R}({\mathcal{S}_{p}})=\mathbf{H}{(\mathcal{S}_{p})}\mathcal{R}\mathbf{H}{(\mathcal{S}_{p})}^{\top}$. The argument ${(\mathcal{S}_{p})}$ will be denoted as subscript $\circ_{\mathcal{S}_{p}}$ for brevity hereinafter.

In our implementation, the matrices $\mathbf{U}$ and $\mathcal{R}$ are generated by modal decomposition of the collected data matrix, in a process known as principal component analysis, or proper orthogonal decomposition  [7, 33]. The collected data $\mathbf{X}=\left[\mathbf{x}_{1},\,\ldots\,,\mathbf{x}_{m}\right]$ are assumed to consist of $n$-point measurements in rows by $m$ instances in columns, where $n\gg m$. $\mathbf{X}$ is decomposed by singular value decomposition into $m$ orthonormal spatial modes $\mathbf{U}_{X}$ and temporal modes $\mathbf{V}_{X}$, and a diagonal matrix of singular values $\mathbf{\Sigma}_{X}$. The approximation mode number $r$ is chosen to retain the covariance matrix for the original data matrix at a high rate:

where $\mathbf{X},\mathbf{U}_{X}\in\mathbb{R}^{n{\times}m}$ and $\mathbf{\Sigma}_{X},\mathbf{V}_{X}\in\mathbb{R}^{m{\times}m}$, and $\mathbf{U}\in\mathbb{R}^{n{\times}r}$, $\mathbf{\Sigma}\in\mathbb{R}^{r{\times}r}$ and $\mathbf{V}\in\mathbb{R}^{m{\times}r}$, respectively. Here, the second term with the subscript $\circ_{N}$ on the right-hand side is the portion that is excluded from $r$ rank representation and thus is regarded as the measurement noise. The $i$-th column of $\mathbf{\Sigma}\mathbf{V}^{\top}$ and $\mathbf{H}_{\mathcal{S}_{p}}\mathbf{X}$ are the variable vector and the measurement in Eq. 1 for the $i$-th instance, respectively. Thus, one can immediately recover the low rank representation of the large-scale measurement as $\tilde{\mathbf{x}}=\mathbf{U}\tilde{\mathbf{z}}$ by obtaining an estimation $\tilde{\mathbf{z}}$ [33]. With respect to $\mathcal{R}$, several approaches are capable of expressing the statistical behavior of the residual between the measurement and the reduced order model, $\mathbf{X}-\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, which are exemplified by kernel functions used in signal processing or data-driven modeling in Ref.  [64, Section 2]. By taking the expectation $\mathcal{R}=\mathbf{E}\left[\mathbf{w}\mathbf{w}^{\top}\right]$, the model of noise in the latter manner is denoted as $\mathcal{R}=\mathbf{U}_{N}\mathbf{\Sigma}_{N}^{2}\mathbf{V}_{N}^{\top}$ of Eq. 3, which is used in section 3. Appendix A shows how the data-driven design of the measurement noise covariance is affected in the standpoint of the correlation and amount of training data.

Several criteria for sensor selection are available for scalar evaluation of the measurement design, like D-, E- or A-optimality mentioned in [62].The performances for D-, E- and A-optimality criteria were previously compared for sensor sets obtained by greedy sensor selection methods suited for these criteria [38]. Sensor selection based on the D-optimality criterion performed well in both of the computational time and the values of other criteria, thus being adopted in the present study. Note that the efficient implementation shown in subsection 2.4 can easily be extended to the A-optimality settings.

Geometrically, this optimization corresponds to the minimization of the volume of an ellipsoid, which represents the expected estimation error variance [23]:

where the operator $E\left[\circ\right]$ means taking the expected value of the argument. Furthermore, this matrix is known to correspond to the inverse of the Fisher information matrix. This equality is easily confirmed under the assumption of Gaussian measurement noise, by differentiating by $\mathbf{z}$ a log likelihood, $L=-(\mathbf{y}-\mathbf{C}_{\mathcal{S}_{p}}\tilde{\mathbf{z}})\mathbf{R}^{-1}(\mathbf{y}-\mathbf{C}_{\mathcal{S}_{p}}\tilde{\mathbf{z}})+\textrm{const.}$, then substituting estimation given by  Eq.(2b). The optimization returns the set of measurement point $\mathcal{S}_{p}$ from all the possible locations, although this is normally an intractable process. Instead, a greedy algorithm is employed with objective functions for both $p\leq r$ and $p>r$. They are derived by generalizing the formulation in Ref. [38] for the correlated measurement noise hereafter. The set of sensors are evaluated only in the observable subspace of $\mathbf{R}_{\mathcal{S}_{p}}^{-1/2}\mathbf{C}_{\mathcal{S}_{p}}$, since the measurement system Eq. 1 is underdetermined when $p<r$. From Eq. 4, the subspace is separated by the projection $\mathbf{z}\rightarrow\mathbf{\xi}=\hat{\mathbf{V}}_{\mathbf{C}}^{\top}\mathbf{z}$ after singular value decomposition of $\mathbf{R}_{\mathcal{S}_{p}}^{-1/2}\mathbf{C}_{\mathcal{S}_{p}}$:

with some matrices of appropriate dimensions,

The evaluation of the error covariance in the observable subspace was recently introduced by Nakai et al. [38], and it is extended to the correlated noise case in the present study for the first time (to the best of our knowledge). One can use various metrics for the projected covariance matrix Eq. 7 like its determinant, trace or minimum eigenvalue, as in Ref. [23, 38]. The determinant of the inverse matrices in Eq. 7 is maximized in the present manuscript:

Note that whitening all candidates with $\mathcal{R}^{-1}$ before sensor selection is based on the assumption of weakly correlated noise, because $\mathcal{R}$ contains noise covariance over sensors that are not selected as pointed out in [32]. In subsection 2.4, an algorithm is presented for achieving Eq. 17b in a greedy manner.

Algorithm 1 shows the procedure implemented in the computation conducted in section 3, which implicitly exploits the one-rank determinant lemma as [32, 50, 64]. It is worth mentioning that the previously presented noise-ignoring algorithm in [50], Algorithm 2, is easily obtained by substituting an identity matrix into $\mathcal{R}$.

The equations are converted by the lemma shown later herein. First, consider an objective function when there are fewer sensors deployed than the number of state variables $(p\leq r)$:

where the subscript ${k(i)}$ represents the component produced by the $i$-th sensor candidate in the $k$-th step:

Here, $\mathbf{u}_{i}$ is the $i$-th row of $\mathbf{U}$, and $\mathbf{s}_{k(i)}$ and ${t}_{i}$ are the noise covariance between the selected sensors given by the previous steps and the $i$-th candidate and noise variance for the $i$-th candidate, respectively. The algorithm avoids expensive computations involving the determinant by separating the components of the obtained sensors from the objective function in the current selection step of Eq. 18:

and then a unit vector $\mathbf{e}_{i_{k}}\in\mathbb{R}^{1{\times}n}$, of which only the $i_{k}$-th entry is unity, is added to the $k$-th row of $\mathbf{H}$. Note that Eq. 19 corresponds to maximization of the difference when an arbitrary sensor is added to the sensor set of the previous step. The numerator of Eq. 19 is the $\ell_{2}$ norm of the vector, and the denominator is positive, because the covariance matrix $\mathbf{R}_{\mathcal{S}_{k-1}}$ is assumed to be positive definite. Subsequently, the objective function is modified for the case in which more sensors than the number of state variables have already been determined:

where $\mathbf{\phi}_{(i)}=\mathbf{s}_{k(i)}\mathbf{R}_{\mathcal{S}_{k-1}}^{-1}\mathbf{C}_{\mathcal{S}_{k-1}}-\mathbf{u}_{i}.$ Eq. 20 is positive, and the objective function, Eq. (17b), increases monotonically. Details of the expansion are found in Ref. [64].

The computational cost of algorithms are listed in Table 1:

The objective function loses submodularity if the measurement noises at different sensor positions are strongly correlated with each other (or when the off-diagonal components in $\mathbf{R}$ are no smaller than the diagonal components.) The following example provides a nonsubmodular and nonsupermodular case for Eq. (17b). For simplicity, the spatial modes $\mathbf{U}\in\mathbb{R}^{3}$ and the noise covariance $\mathcal{R}\in\mathbb{R}^{3{\times}3}$ are set as follows. Here, the noise components $i=2,3$ are strongly correlated, whereas those for $i=1$ are relatively independent.

With these matrices, the values of the determinant function Eq. (17b) are

where $f_{\mathcal{S}}\,(\mathcal{S}\in{2}^{\{1,2,3\}})$ refers to the value of the determinant Eq. (17b) for the power sets of selected sensors. This example immediately shows that the objective function Eq. (17b) has neither submodularity nor supermodularity, whereas submodularity exists for the case with equally distributed uncorrelated measurement noise [50]. Thus, the sensor selection problem Eq. (17b) with a greedy method generally has no performance guarantee based on submodularity or supermodularity.

Generalized results are shown in this subsection. The problem considered here is as follows: a data matrix $\mathbf{X}$ is constructed as $\mathbf{X}=\mathbf{U}_{X}\Sigma_{X}\mathbf{V}_{X}^{\top}$, where $\mathbf{U}_{X}$ and $\mathbf{V}_{X}$ are $5,000\times 100$ and $100\times 100$ orthonormal matrices, respectively, generated from appropriately sized matrices containing numbers from a standard normal distribution, and $\Sigma_{X}$ is a diagonal matrix with $\operatorname{diag}\left(\Sigma_{X}\right)=\left(\begin{array}[]{llllll}1&0.99&\ldots&(101-j)/100&\ldots&0.01\end{array}\right)$. The algorithms for the sensor selection treat these matrices after dividing the first 10 columns as $\mathbf{U}$ and $\mathbf{V}$ and the remaining columns as $\mathbf{U}_{N}$ and $\mathbf{V}_{N}$. Then, the first 10 diagonal components and the remaining are labeled as $\mathbf{\Sigma}$ and $\mathbf{\Sigma}_{N}$, respectively. Note that the measure $e$ in terms of “reconstruction error” is expressed as $e={\|\mathbf{X}-\mathbf{U}\tilde{\mathbf{Z}}\|_{\mathrm{F}}}/{\|\mathbf{X}\|_{\mathrm{F}}}.$ Here, the series for the estimation $\tilde{\mathbf{z}}$ of Eq. 2a is concatenated as $\tilde{\mathbf{Z}}$, and $\|\circ\|_{\mathrm{F}}$ represents the Frobenius norm of $\circ$. Figure 1 shows the result of the reconstruction with the estimate with $p$ sensors and the $r$ dimensional reduced-order model Eq. 3. Here, DG and DG/NC in the legend refer to “determinant-based greedy algorithm” in Ref. [50] and Algorithm 1 considering “noise covariance” in the measurement, respectively, and LS and GLS refer to “linear least squares estimation” and “generalized linear least squares estimation” using noise covariance, respectively. Note that the plots for $p\leq r$ are calculated by the same estimator Eq. (2a), and, therefore, the estimations with a small number of sensors for both GLS and LS are identical for each selected sensor set. First, the GLS estimation reduces the reconstruction error in oversampling cases for sensors for both algorithms. The measurement noise is quite excessive, and thus, sensors of both selection methods exhibit comparable results for the LS estimation. Second, the more sensors are deployed, the lower the reduction becomes thanks to the GLS estimation. This is partly because measurement using a large number of sensors suppresses outliers resulting from the correlated measurement noise. If a much larger number of sensors is available than the number of estimated variables, the importance of correlation in the measurement noise might diminish.

Here, we apply this strategy to pursue sensor selection using large-dimensional climate data. A brief description of the NOAA-SST data is given in Table 2.

Similar to subsection 3.1, orthonormal modes are prepared by conducting SVD on the data matrix with the average being subtracted, then the first 10 of the 520 modes are used to build the reduced-order model of the temperature distribution ($r=10$). The remaining modes are used for the noise covariance matrix. Several sensor positions for which the noise amplitude is extremely low (smaller than 1% of the maximum RMS of noise in this comparison) are eliminated from candidate set $\mathcal{S}_{n}$ beforehand, as conducted in [64]. The results in this subsection can be compared with those in Ref. [33], [64], or [50]. In Fig. 2, the positions of sensors are represented by open circles on the colored maps, which illustrate the fluctuation of estimation error using those sensors, namely $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}-\mathbf{U}\tilde{\mathbf{z}}$. The difference in the sensor positions is remarkable, since the proposed algorithm spreads sensors and avoids neighboring sensors that might be affected by correlated measurement noise. The reduction in the estimation error is also recognizable by comparing the backgrounds.

Figure 3 compares the results of estimation using the noise covariance information. Note that cross-validation is not conducted for this comparison, since it is hard to quantify the estimation error because of the dynamics in SST which is partly extracted by the reduced order modeling. In appendix A, the covariance matrix of measurement noise is characterized by the number of snapshots to form the noise covariance matrix. A horizontal broken line in Fig. 3 shows the modeling error due to the low-rank representation of Eq. 3. The red plots show better performance for sensors using the proposed algorithm than those using the previous DG algorithm [50] owing to the noise covariance matrix in the sensor selection procedure. There are several differences in the trend of plots compared to Fig. 1, e.g.  the contribution from sensors of the proposed algorithm is more significant than that for the GLS estimation. This is perhaps because of the weak amplitude of higher ordered modes in addition to the similarity in the location where the reduced-order phenomena and the measurement noise fluctuates greatly. The proposed algorithm that involves noise covariance evaluates the positions with less measurement noise. Therefore, accurate estimation is achieved even with linear least squares estimation, which contrasts with the errors of sensors of DG algorithm staying relatively high.