SBL-based multi-task algorithms for recovering block-sparse signals with unknown partitions

We consider in this paper the problem of reconstructing block-sparse signals with unknown block partitions. In the first part of this work, we extend the block-sparse Bayesian learning (BSBL) originally developed for recovering a single block-sparse signal in a single compressive sensing (CS) task scenario to the case of multiple CS tasks. A new multi-task signal recovery algorithm, called the extended multi-task block-sparse Bayesian learning (EMBSBL), is proposed. EMBSBL exploits the statistical correlation among multiple signals as well as the intra-block correlation within individual signals to improve performance. Besides, it does not need a priori information on block partition. As the second part of this paper, we develop the EMBSBL-based synthesized multi-task signal recovery algorithm, namely SEMBSBL, to make it applicable to the single CS task case. The idea is to synthesize new CS tasks from the single CS task via circular-shifting operations and utilizes the minimum description length principle to determine the proper set of the synthesized CS tasks for signal reconstruction. SEMBSBL can achieve better signal reconstruction performance over other algorithms that recover block-sparse signals individually. Simulations corroborate the theoretical developments.


Introduction
Compressive sensing (CS) enables reconstructing a signal that is sparse in a certain domain from its measurements obtained at a rate significantly lower than the Nyquist frequency [1]. If in addition to sparsity, the signal representation is also structured in the form of clustered non-zeros, the signal would be referred to as being block-sparse. In practice, block-sparsity can be found in multi-band signals [2] or in the measurements of gene expression levels [3]. It has been shown that exploring the block-sparsity enables robust signal recovery from fewer compressive measurements [4]. We shall consider in this paper the efficient recovery of block-sparse signals.
Several block-sparse signal reconstruction algorithms have been developed in literature. Based on the compressive sampling matching pursuit (CoSaMP) [5], the block compressive sampling matching pursuit (BCoSaMP) was proposed in [4]. It utilizes the knowledge on the number of non-zero blocks to achieve signal recovery. On the basis of the orthogonal matching pursuit (OMP) [6], the block *Correspondence: wyinggui@gmail.com 1 College of Electronic Science and Engineering, National University of Defense Technology, Deya Road, Changsha 410073, People's Republic of China Full list of author information is available at the end of the article orthogonal matching pursuit (BOMP) was developed in [7]. Zou et al. proposed a block fixed-point continuation algorithm in [8] for block-sparse signal recovery. Elhamifar and Vidal approached the problem via the application of convex relaxation and convex optimization [9]. The two methods developed in [8] and [9] require the availability of the information on the block size. In [10], the dictionary optimization for block-sparse signal representation was studied and the work assumed that the maximum block length was known. More recently, CluSS-MCMC [11] and BM-MAP-OMP [12] have been proposed, which require little a priori knowledge on the block partition. On the basis of Bayesian sparse learning for temporally correlated signals [13,14] proposed two block-sparse signal recovery algorithms, the block-sparse Bayesian learning (BSBL) and its extended version, the EBSBL algorithm. The BSBL algorithm utilizes the block partition information while EBSBL handles signals with unknown block partitions. Most techniques reviewed above fall under the category of the single-task CS, where the focus is on recovering a block-sparse signal from its compressive measurements.
The contribution of this paper is twofold. We shall first consider the block-sparse signal reconstruction in a http://asp.eurasipjournals.com/content/2014/1/14 multi-task scenario, where the signals in different CS tasks are statistically correlated. The multi-task compressive sensing (MCS) was originally developed in [15]. Mathematically, we have L CS tasks where y i is the compressive measurement vector of the ith task and i is the M i × N measurement matrix (M i << N). θ i is the original signal in the ith task to be recovered and n i represents the measurement noise.
In MCS, the correlation among θ i is explored so that θ i are reconstructed jointly. MCS outperforms the singletask CS algorithm in terms of the reduced number of compressive measurements needed for efficient signal recovery. However, existing MCS techniques do not take into account the structural information in signals, such as block-sparsity. We shall therefore propose in this paper an extended version of the EBSBL algorithm from [14]. The original EBSBL method does not assume the knowledge on the block partition information, and it exhibits better block-sparse signal reconstruction performance over methods such as CluSS-MCMC and BM-MAP-OMP. We shall generalize EBSBL to the MCS scenario and obtain a new technique, referred to as extended multi-task blocksparse Bayesian learning (EMBSBL). Besides using the statistical correlation among θ i as in MCS, EMBSBL also utilizes the intra-block correlation within each signal to improve performance. Simulations show that the block-sparse signal recovery performance of EMBSBL is superior to that of the benchmark algorithms.
When there is only one CS task, the proposed EMB-SBL algorithm would become inapplicable. To address this problem, in the second part of this work, we shall augment EMBSBL with the concept of the synthesized multi-taskbased signal recovery. The new algorithm is referred to as SEMBSBL in the rest of the paper. SEMBSBL first synthesizes multiple CS tasks from the single-task CS and then applies EMBSBL to recover the block-sparse signal. The multiple CS tasks are produced via simply circularshifting the columns of the measurement matrix of the original CS model, which corresponds to circular-shifting the elements in the original signal vector and creates signals that have overlapping clusters, or equivalently speaking, correlated signals. The number of synthesized tasks is determined by the minimum description length (MDL) principle. With increase in the computational complexity, the newly proposed SEMBSBL technique outperforms the previously developed block-sparse signal recovery methods in terms of significantly reduced reconstruction errors and the removal of the needs for detailed information on the sparsity structure. Computer simulations are provided to demonstrate the good performance of the proposed SEMBSBL method.
The remainder of this paper is organized as follows. Section 2 presents the new EMBSBL algorithm for recovering multiple correlated block-sparse signals jointly. Section 3 illustrates the idea of synthesizing multiple CS tasks from a single one and presents the proposed SEM-BSBL algorithm. Simulation results are given in Section 4 and Section 5 concludes the paper.

