Sparse Signal Recovery with Temporally Correlated Source Vectors Using Sparse Bayesian Learning

I Introduction

Sparse signal recovery, or compressed sensing, is an emerging field in signal processing . The basic mathematical model is

Motivated by many applications such as EEG/MEG source localization and DOA estimation, where a sequence of measurement vectors are available, the basic model (1) has been extended to the multiple measurement vector (MMV) model in , given by

It has been shown that compared to the SMV case, the successful recovery rate can be greatly improved using multiple measurement vectors . For example, Cotter and Rao showed that by taking advantage of the MMV formulation, one can relax the upper bound in the uniqueness condition for the solution. Tang, Eldar and their colleagues showed that under certain mild assumptions the recovery rate increases exponentially with the number of measurement vectors $L$. Jin and Rao analyzed the benefits of increasing $L$ by relating the MMV model to the capacity regions of MIMO communication channels. All these theoretical results reveal the advantages of the MMV model and support increasing $L$ for better recovery performance.

However, under the common sparsity assumption we cannot obtain many measurement vectors in practical applications. The main reason is that the sparsity profile of practical signals is (slowly) time-varying, so the common sparsity assumption is valid for only a small $L$ in the MMV model. For example, in EEG/MEG source localization there is considerable evidence that a given pattern of dipole-source distributions In this application the set of indexes of nonzero rows in $\mathbf{X}$ is called a pattern of dipole-source distribution. may only exist for 10-20 ms. Since the EEG/MEG sampling frequency is generally 250 Hz, a dipole-source pattern may only exist through 5 snapshots (i.e. in the MMV model $L=5$). In DOA estimation , directions of targets In this application the index of a nonzero row in $\mathbf{X}$ indicates a direction. are continuously changing, and thus the source vectors that satisfy the common sparsity assumption are few. Of course, one can increase the measurement vector number at the cost of increasing the source number, but a larger source number can result in degraded recovery performance.

Motivated by applications where signals and other types of data often contain some kind of structures, many algorithms have been proposed , which exploit special structures in the source matrix $\mathbf{X}$. However, most of these works focus on exploiting spatial structures (i.e. the dependency relationship among different sources) and completely ignore temporal structures. Besides, for tractability purposes, almost all the existing MMV algorithms (and theoretical analysis) assume that the sources are independent and identically distributed (i.i.d.) processes. This contradicts the real-world scenarios, since a practical source often has rich temporal structures. For example, the waveform smoothness of biomedical signals has been exploited in signal processing for several decades. Besides, due to high sampling frequency, amplitudes of successive samplings of a source are strongly correlated. Recently, Zdunek and Cichocki proposed the SOB-MFOCUSS algorithm, which exploits the waveform smoothness via a pre-defined smoothness matrix. However, the design of the smoothness matrix is completely subjective and not data-adaptive. In fact, in the task of sparse signal recovery, learning temporal structures of a source is a difficult problem. Generally, such structures are learned via a training dataset (which often contains sufficient data without noise for robust statistical inference) . Although effective for some specific signals, this method is limited. Having noticed that the temporal structures strongly affect the performance of existing algorithms, in we derived the AR-SBL algorithm, which models each source as a first-order autoregressive (AR) process and learns AR coefficients from the data per se. Although the algorithm has superior performance compared to MMV algorithms in the presence of temporal correlation, it is slow, which limits its applications. As such, there is a need for efficient algorithms that can deal more effectively with temporal correlation.

In this work, we present a block sparse Bayesian learning (bSBL) framework, which transforms the MMV model (2) to a SMV model. This framework allows us to easily model the temporal correlation of sources. Based on it, we derive an algorithm, called T-SBL, which is very effective but is slow due to its operation in a higher dimensional parameter space resulting from the MMV-to-SMV transformation. Thus, we make some approximations and derive a fast version, called T-MSBL, which operates in the original parameter space. Similar to T-SBL, T-MSBL is also effective but has much lower computational complexity. Interestingly, when compared to MSBL, the only change of T-MSBL is the replacement of $\|\mathbf{X}_{i\cdot}\|_{2}^{2}$ with the Mahalanobis distance measure, i.e. $\mathbf{X}_{i\cdot}\mathbf{B}^{-1}\mathbf{X}_{i\cdot}^{T}$, where $\mathbf{B}$ is a positive definite matrix estimated from data and can be partially interpreted as a covariance matrix. We analyze the global minimum and the local minima of the two algorithms’ cost function. One of the key results is that in the noiseless case the global minimum is at the sparsest solution. Extensive experiments not only show the superiority of the proposed algorithms, but also provide some interesting (even counter-intuitive) phenomena that may motivate future theoretical study.

The rest of the work is organized as follows. In Section II we present the bSBL framework. In Section III we derive the T-SBL algorithm. Its fast version, the T-MSBL algorithm, is derived in Section IV. Section V provides theoretical analysis on the algorithms. Experimental results are presented in Section VI. Finally, discussions and conclusions are drawn in the last two sections.

We introduce the notations used in this paper:

Bold symbols are reserved for vectors and matrices. Particularly, $\mathbf{I}_{L}$ denotes the identity matrix with size $L\times L$. When the dimension is evident from the context, for simplicity, we just use $\mathbf{I}$;

For a matrix $\mathbf{A}$, $\mathbf{A}_{i\cdot}$ denotes the $i$-th row, $\mathbf{A}_{\cdot i}$ denotes the $i$-th column, and $\mathbf{A}_{i,j}$ denotes the element that lies in the $i$-th row and the $j$-th column;

II Block Sparse Bayesian Learning Framework

To exploit the temporal correlation, we propose another SBL framework, called the block sparse Bayesian learning (bSBL) framework. In this framework, the MMV model is transformed to a block SMV model. In this way, we can easily model the temporal correlation of sources and derive new algorithms.

First, we assume all the sources $\mathbf{X}_{i\cdot}$ ($\forall i$) are mutually independent, and the density of each $\mathbf{X}_{i\cdot}$ is Gaussian, given by

where $\gamma_{i}$ is a nonnegative hyperparameter controlling the row sparsity of $\mathbf{X}$ as in the basic SBL . When $\gamma_{i}=0$, the associated $\mathbf{X}_{i\cdot}$ becomes zeros. $\mathbf{B}_{i}$ is a positive definite matrix that captures the correlation structure of $\mathbf{X}_{i\cdot}$ and needs to be estimated.

Assume elements in the noise vector $\mathbf{v}$ are independent and each has a Gaussian distribution, i.e. $p(v_{i})\sim\mathcal{N}(0,\lambda)$, where $v_{i}$ is the $i$-th element in $\mathbf{v}$ and $\lambda$ is the variance. For the block model (3), the Gaussian likelihood is

Using the Bayes rule we obtain the posterior density of $\mathbf{x}$, which is also Gaussian,

So given all the hyperparameters $\lambda,\gamma_{i},\mathbf{B}_{i},\forall i$, the MAP estimate of $\mathbf{x}$ is given by:

where the last equation follows the matrix identity $(\mathbf{I}+\mathbf{AB})^{-1}\mathbf{A}\equiv\mathbf{A}(\mathbf{I}+\mathbf{BA})^{-1}$, and $\mathbf{\Sigma}_{0}$ is the block diagonal matrix given by (4) with many diagonal block matrices being zeros. Clearly, the block sparsity of $\mathbf{x}^{*}$ is controlled by the $\gamma_{i}$’s in $\mathbf{\Sigma}_{0}$: during the learning procedure, when $\gamma_{k}=0$, the associated $k$-th block in $\mathbf{x}^{*}$ becomes zeros, and the associated dictionary vectors $\boldsymbol{\phi}_{k}\otimes\mathbf{I}_{L}$ are pruned out In practice, we judge whether $\gamma_{k}$ is less than a small threshold, e.g. $10^{-5}$. If it is, then the associated dictionary vectors are pruned out from the learning procedure and the associated block in $\mathbf{x}$ is set to zeros..

To estimate the hyperparameters we can use evidence maximization or Type-II maximum likelihood . This involves marginalizing over the weights $\mathbf{x}$ and then performing maximum likelihood estimation. We refer to the whole framework including the solution (7) and the hyperparameter estimation as the block sparse Bayesian learning (bSBL) framework. Note that in contrast to the original SBL framework, the bSBL framework models the temporal structures of sources in the prior density via the matrices $\mathbf{B}_{i}$ ($i=1,\cdots,M$). Different ways to learn the matrices result in different algorithms. We will discuss the learning of these matrices and other hyperparameters in the following sections.

III Estimation of Hyperparameters

To find the hyperparameters $\Theta=\{\gamma_{1},\cdots,\gamma_{M},\mathbf{B},\lambda\}$, we employ the Expectation-Maximization (EM) method to maximize $p(\mathbf{y};\Theta)$. This is equivalent to minimizing $-\log p(\mathbf{y};\Theta)$, yielding the effective cost function:

where $\mathbf{\Sigma}_{y}\triangleq\lambda\mathbf{I}+\mathbf{D}\mathbf{\Sigma}_{0}\mathbf{D}^{T}$. The EM formulation proceeds by treating $\mathbf{x}$ as hidden variables and then maximizing:

To estimate $\boldsymbol{\gamma}\triangleq[\gamma_{1},\cdots,\gamma_{M}]$ and $\mathbf{B}$, we notice that the first term in (9) is unrelated to $\boldsymbol{\gamma}$ and $\mathbf{B}$. So, the Q function (9) can be simplified to:

It can be shown that The $\propto$ notation is used to indicate that terms that do not contribute to the subsequent optimization of the parameters have been dropped. This convention will be followed through out the paper.

The derivative of (10) with respect to $\gamma_{i}\;(i=1,\cdots,M)$ is given by

where we define (using the MATLAB notations)

So the learning rule for $\gamma_{i}\;(i=1,\cdots,M)$ is given by

On the other hand, the gradient of (10) over $\mathbf{B}$ is given by

Thus we obtain the learning rule for $\mathbf{B}$:

To estimate $\lambda$, the Q function (9) can be simplified to

where (14) follows from the first equation in (6), and $\widehat{\lambda}$ denotes the estimated $\lambda$ in the previous iteration. The $\lambda$ learning rule is obtained by setting the derivative of (15) over $\lambda$ to zero, leading to

where the $\lambda$ on the right-hand side is the $\widehat{\lambda}$ in (15). There are some challenges to estimate $\lambda$ in SMV models. This, however, is alleviated in MMV models when considering temporal correlation. We elaborate on this next.

However, our learning rule (16) does not have such ambiguity problem. To see this, we now examine the covariance matrix $\mathbf{\Sigma}_{y}$ in our cost function (8). Noting that $\mathbf{D}=\mathbf{\Phi}^{\prime}\otimes\mathbf{I}$, we have

Obviously, since $\mathbf{B}$ is not an identity matrix Note that even all the sources are i.i.d. processes, the estimated $\mathbf{B}$ in practice is not an exact identity matrix., $\lambda$ and $\{\gamma_{M-N+1},\cdots,\gamma_{M}\}$ cannot identically contribute to $\mathbf{\Sigma}_{y}$.

The SBL algorithm using the learning rules (6), (7), (12), (13) and (16) is denoted by T-SBL.

IV An Efficient Algorithm Processing in the Original Problem Space

The proposed T-SBL algorithm has excellent performance in terms of recovery performance (see Section VI). But it is not fast because it learns the parameters in a higher dimensional space instead of the original problem space T-SBL can be directly used to solve the block sparsity models . In this case, the algorithm directly performs in the original parameter space and thus it is not slow (compared to the speed of some other algorithms for the block sparsity models).. For example, the dictionary matrix is of the size $NL\times ML$ in the bSBL framework, while it is only of the size $N\times M$ in the original MMV model. Interestingly, the MSBL developed for i.i.d. sources has complexity $\mathcal{O}(N^{2}M)$ and does not exhibit this drawback . Motivated by this, we make a reasonable approximation and back-map T-SBL to the original space By back-mapping, we mean we use some approximation to simplify the algorithm such that the simplified version directly operates in the parameter space of the original MMV model..

For convenience, we first list the MSBL algorithm derived in :

An important observation is the lower dimension of the matrix operations involved in this algorithm. We attempt to achieve similar complexity for the T-SBL algorithm by adopting the following approximation:

which is exact when $\lambda=0$ or $\mathbf{B}=\mathbf{I}_{L}$. For high signal-to-noise ratio (SNR) or low correlation the approximation is quite reasonable. But our experiments will show that our algorithm adopting this approximation performs quite well over a broader range of conditions (see Section VI).

Now we use the approximation to simplify the $\gamma_{i}$ learning rule (12). First, we consider the following term in (12):

where (21) follows the second equation in (6), and $\mathbf{\Xi}_{x}$ is given in (17). Using the same approximation (), the $\boldsymbol{\mu}_{x}$ in (12) can be expressed as

where (23) follows (5) and the approximation (), and $\mathbf{X}$ is given in (18). Therefore, based on (22) and (24), we can transform the $\gamma_{i}$ learning rule (12) to the following form:

To simplify the $\mathbf{B}$ learning rule (13), we note that

where (26) uses the approximation (). Using the definition (11), we have $\mathbf{\Sigma}_{x}^{i}=(\mathbf{\Xi}_{x})_{ii}\mathbf{B}$. Therefore, the learning rule (13) becomes:

From the learning rule above, we can directly construct a fixed-point learning rule, given by

where $\rho=\frac{1}{M}\sum_{i=1}^{M}\gamma_{i}^{-1}(\mathbf{\Xi}_{x})_{ii}$. To increase the robustness, however, we suggest using the rule below:

where $\eta$ is a positive scalar. This regularized form (30) ensures that $\widetilde{\mathbf{B}}$ is positive definite.

Similarly, we simplify the $\lambda$ learning rule (16) as follows:

We denote the algorithm using the learning rules (17), (18), (25), (28), (29) (or (30)), and () by T-MSBL (the name emphasizes the algorithm is a temporal extension of MSBL). Note that T-MSBL cannot be derived by modifying the cost function of MSBL.

Comparing the $\gamma_{i}$ learning rule of T-MSBL (Eq.(25)) with the one of MSBL (Eq.(19)), we observe that the only change is the replacement of $\|\mathbf{X}_{i\cdot}\|_{2}^{2}$ with $\mathbf{X}_{i\cdot}\mathbf{B}^{-1}\mathbf{X}_{i\cdot}^{T}$, which incorporates the temporal correlation of the sources. Hence, T-MSBL has only extra computational load for calculating the matrix $\mathbf{B}$ and the item $\mathbf{X}_{i\cdot}\mathbf{B}^{-1}\mathbf{X}_{i\cdot}^{T}$ Here we do not compare the $\lambda$ learning rules of both algorithms, since in some cases one can feed the algorithms with suitable fixed values of $\lambda$, instead of using the $\lambda$ learning rules. However, the computational load of the simplified $\lambda$ learning rule of T-MSBL is also not high.. Since the matrix $\mathbf{B}$ has a small size and is positive definite and symmetric, the extra computational load is low.

Note that $\mathbf{X}_{i\cdot}\mathbf{B}^{-1}\mathbf{X}_{i\cdot}^{T}$ is the quadratic Mahalanobis distance between $\mathbf{X}_{i\cdot}$ and its mean (a vector of zeros). In the following section we will get more insight into this change.

V Analysis of Global Minimum and Local Minima

Since our bSBL framework generalizes the basic SBL framework, many proofs below are rooted in the theoretic work on the basic SBL . However, some essential modifications are necessary in order to adapt the results to the bSBL model. Due to the equivalence of the original MMV model (2) and the transformed block sparsity model (3), in the following discussions we use (2) or (3) interchangeably and per convenience.

We have the following result on the global minimum of the cost function (8) For convenience, in this theorem we consider the cost function with $\mathbf{\Sigma}_{0}$ given by (4), i.e. the one before we use our strategy to avoid the overfitting.:

Note that $\widehat{\boldsymbol{\gamma}}$ is a function of the estimated $\widehat{\mathbf{B}}_{i}$ ($\forall i$). However, the theorem implies that even when the estimated $\widehat{\mathbf{B}}_{i}$ is different from the true $\mathbf{B}_{i}$, the estimated sources are the true sources at the global minimum of the cost function. As a reminder, in deriving our algorithms, we assumed $\mathbf{B}_{i}=\mathbf{B}$ ($\forall i$) to avoid overfitting. Theorem 1 ensures our algorithms using this strategy also have the global minimum property. Also, the theorem explains why MSBL has the ability to exactly recover true sources in noiseless cases even when sources are temporally correlated. But we hasten to add that this does not mean $\mathbf{B}$ is not important for the performance of the algorithms. For instance, MSBL is more frequently attracted to local minima than our proposed algorithms, as experiments show later.

V-B Analysis of the Local Minima

In this subsection we discuss the local minimum property of the cost function $\mathcal{L}$ in (8) with respect to $\boldsymbol{\gamma}\triangleq[\gamma_{1},\cdots,\gamma_{M}]$, in which $\mathbf{\Sigma}_{0}=\mathbf{\Gamma}\otimes\mathbf{B}$ for fixed $\mathbf{B}$. Before presenting our results, we provide two lemmas needed to prove the results.

$\log|\mathbf{\Sigma}_{y}|\triangleq\log|\lambda\mathbf{I}+\mathbf{D}\mathbf{\Sigma}_{0}\mathbf{D}^{T}|$ is concave with respect to $\boldsymbol{\gamma}$.

This can be shown using the composition property of concave functions .

$\mathbf{y}^{T}\mathbf{\Sigma}_{y}^{-1}\mathbf{y}$ equals a constant $C$ when $\mathbf{\boldsymbol{\gamma}}$ satisfies the linear constraints