EMBSBL algorithm
The development of EMBSBL starts with extending BSBL-BO in [14] to the case of multiple CS tasks. The resulting algorithm, called MBSBL, can jointly recover block-sparse signals when their non-zero blocks all have the same size. We next generalize MBSBL to obtain EMBSBL that does not need the information on the signal sparsity structure.

MBSBL
Let S be the block size and K be the number of blocks in every signal to be recovered. If the measurement noise n i in (1) follows an i.i.d. Gaussian distribution with zero mean and covariance matrix β −1 I, the conditional likelihood function of y i is where N y i | i θ i , β −1 I represents a Gaussian distribution with mean i θ i and covariance matrix β −1 I. In BSBL, each block θ i,j ∈ R S×1 is assumed to satisfy a zero-mean multivariate Gaussian distribution If we further assume that blocks are mutually uncorrelated, the prior for θ i is given by Here, B j is a positive definite matrix, capturing the correlation structure within the jth block, and γ j is a nonnegative parameter controlling the block-sparsity of θ i,j . When γ j = 0, the jth block becomes zero. During the learning process, most γ j tend to be zero, due to the mechanism of automatic relevance determination [13].
To avoid overfitting, we set B j = B, j = 1, . . . , K. Thus, 0 = ⊗ B, where = diag (γ 1 , . . . , γ K ) and ⊗ denotes the Kronecker product. The posterior distribution of θ i is then given by where From (5), we note that β, γ , and B are the sharing parameters of all CS tasks. To estimate them, let Y= y 1 , . . . , y L be the measurement set of the L CS tasks. We have The where would yield the estimates of the sharing parameters β, γ and B. We shall adopt the approach used in [14] to identify γ via the bound-optimization method, and find β and B via expectation maximization (EM).

Estimating γ
Maximizing (9) is equivalent to the minimization of For this purpose, we replace the term log (det (C i )) with an upper bound and apply a surrogate function for the term y T i C −1 i y i and then minimize their summation.
The upper bound of log (det (C i )) depends on its supporting hyperplane. Let γ * be a given point in the γ -space and we have , which corresponds to the jth block of θ i . We next introduce the surrogate function for the term y T i C −1 i y i . The purpose is to facilitate evaluating the partial derivatives of y T i C −1 i y i with respect to the sharing parameters β and γ . Originally, the sharing parameters appear in the inverse of the matrix C i (see the definition of C i under (9)). The surrogate function for It can be easily verified that the cost function to be minimized on the rightmost of (11) is the logarithm of the numerator of (5), p y i |θ i , β p (θ i |γ , B ), and the solution to the minimization problem is μ θ i defined in (6).
Putting (10) and (11) into (9), we have Let = {θ 1 , . . . , θ L } be the set of original signals from the L CS tasks. We can express the upper bound of Taking the partial derivative of G (γ , ) with respect to γ j and setting the result to zero yield the desired estimate of γ j , the jth element in γ , which is given by