where $\mathbf{A}$ is full row rank, $\mathbf{1}_{L}$ is an $L\times 1$ vector of ones, and $\mathbf{u}$ is any fixed vector such that $\mathbf{y}^{T}\mathbf{u}=C$.

The proof is given in the Appendix. According to the definition of basic feasible solution (BFS) , we know that if $\boldsymbol{\gamma}$ satisfies Equation (34), then it is a BFS to (34) if $\|\boldsymbol{\gamma}\|_{0}\leq NL$, or a degenerate BFS to (34) if $\|\boldsymbol{\gamma}\|_{0}<NL$. Now we give the following result:

Every local minimum of the cost function $\mathcal{L}$ with respect to $\boldsymbol{\gamma}$ is achieved at a solution with $\|\widehat{\boldsymbol{\gamma}}\|_{0}\leq NL$, regardless of the values of $\lambda$ and $\mathbf{B}$.

Admittedly, the bound on the local minima $\|\widehat{\boldsymbol{\gamma}}\|_{0}$ is loose, and it is not meaningful when $NL>M$. However, we empirically found that $\|\widehat{\boldsymbol{\gamma}}\|_{0}$ actually is very smaller than $NL$.

Now, we calculate the local minima of the cost function $\mathcal{L}$. The result can provide some insights to the role of $\mathbf{B}$. Particularly, we are more interested in the local minima satisfying $\|\widehat{\boldsymbol{\gamma}}\|_{0}\leq N$, since the global minimum satisfies $\|\widehat{\boldsymbol{\gamma}}\|_{0}<N$. For these local minima, we have the following result:

In noiseless cases ($\lambda\rightarrow 0$), for every local minimum of $\mathcal{L}$ that satisfies $\|\widehat{\boldsymbol{\gamma}}\|_{0}\triangleq K\leq N$, its $i$-th nonzero element is given by $\widehat{\gamma}_{(i)}=\frac{1}{L}\widetilde{\mathbf{X}}_{i\cdot}{\mathbf{B}}^{-1}\widetilde{\mathbf{X}}^{T}_{i\cdot}\,(i=1,\cdots,K)$, where $\widetilde{\mathbf{X}}_{i\cdot}$ is the $i$-th nonzero row of $\widehat{\mathbf{X}}$ and $\widehat{\mathbf{X}}$ is the basic feasible solution to $\mathbf{Y}=\mathbf{\Phi}\mathbf{X}$.

From this lemma we immediately have the closed form of the global minimum.

$\mathbf{B}$ actually plays a role of temporally whitening the sources during the learning of $\boldsymbol{\gamma}$. To see this, assume all the sources have the same correlation structure, i.e. share the same matrix $\mathbf{B}$. Let $\mathbf{Z}_{i\cdot}\triangleq\widetilde{\mathbf{X}}_{i\cdot}\mathbf{B}^{-1/2}$. From Lemma 3, at the global minimum we have $\widehat{\gamma}_{(i)}=\frac{1}{L}\mathbf{Z}_{i\cdot}\mathbf{Z}^{T}_{i\cdot}\,(i=1,\cdots,K_{0})$. On the other hand, in the case of i.i.d. sources, at the global minimum we have $\widehat{\gamma}_{(i)}=\frac{1}{L}\widetilde{\mathbf{X}}_{i\cdot}\widetilde{\mathbf{X}}^{T}_{i\cdot}\,(i=1,\cdots,K_{0})$. So the results for the two cases have the same form. Since $E\{\mathbf{Z}_{i\cdot}^{T}\mathbf{Z}_{i\cdot}\}=\gamma_{i}\mathbf{I}$, we can see in the learning of $\boldsymbol{\gamma}$, $\mathbf{B}$ plays the role of whitening each source. This gives us a motivation to modify most state-of-the-art iterative reweighted algorithms by temporally whitening the estimated sources during iterations .

VI Computer Experiments

In our experiments we compared our T-SBL and T-MSBL with the following algorithms:

MSBL, proposed in The MATLAB code was downloaded at http://dsp.ucsd.edu/~zhilin/MSBL_code.zip.;

MFOCUSS, the regularized M-FOCUSS proposed in . In all the experiments, we set its p-norm $p=0.8$, as suggested by the authors The MATLAB code was downloaded at http://dsp.ucsd.edu/~zhilin/MFOCUSS.m.;

Set the iteration count $k$ to zero and $w_{i}^{(0)}=1,i=1,\cdots,M$

Update the weights for each $i=1,\cdots,M$

where $\epsilon^{(k)}$ is adaptively adjusted as in ;