Estimating B and β
The EM technique is used to find = {B, β}. We proceed by first treating θ i as hidden variables and then maximizing where (old) denotes the evaluated parameters in the previous iteration, Here, we only consider the terms relating to B in W (β, B) and use the notation The partial derivative of (17) with respective to B, which is symmetric and positive definite because it characterizes the covariance matrix of every signal block (see the definition given above (5)), is where μ Similar to [14], we improve the performance of the algorithm by restraining the matrix B. Specifically, we attempt to find a positive definite and symmetric matrixB to approximate B. Mathematically, we setB to be a Toeplitz matrix equal tô where r = m 1 m 0 , m 0 , and m 1 are obtained by averaging the elements along the main diagonal and the main sub-diagonal of B in (19). As a result, the approximated version of B is fully characterized by r. This method can also be applied with some modifications to the case where signal blocks have different sizes. In particular, in this case, we first computer =m 1 m 0 , wherem 0 = . B j are approximated withB j = Toeplitz 1,r, · · · ,r S j −1 such that again,B j depend on the value ofr only. Here, S j is the size of block j. http://asp.eurasipjournals.com/content/2014/1/14 We next evaluate β. Consider the terms relating to β in W (β, B) and use the notation Differentiating W 2 (β) with respect to β and then setting the result to zero, we obtain We have, after some manipulations, The iterative process for estimating β, γ , and B starts with initial solution guesses of μ θ i , θ i , γ j , B, and β. We then evaluate sequentially (6), (7), and (14) to find γ j and proceed to find the updated estimates of B and β using (19) and (22). With the obtained estimates of the sharing parameters, the original signals θ i of the L CS tasks can be reconstructed by following the MCS technique [15]. This completes the development of MBSBL.

EMBSBL
We shall present the new EMBSBL algorithm that is based on the developed MBSBL technique. Similar to [14], we first assume that all the blocks are of equal size h and the non-zero blocks are arbitrarily located. We will show via simulations that EMBSBL is not sensitive to the choice of h. There are p = N − h + 1 possible blocks in every signal θ i . The jth block starts at the jth element of θ i and continues until the j + h − 1 th element. All the non-zero elements of θ i lie within a subset of these blocks. From the analysis above, we can have the decomposition of θ i where N×h is a zero matrix except that the submatrix composed of its jth row to j + h − 1 th row is replaced by the identity matrix I, and E j is the same for every θ i . The CS model (1) can then be re-expressed as and The new CS model (24) has its signals with the property of block-sparsity and the intra-block correlation is explicit. z i can be recovered using MBSBL and by utilizing (23), the original signals θ i of the CS tasks can then be found, which finishes the development of EMB-SBL for recovering block-sparse signals under the MCS framework.

SEMBSBL algorithm
The EMBSBL cannot be directly applied to recover a single block-sparse signal in the single CS task scenario, due to its nature of being an MCS technique. We shall augment it with the idea of synthesized MCS to address this difficulty. CS task synthesis via circular-shifting operation is developed below. This section ends with the improved EMBSBL algorithm, namely synthesized EMBSBL (SEM-BSBL), which utilizes the MDL principle to determine the optimal number of synthesized CS tasks to achieve satisfactory signal recovery performance. Figure 1 illustrates synthesizing multiple CS tasks from a single one. The absence of measurement noise is assumed here to improve clarity. The original CS task is y 1 = 1 θ 1 , where θ 1 is the block-sparse signal to be recovered and it has two non-zero clusters (shadowed). The columns of the measurement matrix 1 corresponding to the nonzero elements in θ 1 are also shadowed for illustration. Figure 1 indicates that a new CS task can be synthesized from the original one by circularly shifting the columns of http://asp.eurasipjournals.com/content/2014/1/14 Figure 1 Synthesis of a new CS task via circular shifting. 1 to the right by one column. In this way, the new CS task has a new measurement matrix 2 and a new signal θ 2 whose the elements are generated by circularly shifting θ 1 downward by one sample. The new CS task has the same compressive measurements as the original one. We assume that this observation holds also for the case where measurement noise is present. Comparing θ 1 with θ 2 reveals that the locations of their non-zero elements have overlaps. This implies that θ 1 and θ 2 are correlated. This forms the basis for utilizing EMBSBL in blocksparse signal recovery. Additional CS tasks can be synthesized by following a similar approach but with different directions and shifting amounts of the circular-shifting operations.

Synthesis of multiple CS tasks
It can be expected that due to the block-sparsity of the signal to be recovered, the signals of some synthesized CS tasks may not be well correlated with others. In other words, they only have few overlapping non-zero elements. The utilization of these CS tasks in recovering the original signal via EMBSBL would lead to poor signal reconstruction performance. To address this problem, we propose to utilize the MDL principle to determine the number of synthesized CS tasks for the block-sparse signal reconstruction, as will be detailed in the following subsection.

SEMBSBL
This section presents the proposed SEMBSBL algorithm. We shall first provide a method for evaluating the signal recovery quality of EMBSBL for a given set of synthesized CS tasks. This is essential for selecting the optimal set of synthesized CS tasks for block-sparse signal recovery. For this purpose, we apply the MDL principle. Basically, it states that among a set of competing statistical models, the best model is the one having the minimum code length for the given data [16,17]. This is mathematically equivalent to solvingQ = arg min Q∈M CL y, Q , where M denotes the set of possible models and CL y, Q is the code length function. We set CL y, Q to be the Shannon code length [18], i.e., CL y, Q = − log 2 p y, Q , where p y, Q is the probability density function of y under the model Q.
For the problem of applying EMBSBL to recovering the block-sparse signal in a single CS task scenario (without loss of generality, we assume the task is y 1 = 1 + n 1 ), we denote the estimates of the sharing parameters β,γ , B as β, γ , B. They are output by the EMBSBL algorithm for a given set of synthesized CS tasks. The description length for y 1 , CL y 1 , can be then expressed as, after using (9) and setting L = 1, where CL y 1 We are now ready to present the proposed SEMBSBL algorithm. It is an iterative method that improves the signal recovery quality gradually. In each iteration, a new CS tasks is synthesized using the circular-shifting operation as illustrated in Figure 1. The newly produced task is applied together with the previously synthesized CS tasks as well as the original CS task in EMBSBL for jointly reconstructing the block-sparse signal. The above process continues until the number of synthesized CS tasks reaches a pre-specified value or including the newly synthesized CS task does not lead to better signal reconstruction quality (or equivalently, reduced code length for describing the data, which is given in (27)). http://asp.eurasipjournals.com/content/2014/1/14 The algorithm is summarized in Algorithm 1. l max is the user-specified maximum number of the synthesized CS tasks. EMBSBL l (Y,A) represents the application of EMBSBL for signal reconstruction in the lth iteration and it uses l CS tasks. Y and A collect the compressive measurements and their associated measurement matrices of the l CS tasks. The output of EMBSBL l (Y,A) is β l , γ l , B l , θ 1 l , which are the estimates of the sharing parameters β,γ , B and the original signal θ 1 (25) and (26), A); calculate CL l y 1 using (27).