For reproducibility, the experiment codes can be downloaded at http://dsp.ucsd.edu/~zhilin/TSBL_code.zip.

Figure 1 shows that with $L$ increasing, all the algorithms had better performance. But as $|\beta|\rightarrow 1$, for all the compared algorithms the benefit from multiple measurement vectors diminished. One surprising observation is that our T-MSBL and T-SBL had excellent performance in all cases, no matter what the temporal correlation was. Notice that even sources had no temporal correlation ($\beta=0$), T-MSBL and T-SBL still had better performance than MSBL.

Since the performance of all the algorithms at a given correlation level $\beta$ is the same as their performance at the correlation level $-\beta$, in the following we mainly show their performance at positive correlation levels.

VI-B Recovered Source Number at Different Temporal Correlation Levels

In this experiment we study the effects of temporal correlation on the number of accurately recovered sources in a noiseless case. The dictionary matrix $\mathbf{\Phi}$ was of the size $25\times 125$. $L$ was 4. $K$ varied from 10 to 18. The sources were generated in the same manner as before. Algorithms were compared at four different temporal correlation levels, i.e. $\beta=0$, $0.5$, $0.9$, and $0.99$. Results (Fig.3) show that T-MSBL and T-SBL accurately recovered much more sources than other algorithms, especially at high temporal correlation levels. This indicates that our proposed algorithms are very advantageous in the cases when the source number is large.

VI-C Ability to Handle Highly Underdetermined Problem

Most published works only compared algorithms in mildly underdetermined cases, namely, the ratio of $M/N$ was about $2\sim 5$. However, in some applications such as neuroimaging, one can easily have $N\approx 100$ and $M\approx 100000$. So, in this experiment we compare the algorithms in the highly underdetermined cases when $N$ was fixed at 25 and $M/N$ varied from 1 to 25. The source number $K$ was 12, and the measurement vector number $L$ was 4. SNR was 25 dB. Different to previous experiments, all the sources were AR(1) processes but with different AR coefficients. Their AR coefficients were uniformly chosen from $(0.5,1)$ at random. Results are presented in Fig.4, from which we can see that when $M/N\geq 10,$ all the compared algorithms had large errors. In contrast, our proposed algorithms had much lower errors. Note that due to the performance trade-off between $N$ and $M$, if one increases $N$, algorithms can keep the same recovery performance for larger $M/N$.

VI-D Recovery Performance for Different Kinds of Sources

In previous experiments all the sources were AR(1) processes. Although we have pointed out that for small $L$ modeling sources by AR(1) processes is sufficient, here we carry out an experiment to show our algorithms maintaining the same superiority for various time series. Since from previous experiments we have seen that T-SBL has similar performance to T-MSBL, and that MSBL has the best performance among the compared algorithms, in this experiment we only compare T-MSBL with MSBL.

VI-E Recovery Ability at Different Noise Levels

From previous experiments we have seen that the proposed algorithms significantly outperformed all the compared algorithms in noiseless scenarios and mildly noisy cases, even though to derive T-MSBL we used the approximation () which takes the equal sign only when $\mathbf{B}=\mathbf{I}$ (no temporal correlation) or $\lambda=0$ (no noise). Some natural questions may be raised: What is the performance of T-SBL and T-MSBL in strongly noisy cases? Is it still beneficial to exploit temporal correlation in these cases? To answer these questions, we carry out the following experiment.

Note that in low SNR cases, the estimated $\mathbf{B}$ of T-SBL and T-MSBL can include large errors, and thus the estimated amplitudes of sources are distorted. To reduce the distortion, we set $\mathbf{B}=\mathbf{I}$ once the number of nonzero $\gamma_{i}$ was less than $N$ during the learning procedure. The reason is that the role of $\mathbf{B}$ is to prevent T-SBL/T-MSBL from arriving at local minima; once the algorithms approach global minima very closely, $\mathbf{B}$ is no longer useful.

Also note that the $\lambda$ learning rules of T-SBL, T-MSBL and MSBL may not lead to optimal performance in low SNR cases. To avoid the potential disturbance of these $\lambda$ learning rules, we provided the three SBL algorithms with the optimal $\lambda^{*}$’s, which were obtained by the exhaustive search method stated previously.

Figure 6 shows that T-SBL and T-MSBL outperformed other algorithms in all the noise levels. This implies that even in low SNR cases exploiting temporal correlation of sources is beneficial.

But we want to emphasize that although the $\lambda$ learning rules of the three SBL algorithms may not be optimal in low SNR cases, our proposed $\lambda$ learning rules can lead to near-optimal performance, compared to the one of MSBL. To see this, we ran T-MSBL and MSBL again, but this time both algorithms used their $\lambda$ learning rules. T-MSBL used the modified version of the $\lambda$ learning rule (), i.e. setting the off-diagonal elements of $\mathbf{\Phi}\mathbf{\Gamma}\mathbf{\Phi}^{T}$ to zeros. The results (Fig. 6) show that MSBL had very poor performance when using its $\lambda$ learning rule. In contrast, T-MSBL’s performance was very close to its performance when using its optimal $\lambda^{*}$ T-SBL had the same behavior. But for clarity we do not present its performance curve.. The results indicate our proposed algorithms are advantageous in practical applications, since in practice the optimal $\lambda^{*}$’s are difficult to obtain, if not impossible.

VI-F Temporal Correlation: Beneficial or Detrimental?

From previous experiments one may think that temporal correlation is always harmful to algorithms’ performance, at least not helpful. However, in this experiment we will show that when SNR is high, the performance of our proposed algorithms increases with increasing temporal correlation.

The results indicating that temporal correlation is helpful may appear counterintuitive at first glance. A closer examination of the sparse recovery problems indicates a plausible explanation. There are two elements to the sparse recovery task; one is the location of the nonzero entries and the other is the value for the nonzero entries. Both tasks interact and combine to determine the overall performance. Correlation helps the estimation of the values for the nonzero entries and this may be important for the problem when dealing with finite matrices and may be lost when dealing with limiting results as the matrix dimension go to infinity. A more rigorous study of the interplay between estimation of the values and estimation of the locations is an interesting topic.

VI-G An Extreme Experiment on the Importance of Exploiting Temporal Correlation

It may be natural to take for granted that in noiseless cases, when source vectors are almost identical, algorithms have almost the same performance as in the case when only one measurement vector is available. In the following we show that it is not the case.

These results emphasize the importance of exploiting the temporal correlation, and also motivate future theoretical studies on the temporal correlation and the ill-condition issue of source matrices.

VII Discussions

Although there are a few works trying to exploit temporal correlation in the MMV model, based on our knowledge no works have explicitly studied the effects of temporal correlation, and no existing algorithms are effective in the presence of such correlation. Our work is a starting point in the direction of considering temporal correlation in the MMV model. However, there are many issues that are unclear so far. In this section we discuss some of them.

In our algorithm development we used one single matrix $\mathbf{B}$ as the covariance matrix (up to a scalar) for each source model in order to avoid overfitting. Mathematically, it is straightforward to extend our algorithms to use multiple matrices to capture the covariance structures of sources. For example, one can classify sources into several groups, say $G$ groups, and the sources in a group are all assigned by a common matrix $\mathbf{B}_{i}$ ($i=1,\cdots,G,\;G\ll M$) as the covariance matrix (up to a scalar). It seems that this extension can better capture the covariance structures of sources while still avoiding overfitting. However, we find that this extension (even for $G=2$) has much poorer performance than our proposed algorithms and MSBL. One possible reason is that during the early stage of the learning procedure of our algorithms, the estimated sources from each iteration are far from the true sources, and thus grouping them based on their covariance structures is difficult, if not impossible. The grouping error may cause avalanche effect, leading to the noted poor performance. Reducing the grouping error and more accurately capturing the temporal correlation structures is an area for future work.

On the other hand, there may be many ways to parameterize and estimate $\mathbf{B}$. In this work we provide a general method to estimate $\mathbf{B}$. In we proposed a method to parameterize $\mathbf{B}$ by a hyperparameter $\beta$, i.e.,

which equivalently assumes the sources are AR(1) processes with the common AR coefficient $\beta$. The resulting algorithms have good performance as well. Also, for low SNR cases in our experiments, we added an identity matrix (with a scalar) to the estimated $\mathbf{B}$ in T-MSBL, and achieved satisfying performance. All these imply that $\mathbf{B}$ could have many forms. Finding the forms that are advantageous in strongly noisy environments is an important issue and needs further study.

VII-B The Parameter λ\lambda: Noise Variance or Regularization Parameter?

Some works considered alternative noise covariance models. In the authors assumed that the covariance matrix of multi-channel noise is $\lambda\mathbf{C}$, instead of $\lambda\mathbf{I}_{N}$, where $\mathbf{C}$ is a known positive definite and symmetric matrix and $\lambda$ is an unknown noise-variance parameter. This model may better capture the noise covariance structures, but generally one does not know the exact value of $\mathbf{C}$. Thus there is no clear benefit from this covariance model. In , instead of deriving a learning rule for the noise covariance inside the SBL framework, the authors estimated the noise covariance by a method independent of the SBL framework. But this method is based on a large number of measurement vectors, and has a high computational load.