Simulations
We shall provide simulation results to demonstrate the performance of the EMBSBL algorithm proposed in Section 2 and the SEMBSBL algorithm developed in Section 3. The signal reconstruction error is quantified using θ i −θ i 2 θ i 2 , where θ i andθ i are the true and the estimated signals. The elements of the measurement matrix i are initially drawn from the standard normal distribution N (0, 1) and each row of i is then normalized to have a unit norm.
In the first experiment, we simulate a two CS task scenario (L = 2) where the original signals both have a length of N = 500 and each contains 50 spikes with different amplitudes at random locations. The two signals also have six non-zero blocks with random sizes and they are at non-overlapping random locations. We consider two cases where 80% and 100% of the spikes of the two signals are at the same positions.
The signal-to-noise ratio (SNR) in log scale is defined as SNR = 20log 10 √ β i θ i 2 and zero-mean Gaussian noise is added to every CS measurement vector. We set SNR = 15 dB and use MCS, BSBL-BO, and EBSBL-BO as benchmark techniques for comparison. When implementing BSBL-BO, EBSBL-BO, and EMBSBL, we set the block size parameter h to be h = 4 and h = 8 to illustrate the impact of different choice of h on their performance. The signal reconstruction error results shown are obtained via averaging over 50 ensemble runs. Figure 2 compares MCS, BSBL-BO, EBSBL-BO, and EMBSBL in terms of their signal reconstruction errors as a function of the number of the compressive measurements. In this simulation, the intra-block correlation coefficient for each block is uniformly distributed between 0 and 0.1. Simulation results indicate that for BSBL-BO and EBSBL-BO, the performance curves when 80% and 100% of the spikes of the two original signals are at the same locations are very similar to each other. Therefore, to improve the clarity of the figures, we only provide in Figure 2 and the following Figures 3 and 4 the results when 80% of the spikes have the same locations.
We can see from Figure 2a that the proposed EMB-SBL has the least signal recovery error and its performance improves as more spikes of the original signals share the same locations. More importantly, EMBSBL is less sensitive to the choice of the block size parameter h than BSBL-BO and EBSBL-BO. This is because the new EMBSBL technique is an MCS algorithm that recovers the original signals of multiple CS tasks jointly. Compared with BSBL-BO and EBSBL-BO that are single-task CS algorithms, EMBSBL explores the inter-correlation among original signals to improve performance, besides the intra-block correlation. The use of this additional information improves the robustness of EMBSBL to the deviation of the presumed block size from the true value. Finally, as shown in Figure 2b, EMBSBL has the comparable complexity as EBSBL-BO, despite of being an MCS technique.
We repeat the simulation that produced Figure 2 but this time, we allow the intra-block correlation coefficient for each block to be uniformly distributed between 0.4 and 0.5. The obtained results are summarized in Figure 3. It can be seen that the proposed EMBSBL continues to offer the best signal recovery performance. Besides, the increase in the intra-block correlation improves the performance of the EBSBL-BO.
In producing Figure 4, we further increase the intrablock correlation coefficient for each block to be uniformly distributed from 0.8 to 0.9. The performance curve of MCS is excluded from the figure for sake of clarity, because in this case, the signal reconstruction error of MCS is large. We can find in Figure 4 that under significant intra-block correlation, the proposed EMBSBL still has the best signal recovery performance. Besides, comparing   it explicitly utilizes the intra-block correlation for better signal recovery (see Section 2).
We next study the impact of different signal intercorrelation levels on the signal recovery performance of EMBSBL. The simulation setup is the same as that leading to Figure 2 except that the percentage of overlapping nonzeros of the two original signals is set to be 40%, 60%, 80%, and 100%. Three sets of simulation results are produced. They correspond to the intra-block correlation coefficient being uniformly distributed within [0, 0.1], [0.4, 0.5], and [0.8, 0.9]. The obtained simulation results are summarized in Figure 5. We find that EMBSBL can recover the original signals with reduced signal reconstruction error as the inter-correlation among original signals increases. This observation holds for different intra-block correlation coefficients. The performance improvement is somewhat expected, because the EMBSBL algorithm is an MCS technique and its signal recovery performance would benefit from increased inter-correlation among original signals.
To validate the development of SEMBSBL (see Section 3), we consider a single CS task scenario. Again, the original signal has a length of N = 500 and it contains 50 spikes at random locations. Besides, the signal include six non-zero blocks with random sizes and random but non-overlapping locations. We set SNR = 15 dB in the simulation and use BSBL-BO and EBSBL-BO as benchmarking techniques. For BSBL-BO, EBSBL-BO, and the proposed SEMBSBL, we generate signal recovery error performance curves with the block size parameter setting to be h = 4 and h = 8. The results are averaged over 50 runs. For the proposed SEMBSBL, we set the prespecified maximum number of the synthesized CS tasks   Figure 6. It can be observed that the proposed SEMBSBL outperforms benchmark algorithms and SEMBSBL is less sensitive to the choice of the parameter h than EBSBL-BO. But the running time of SEMBSBL is high, since it executes the EMBSBL algorithm a couple of times before it finds an optimal set of synthesized CS tasks for signal recovery.

Conclusion
In this paper, a novel algorithm for jointly recovering multiple block-sparse signals from their compressive measurements, termed as the EMBSBL algorithm, was developed. EMBSBL exploits both the statistical correlation among signals and signals' intra-block correlation to achieve superior signal recovery performance. Moreover, the new algorithm eliminates the requirement on the availability of a priori information on the sparsity structure of the original signal. We also developed in this paper SEMBSBL that applies EMBSBL to the single CS task case. It synthesizes new CS tasks from the single CS task via simple circular-shifting operations to make EMBSBL applicable. The MDL principle was adopted to determine the proper set of the synthesized CS tasks for reconstructing the block-sparse signal. Computer simulations were carried out and revealed that the proposed EMBSEBL and SEMBSBL are able to outperform existing techniques in providing greatly enhanced block-sparse signal reconstruction performance at the cost of increased computational complexity.