On the other hand, due to the works in , which connected SBL algorithms to traditional convex relaxation methods such as Lasso and Basis Pursuit Denoising , it was found that $\lambda$ is functionally the same as the regularization parameters of those convex relaxation algorithms. This suggests the use of methods such as the modified L-curve procedure or the cross-validation to choose $\lambda$ especially in strongly noisy environments. It is also interesting to see that SBL algorithms could adopt the continuation strategies , used in Lasso-type algorithms, to adjust the value of $\lambda$ for better recovery performance or faster speed.

However, if some channels contain very large noise (e.g. outliers) and the number of such channels is very small, then as suggested in , we can extend the dictionary matrix $\mathbf{\Phi}$ to $[\mathbf{\Phi},\mathbf{I}]$ and perform any sparse signal recovery algorithms without modification. The estimated ‘sources’ associated with the identity dictionary matrix are these large noise components.

VII-C Connections to Other Models

In fact, our bSBL framework is a block sparsity model , and thus the derived T-SBL algorithm can be directly used for this model. Compared to most existing algorithms derived in this model , an important difference is that T-SBL considers the correlation within each block.

The time-varying sparsity model is another related model. Different to our MMV model that assumes the support of each source vector is the same, the time-varying sparsity model assumes the support is slowly time-varying. It is interesting to note that this model can be approximated by concatenation of several MMV models, where in each MMV model the support does not change. Thus our proposed T-SBL and T-MSBL can be used for this model. The results are appealing, as shown in our recent work .

It should be noted that the proposed algorithms can be directly used for the SMV model. In this case the matrix $\mathbf{B}$ reduces to a scalar, and the $\gamma_{i}$ learning rules are the same as the one in the basic SBL algorithm . But due to the effective $\lambda$ learning rules, our algorithms are superior to the basic SBL algorithm, especially in noisy cases.

VIII Conclusions

Acknowledgement

Z.Z would like to thank Dr. David Wipf for his considerable help with the study of SBL, Ms. Jing Wan for kind help in performing some experiments, Mr. Tim Mullen for kind help in the paper writing, Dr. Rafal Zdunek for providing the code of SOB-MFOCUSS, and Mr. Md Mashud Hyder for providing the code of ISL0. The authors thank the reviewers for their helpful comments and especially thank a reviewer for the idea of using multiple covariance matrices, which is discussed in Section VII.A.

Appendix

Since the proof is a generalization of the Theorem 1 in , we only give an outline.

It can be shown that when $\lambda\rightarrow 0$ (noiseless case), the above problem is equivalent to

So we only need to show the global minimizer of (37) satisfies the property stated in the theorem.

Assume in the noiseless problem $\mathbf{Y}=\mathbf{\Phi}\mathbf{X}$, $\mathbf{\Phi}$ satisfies the URP condition . For its any solution $\widehat{\mathbf{X}}$, denote the number of nonzero rows by $K$. Thus following the method in , we can show the above $g(\mathbf{x})$ satisfies

VIII-B Proof of Lemma 2

VIII-C Proof of Theorem 2

The proof follows along the lines of Theorem 2 in using our Lemma 1 and Lemma 2. Consider the optimization problem:

which indicates $\|\boldsymbol{\gamma}\|_{0}\leq NL$.

VIII-D Proof of Lemma 3

Letting $\frac{\partial\mathcal{L}(\boldsymbol{\gamma})}{\partial\widetilde{\gamma}_{i}}=0$ gives

The second derivative of $\mathcal{L}$ at $\widetilde{\gamma}_{i}=\frac{1}{L}\widetilde{\mathbf{x}}_{i}^{T}{\mathbf{B}}^{-1}\widetilde{\mathbf{x}}_{i}$ is given by

Since $\mathbf{B}$ is positive definite and $\widetilde{\mathbf{x}}_{i}\neq\mathbf{0}$, $\frac{\widetilde{\mathbf{x}}_{i}^{T}{\mathbf{B}}^{-1}\widetilde{\mathbf{x}}_{i}}{\widetilde{\gamma}_{i}^{3}}>0$. So $\widetilde{\gamma}_{i}=\frac{1}{L}\widetilde{\mathbf{x}}_{i}^{T}{\widehat{\mathbf{B}}}^{-1}\widetilde{\mathbf{x}}_{i}\;(i=1,\cdots,K)$ is a local minimum.

References