# A Hybrid Orthogonal Forward-Backward Pursuit Algorithm for Partial Fourier Multiple Measurement Vectors Problem

Authors: Haiqiang Liu; Gang Hua; Aichun Zhu; Hongsheng Yin; Yonggang Xu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In solving the partial Fourier Multiple Measurement Vectors (FMMV) problem, existing greedy pursuit algorithms such as Simultaneous Orthogonal Matching Pursuit (SOMP), Simultaneous Subspace Pursuit (SSP), Hybrid Matching Pursuit (HMP), and Forward-Backward Pursuit (FBP) suffer from low recovery ability or need sparsity as a prior information. This paper combines SOMP and FBP to propose a Hybrid Orthogonal Forward-Backward Pursuit (HOFBP) algorithm. As an iterative algorithm, each iteration of HOFBP consists of two stages. In the first stage,α indices selected by SOMP are added to the support set. In the second stage, the support set is shrank by removing β indices. The choice of α and β is critical to the performance of this algorithm. The simulation results showed that, by using proper parameters, HOFBP has better performance than other greedy pursuit algorithms at the expense of more computing time in some cases. HOFBP does not need sparsity as a prior knowledge.


---

## Body

## 1. Introduction

Compressed Sensing (CS) [1–3] is a new information collection theory which has broken through the Nyquist sampling theorem. The Multiple Measurement Vector (MMV) problem [4] is a natural extension of CS to deal with multiple signals. It has potentials in many practical problems, such as underwater acoustic channel estimation [5], source localization [6], Direction-Of-Arrival (DOA) [7] as well as image super-resolution [8], face recognition [9], and image reconstruction [10]. In MMV problems, if the partial Fourier measurement matrix is used to measure signals, we have a partial Fourier MMV problem.Many algorithms were proposed to solve the MMV problem, such as Multiple Signal Classification algorithms [10–13], Bayesian learning algorithms [14–16], and greedy pursuit (GP) algorithms [17–20] as well as l1-SVD  [6] and q-thresholding algorithm [21]. Generally, GP algorithms are easily implemented and have low computational complexity. It is known that the GP algorithms are based on the iterative selection of the strongest components in the measurements projected onto the measurement matrix columns. Different from Gaussian random measurement matrix, partial Fourier matrix has the property that its adjoining atoms (columns of the measurement matrix) are strongly correlated [17, 22]. Because of this, greedy pursuit (GP) algorithms have some problems in solving the partial Fourier MMV problem. Simultaneous Orthogonal Matching Pursuit (SOMP) [18] lacks a mechanism to revaluate the selected indices. Once a wrong index is selected, there is no chance to remove it. Simultaneous Subspace Pursuit (SSP) [19] and Simultaneous Compressive Sampling Matching Pursuit (SCoSaMP) [20] select multiple atoms in each selection. This makes them weak in distinguishing adjoining atoms. G. Li and R. J. Burkholder [17] combine SOMP and SSP to propose the Hybrid Matching Pursuit (HMP) algorithm. However, in our simulations, HMP’s performance is not always good enough. Moreover, SSP, SCoSaMP, and HMP need sparsity as a prior information, which is impossible to know a prior for most practical applications. Forward-Backward Pursuit (FBP) [23] does not need the sparsity constraint as a prior information. However, it is also weak in distinguishing adjoining atoms. As stated in [24], the Bayesian-based CS reconstruction algorithms run in the opposite direction to GP algorithms. In the Bayesian reconstruction, the sparsity is reduced from its maximal value by using hyperparameters. In this case, if a component is missed, then the signal cannot be recovered. While the inclusion of a nonexisting component will only influence the noise in the resulting reconstruction, the omission of such a component will prevent correct reconstruction.In this paper, we combine SOMP and FBP to propose a Hybrid Orthogonal Forward-Backward Pursuit (HOFBP) algorithm. Inheriting the advantages of both SOMP and FBP, HOFBP has high-resolution in distinguishing adjoining atoms, as well as a mechanism to revaluate the selected indices, and does not require sparsity as a prior information. Simulation results show that, by using proper parameters, HOFBP has better performance than other GP algorithms at the expense of more computing time in some cases.The rest of this paper is organized as follows. In Section2, we introduce some notation and concepts. In Section 3, we introduce related GP algorithms, including SOMP, SSP, HMP, and FBP. In Section 4, we propose the HOFBP algorithm. In Section 5, the computational complexity of HOFBP is analyzed. In Section 6, we provide the simulation results. Conclusions are stated in Section 7.


## 2. Notation and Concepts

In this paper, we use the following notations.L denotes the number of signals. S≜1,2,…,N denotes the set of indices which correspond to the rows of the joint signal X. For the matrix A∈RM×N, Aj represents the j-th column of A, and A(i) represents the i-th row of A. For the set G={1,2,…,q}, G denotes the cardinality of G, AG represents the matrix composed of the columns Ajj∈G, and X(G) represents the matrix composed of the rows A(i)i∈G. Λ(X)≜1≤i≤N∣X(i)≠0 denotes the support set, which records the indices (positions of nonzero rows in X). If Λ(X)≤K≪M, X is called jointly sparse, and K is called the sparsity. In this paper, we assume Λ(X)=K. AH represents the conjugate transpose matrix of A. A†=(AHA)-1AH represents the pseudoinverse matrix of A. AF denotes the Frobenius norm of A. span(AG) denotes the space spanned by the columns of AG.


### 2.1. Single Measurement Vector (SMV)

The measurement process of CS can be denoted asy=Ax+n, where x=x1,x2,…,xNT∈RN is a sparse signal, A∈RM×N(M<N) is the measurement matrix, y=y1,y2,…,yMT∈RM is the measurement vector, and n=n1,n2,…,nMT∈RM is the additive noise. Given A and y, estimating x from y is named as the Single Measurement Vector (SMV) problem [1–3, 25]. l1 minimization algorithms [26–28] and GP algorithms are proposed to solve the SMV problem. There are many GP algorithms, such as Orthogonal Matching Pursuit (OMP) [29], Subspace Pursuit (SP) [30], Compressive Sampling Matching Pursuit (CoSaMP) [31], Forward-Backward Pursuit (FBP) [23] as well as Regularized Orthogonal Matching Pursuit (ROMP) [32], Stagewise Orthogonal Matching Pursuit (StOMP) [33], and Generalized Orthogonal Matching Pursuit (gOMP) [34].


### 2.2. Multiple Measurement Vector (MMV)

The MMV problem is an extension of the SMV problem. In the MMV problem, the measurement matrixA measures L signals, which can be denoted as(1)yl=Axl+nl,l=1,2,…,Lwhere xl∈RN represents the l-th signal, yl∈RM denotes the l-th measurement vector, A∈RM×N(M<N) denotes the measurement matrix, and nl denotes the additive noise. The measurement process can be rewritten as Y=AX+Noise, where Y=[y1,y2,…,yL], X=[x1,x2,…,xL], and Noise=[n1,n2,…,nL]. If only K rows of X are nonzero vectors, X is called jointly sparse and the sparsity is K. Supposing X is jointly sparse, the problem of recovering X from Y, which is also called the MMV problem, can be formulated as (2)minXΛXs.t.Y-AXF<ε,where ε≥0 is the error bound.


### 2.3. Partial Fourier Measurement Matrix

In MMV, many measurement matrices can be used to measure signals, such as Gaussian random matrix [1], partial Fourier matrix [2, 3], and Bernoulli random matrix [2] as well as deterministic measurement matrices [35, 36].If Gaussian random matrix is used to measure signals, we call the MMV problem as Gaussian random MMV problem. If the partial Fourier measurement matrix is used to measure signals, we call the MMV problem as partial Fourier MMV problem.The elements of the Gaussian random measurement matrix follow a Gaussian distribution with meanμ=0 and variance σ2=1/M. The Gaussian random measurement matrix can be defined as(3)Ap,s~N0,1M,p∈0,1,…,M-1,s∈0,1,…,N-1The Gaussian random measurement matrix is a random projection matrix [37–39], which is widely used in machine learning and dimensionality reduction.The discrete Fourier transform (DFT) matrix is defined as(4)Ψw,z=1Nexp⁡-i2πwzN,w,z∈0,1,…,N-1The partial Fourier measurement matrix is formed by randomly choosingM rows from Ψ, where M<N.The elements of Gaussian random measurement matrix follow Gaussian distribution. Its adjoining atoms are weakly correlated. The measurement matrix used for radar imaging is typically a Fourier-like measurement matrix. As illustrated in [17, 22], compared with Gaussian random measurement matrix, the adjoin atoms of the Fourier-like measurement matrix are strongly correlated.As GP algorithms iteratively select the atoms with the largest correlation to the measurement signal, it is harder for them to distinguish adjoining atoms of partial Fourier measurement matrix than that of Gaussian random measurement matrix. Our algorithm solves this problem by combining SOMP and FBP.What should be pointed out is that because Gaussian random measurement matrix has not the problem as partial Fourier measurement matrix, our algorithm has no advantage over some GP algorithms in solving the Gaussian random MMV problem. The experimental results in the Appendix illustrate this aspect.


## 2.1. Single Measurement Vector (SMV)

The measurement process of CS can be denoted asy=Ax+n, where x=x1,x2,…,xNT∈RN is a sparse signal, A∈RM×N(M<N) is the measurement matrix, y=y1,y2,…,yMT∈RM is the measurement vector, and n=n1,n2,…,nMT∈RM is the additive noise. Given A and y, estimating x from y is named as the Single Measurement Vector (SMV) problem [1–3, 25]. l1 minimization algorithms [26–28] and GP algorithms are proposed to solve the SMV problem. There are many GP algorithms, such as Orthogonal Matching Pursuit (OMP) [29], Subspace Pursuit (SP) [30], Compressive Sampling Matching Pursuit (CoSaMP) [31], Forward-Backward Pursuit (FBP) [23] as well as Regularized Orthogonal Matching Pursuit (ROMP) [32], Stagewise Orthogonal Matching Pursuit (StOMP) [33], and Generalized Orthogonal Matching Pursuit (gOMP) [34].


## 2.2. Multiple Measurement Vector (MMV)

The MMV problem is an extension of the SMV problem. In the MMV problem, the measurement matrixA measures L signals, which can be denoted as(1)yl=Axl+nl,l=1,2,…,Lwhere xl∈RN represents the l-th signal, yl∈RM denotes the l-th measurement vector, A∈RM×N(M<N) denotes the measurement matrix, and nl denotes the additive noise. The measurement process can be rewritten as Y=AX+Noise, where Y=[y1,y2,…,yL], X=[x1,x2,…,xL], and Noise=[n1,n2,…,nL]. If only K rows of X are nonzero vectors, X is called jointly sparse and the sparsity is K. Supposing X is jointly sparse, the problem of recovering X from Y, which is also called the MMV problem, can be formulated as (2)minXΛXs.t.Y-AXF<ε,where ε≥0 is the error bound.


## 2.3. Partial Fourier Measurement Matrix

In MMV, many measurement matrices can be used to measure signals, such as Gaussian random matrix [1], partial Fourier matrix [2, 3], and Bernoulli random matrix [2] as well as deterministic measurement matrices [35, 36].If Gaussian random matrix is used to measure signals, we call the MMV problem as Gaussian random MMV problem. If the partial Fourier measurement matrix is used to measure signals, we call the MMV problem as partial Fourier MMV problem.The elements of the Gaussian random measurement matrix follow a Gaussian distribution with meanμ=0 and variance σ2=1/M. The Gaussian random measurement matrix can be defined as(3)Ap,s~N0,1M,p∈0,1,…,M-1,s∈0,1,…,N-1The Gaussian random measurement matrix is a random projection matrix [37–39], which is widely used in machine learning and dimensionality reduction.The discrete Fourier transform (DFT) matrix is defined as(4)Ψw,z=1Nexp⁡-i2πwzN,w,z∈0,1,…,N-1The partial Fourier measurement matrix is formed by randomly choosingM rows from Ψ, where M<N.The elements of Gaussian random measurement matrix follow Gaussian distribution. Its adjoining atoms are weakly correlated. The measurement matrix used for radar imaging is typically a Fourier-like measurement matrix. As illustrated in [17, 22], compared with Gaussian random measurement matrix, the adjoin atoms of the Fourier-like measurement matrix are strongly correlated.As GP algorithms iteratively select the atoms with the largest correlation to the measurement signal, it is harder for them to distinguish adjoining atoms of partial Fourier measurement matrix than that of Gaussian random measurement matrix. Our algorithm solves this problem by combining SOMP and FBP.What should be pointed out is that because Gaussian random measurement matrix has not the problem as partial Fourier measurement matrix, our algorithm has no advantage over some GP algorithms in solving the Gaussian random MMV problem. The experimental results in the Appendix illustrate this aspect.


## 3. Greedy Pursuit Algorithms

In all greedy pursuit algorithms, a two-step strategy is adopted to solve problem (2). Firstly, the support set Λ is estimated. Then, the nonzero elements of the signal are estimated by using the least squares method, and the solution of problem (2) can be calculated as(5)xl^Λ=AΛ†ylandxl^S-Λ=0,where l=1,2,…,L and S≜1,2,…,N.To simplify the description of the following algorithms, we define two helping functions. The functionmax_indices() selects a certain number of indices which correspond to the largest amplitude entries of a vector. The function min_indices() selects a certain number of indices which correspond to the smallest amplitude entries of a vector.(6)max_indicesa,α≜αindiceswhichcorrespondtothelargestamplitudeentriesofthevectora.(7)min_indicesa,β≜βindiceswhichcorrespondtothesmallestamplitudeentriesofthevectora.SOMP is summarized in Algorithm1. In SOMP, the support set expands iteratively. In each iteration, the index corresponding to the atom which is most related to the residual vectors is added to the support set. If the sum of all the residual vectors’ l2-norms is larger than the termination threshold, we will terminate the iteration process, output the support set Λ, and then estimate the signals by using the least squares method as defined in (5); otherwise, the iteration process continues.Algorithm 1:SOMP.
Input:the measurement vectors yl where l = 1,2 , … , L, the measurement matrix A, the termination thresholdε. Initialization:Let the residual vectors rl=yl wherel=1,2,…,L, the support set Λ=⌀ and the number of iterationt=0. Iteration: (1) Select i~=argmaxj=1,2,…,N∑l=1LAjHrl and let Λ = Λ ∪ { i ~ }. (2) Update the residual vector rl=yl-AΛAΛ†yl wherel=1,2,…,L. (3) t=t+1. If ∑l=1Lrl2<ε, terminate the iteration; otherwise, go back to step (1). Output: Compute the estimated signals xl^ by (5).SSP is summarized in Algorithm2; it searches the support set iteratively. Each iteration is divided into index selection stage and backtracking stage. In the index selection stage, K indices corresponding to the atoms which are most related to all the residual vectors are added to the support set. In this way, the size of the support set expands to 2K (except on the first iteration). In the backtracking stage, L measurement vectors are individually projected onto the subspace spanned by the 2K atoms. Then, the support set is refined by selecting K indices which correspond to the largest entries in the fusion result of L projection coefficients. Once the sum of all the residual vectors’ l2-norms does not decrease, the iteration process is terminated and the signals can be estimated by (5).Algorithm 2:SSP.
Input:The measurement matrix A, the measurement vectorsyl where l=1,2,…,L, the sparsityK. Initialization: Let the residual vectors rl_old=yl wherel=1,2,…,L and the support set Λold=⌀. Iteration: (1) Let the temporary support set Λ t e m p = Λ o l d ∪ m a x _ i n d i c e s ( ∑ l = 1 L a b s ( A H r l _ o l d ) , K ). (2) Let Λnew=max_indices(∑l=1Labs(AΛtemp†yl),K). (3) The residual vector is updated as r l _ n e w = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) If ∑l=1Lrl_old2>∑l=1Lrl_new2, let Λold=Λnew, r l _ o l d = r l _ n e w where l=1,2,…,L and go back to step (1); otherwise, we terminate the iteration and let the final support setΛ=Λold. Output:Compute the estimated signals xl^ by (5).Similar to SSP, FBP is also a two-stage algorithm. It has two parameters: forward step sizeα and backward step size β, where β<α. In the first stage, α indices are added to the support set. In the second stage, β indices corresponding to the atoms which have the smallest contributions to projection coefficients are removed from the support set. Because β<α, the support set expands with the step size α-β. When the residual energy is less than a threshold, the iteration process is terminated. The choices of α and β are critical to the performance of FBP. According to the simulation results in [23], α∈[0.2M,0.3M] and α-β=1 lead to the best performance. In case of a failure, in order to avoid too many iterations, the size of the estimated support set is controlled by the parameter T, which is usually set as M.HMP combines the strengths of SOMP and SSP. It also searches the support set iteratively. Each iteration of HMP can be divided into index selection stage and backtracking stage. In the index selection stage, it utilizes the strategy of SOMP to select indices. Firstly,L solutions are obtained by individually applying OMP on L residual vectors. Then, K indices corresponding to the largest amplitude entries of L solutions’ fusion result are added to the support set. In the backtracking stage, HMP utilizes the revaluation mechanism of SSP to revaluate the selected indices. Firstly, L projection coefficients are obtained by individually projecting measurement vectors onto the subspace spanned by AΛtemp, where Λtemp is the support set. Then, the support set is refined by selecting K indices which correspond to the largest entries in the fusion result of the L projection coefficients. Once the sum of all the residual vectors’ l2-norms no longer decreases, the iteration process is terminated and the signals can be estimated by (5); otherwise, the iteration process continues.


## 4. HOFBP

### 4.1. Strengths and Drawbacks of SOMP, SSP, FBP, and HMP

From the description of SOMP and SSP, we can see that SOMP and SSP have both strengths and drawbacks in solving the partial Fourier MMV problem. SOMP selects one index in each iteration, which ensures the strong ability in distinguishing adjoining atoms. However, it lacks a mechanism to revaluate the selected indices. Once a wrong index is selected, it will lie in the support set forever.SSP has a mechanism to revaluate the selected indices. If the selected indices are incorrect, they will be removed from the support set in the next iteration. However, because the adjoining atoms of partial Fourier matrix are strongly correlated and SSP selects multiple atoms at once, it has a weak ability in distinguishing adjoining atoms. We observe that SSP performs worse than SOMP in solving the partial Fourier MMV problem.HMP combines the strengths of SOMP and SSP. On the one hand, HMP selects atoms by individually applying OMP on local residual vectors (Algorithm4). This ensures the atoms are selected one by one. On the other hand, in the backtracking stage, it utilizes the mechanism of SSP to revaluate the selected indices, which will remove the wrong selected indices from the support set. However, in this paper, we summarized the drawbacks of HMP. Firstly, our simulation results demonstrate that it does not always performs better than other algorithms; e.g., it performs worse than SOMP in the case of relative few signals. Moreover, it needs sparsity as a prior information which is almost impossible to be known in most practical applications. Although G. Li and R. J. Burkholder [17] propose a method to estimate sparsity, it is not always effective. Once the estimated sparsity is smaller than the real sparsity, the recovery process will fail.Being a two-step algorithm, FBP has mechanism to revaluate the selected atoms, and it does not need sparsity as a prior information (Algorithm3). However, in the forward step, it selects multiple atoms at once, then it is still weak in distinguishing adjoining atoms.Algorithm 3:FBP.
Input:The measurement vectors yl where l = 1,2 , … , L, the measurement matrix A, the forward step sizeα, the backward step size β, the terminate thresholdε, the limit of the size of the estimated support setT. Initialization:Letresidual vectors rl=yl where l = 1,2 , … , L, the support set Λold=⌀ and the number of iterationt=0. Iteration: (1) Let the temporary support set Λ t e m p = Λ o l d ∪ m a x _ i n d i c e s ( ∑ l = 1 L a b s ( A H r l _ o l d ) , α ). (2) Let Λ n e w = Λ t e m p - m i n _ i n d i c e s ( ∑ l = 1 L a b s ( A Λ t e m p † y l ) , β ). (3) Update the residual vector as r l = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) t=t+1. If ∑l=1Lrl2<ε or t>T, terminate the iteration and let the final support setΛ=Λold; otherwise, letΛold=Λnew and go back to step (1). Output:Compute the estimated signals xl^ by (5).Algorithm 4:HMP.
Input:The measurement matrix A, the measurement vectorsyl where l=1,2,…,L, the sparsityK. Initialization:Let the residual vectors rl_old=yl wherel=1,2,…,L and the support set Λold=⌀. Iteration: (1) Let the temporary support setΛtemp=Λold∪max_indices(∑l=1Labs(σl_omp),K) whereσl_omp=OMPrl,A,K is the recovered signal by applying OMP on thel-th residual vector. (2) Let Λnew=max_indices(∑l=1Labs(AΛtemp†yl),K). (3) The residual vector is updated as r l _ n e w = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) If ∑l=1Lrl_old2>∑l=1Lrl_new2, let Λold=Λnew, r l _ o l d = r l _ n e w and go back to step (1); otherwise, terminate the iteration and let the final support set Λ = Λ o l d. Output:Compute the estimated signals xl^ by (5).


### 4.2. Description of HOFBP

Inspired by the GP algorithms described above, we combine the strengths of SOMP and FBP to propose a Hybrid Orthogonal Forward-Backward Pursuit (HOFBP) algorithm.HOFBP is summarized in Algorithm5. In the initialization, the residual vectors are initialized as the measurement vectors and the support set is initialized as the empty set. In HOFBP, the support set expands iteratively. In the first step, as in HMP, L recovered signals are obtained by individually applying OMP on L residual vectors, then, α indices corresponding to the largest amplitude entries of the L recovered signals’ fusion result are added to the support set. In the second step, L projection coefficients are obtained by individually projecting measurement vectors onto the space spanned by AΛtemp, where Λtemp is the support set. Then, β indices which correspond to the minimum entries in the fusion result of the L projection coefficients are removed from the support set. In the third step, the residual vectors are updated. In the fourth step, the termination condition is checked. If the sum of all residual vectors’ l2-norms is smaller than the termination threshold or the size of the estimated support set exceeds the maximum allowed size of the support set, the iteration process is terminated, the support set obtained in the last iteration is taken as the final support set and the signals are estimated by using the least squares method as stated in (5). We should point out that this is because sparsity is assumed unknown in HOFBP, we take α as sparsity in the OMP algorithm.Algorithm 5:HOFBP.
Input:The measurement vectors yl where l = 1,2 , … , L, the measurement matrix A, the forward step sizeα, the backward step size β, the terminate thresholdε, the limit of the size of the estimated support setT. Initialization:Letresidual vectors rl=yl where l = 1,2 , … , L, the support set Λold=⌀ and the number of iterationt=0. Iteration: (1) Let the temporary support setΛtemp=Λold∪max_indices(∑l=1Labs(σl_omp),α) whereσl_omp=OMPrl,A,α is the recovered signal by applying the OMP algorithm on thel-th residual vector. (2) Let Λ n e w = Λ t e m p - m i n _ i n d i c e s ( ∑ l = 1 L a b s ( A Λ t e m p † y l ) , β ). (3) Update the residual vector as r l = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) t=t+1. If ∑l=1Lrl2<ε or t>T, terminate the iteration and let the final support setΛ=Λold; otherwise, letΛold=Λnew and go back to step (1). Output:Compute the estimated signals xl^ by (5).The index selection stage of HOFBP utilizes SOMP. Because OMP is a simplified version of SOMP for one signal, then HOFBP also seems closely related to OMP. HMP combines SOMP and SSP, and HOFBP combines SOMP and FBP. Because SSP and FBP are similar to each other, then, HOFBP also seems closely related to HMP. Although HOFPB seems closely related to OMP and HMP, because SOMP and FBP are basic algorithms of HOFPB, we still say that HOFBP is an algorithm combining SOMP and FBP.The Halting Conditions. The terminate threshold ε should be set as a very small positive number. In the noiseless case, we set ε=10-3. In noisy cases, it is set as ε=10-SNR/20, where SNR=10log10YF2/NoiseF2 is the signal-to-noise ratio (SNR). Moreover, in case of a failure, in order to avoid too many iterations, T is used to limit the size of the support set. In order to ensure all correct atoms are found before the support set size reaches T, it should be set larger than the real sparsity K. Because K is unknown, T can be set large enough or simply as M. If the sum of all residual vectors’ l2-norms is smaller than the terminate threshold ε or the number of iterations exceeds the maximum allowed number of iterations, the iteration process is terminated.The Choice of Forward and Backward Step Size. The forward step size α and backward step size β are two important parameters for HOFBP. The choices are critical to the performances of the algorithm. α should be set larger than 1 and smaller than the measurement number M. According to the research result of the paper [18], α∈[0.2M,0.3M] leads to the best performance. According to the definition in HOFBP, β should be set smaller than α. The smaller the value of α-β, the more carefully the indices are selected. When β=α-1, HOFBP has the best performance.


### 4.3. Relations to Other GP Algorithms

We can see that HOFBP has more strengths over other greedy pursuit algorithms. Firstly, compared with SOMP, it has a mechanism to revaluate the selected indices. Secondly, compared with SSP, it can distinguish adjoining atoms well by individually applying OMP on local residual vectors. Moreover, it does not need sparsity as a prior information. Thirdly, compared with HMP, the support set expands iteratively, which makes it does not need sparsity as a prior information. Lastly, compared with FBP, in the forward step, it selects atoms one by one, then it can distinguish adjoining atoms well.We noted the existence of two improved versions of FBP: Iterative Forward-Backward Pursuit (IFBP) [40] and Adaptive Forward-Backward Orthogonal Matching Pursuit (AFBOMP) [41]. However, our algorithm is very different from them. IFBP is an algorithm embedding FBP into the iterative framework. The performance improvement is mainly due to the iterative framework. In AFBOMP, an adaptive method of selecting α and β is proposed. By using this method, FBP’s performance is improved. Different from them, the main innovation of our algorithm is the atom selection strategy. The innovations of these three algorithms come from different aspects.


## 4.1. Strengths and Drawbacks of SOMP, SSP, FBP, and HMP

From the description of SOMP and SSP, we can see that SOMP and SSP have both strengths and drawbacks in solving the partial Fourier MMV problem. SOMP selects one index in each iteration, which ensures the strong ability in distinguishing adjoining atoms. However, it lacks a mechanism to revaluate the selected indices. Once a wrong index is selected, it will lie in the support set forever.SSP has a mechanism to revaluate the selected indices. If the selected indices are incorrect, they will be removed from the support set in the next iteration. However, because the adjoining atoms of partial Fourier matrix are strongly correlated and SSP selects multiple atoms at once, it has a weak ability in distinguishing adjoining atoms. We observe that SSP performs worse than SOMP in solving the partial Fourier MMV problem.HMP combines the strengths of SOMP and SSP. On the one hand, HMP selects atoms by individually applying OMP on local residual vectors (Algorithm4). This ensures the atoms are selected one by one. On the other hand, in the backtracking stage, it utilizes the mechanism of SSP to revaluate the selected indices, which will remove the wrong selected indices from the support set. However, in this paper, we summarized the drawbacks of HMP. Firstly, our simulation results demonstrate that it does not always performs better than other algorithms; e.g., it performs worse than SOMP in the case of relative few signals. Moreover, it needs sparsity as a prior information which is almost impossible to be known in most practical applications. Although G. Li and R. J. Burkholder [17] propose a method to estimate sparsity, it is not always effective. Once the estimated sparsity is smaller than the real sparsity, the recovery process will fail.Being a two-step algorithm, FBP has mechanism to revaluate the selected atoms, and it does not need sparsity as a prior information (Algorithm3). However, in the forward step, it selects multiple atoms at once, then it is still weak in distinguishing adjoining atoms.Algorithm 3:FBP.
Input:The measurement vectors yl where l = 1,2 , … , L, the measurement matrix A, the forward step sizeα, the backward step size β, the terminate thresholdε, the limit of the size of the estimated support setT. Initialization:Letresidual vectors rl=yl where l = 1,2 , … , L, the support set Λold=⌀ and the number of iterationt=0. Iteration: (1) Let the temporary support set Λ t e m p = Λ o l d ∪ m a x _ i n d i c e s ( ∑ l = 1 L a b s ( A H r l _ o l d ) , α ). (2) Let Λ n e w = Λ t e m p - m i n _ i n d i c e s ( ∑ l = 1 L a b s ( A Λ t e m p † y l ) , β ). (3) Update the residual vector as r l = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) t=t+1. If ∑l=1Lrl2<ε or t>T, terminate the iteration and let the final support setΛ=Λold; otherwise, letΛold=Λnew and go back to step (1). Output:Compute the estimated signals xl^ by (5).Algorithm 4:HMP.
Input:The measurement matrix A, the measurement vectorsyl where l=1,2,…,L, the sparsityK. Initialization:Let the residual vectors rl_old=yl wherel=1,2,…,L and the support set Λold=⌀. Iteration: (1) Let the temporary support setΛtemp=Λold∪max_indices(∑l=1Labs(σl_omp),K) whereσl_omp=OMPrl,A,K is the recovered signal by applying OMP on thel-th residual vector. (2) Let Λnew=max_indices(∑l=1Labs(AΛtemp†yl),K). (3) The residual vector is updated as r l _ n e w = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) If ∑l=1Lrl_old2>∑l=1Lrl_new2, let Λold=Λnew, r l _ o l d = r l _ n e w and go back to step (1); otherwise, terminate the iteration and let the final support set Λ = Λ o l d. Output:Compute the estimated signals xl^ by (5).


## 4.2. Description of HOFBP

Inspired by the GP algorithms described above, we combine the strengths of SOMP and FBP to propose a Hybrid Orthogonal Forward-Backward Pursuit (HOFBP) algorithm.HOFBP is summarized in Algorithm5. In the initialization, the residual vectors are initialized as the measurement vectors and the support set is initialized as the empty set. In HOFBP, the support set expands iteratively. In the first step, as in HMP, L recovered signals are obtained by individually applying OMP on L residual vectors, then, α indices corresponding to the largest amplitude entries of the L recovered signals’ fusion result are added to the support set. In the second step, L projection coefficients are obtained by individually projecting measurement vectors onto the space spanned by AΛtemp, where Λtemp is the support set. Then, β indices which correspond to the minimum entries in the fusion result of the L projection coefficients are removed from the support set. In the third step, the residual vectors are updated. In the fourth step, the termination condition is checked. If the sum of all residual vectors’ l2-norms is smaller than the termination threshold or the size of the estimated support set exceeds the maximum allowed size of the support set, the iteration process is terminated, the support set obtained in the last iteration is taken as the final support set and the signals are estimated by using the least squares method as stated in (5). We should point out that this is because sparsity is assumed unknown in HOFBP, we take α as sparsity in the OMP algorithm.Algorithm 5:HOFBP.
Input:The measurement vectors yl where l = 1,2 , … , L, the measurement matrix A, the forward step sizeα, the backward step size β, the terminate thresholdε, the limit of the size of the estimated support setT. Initialization:Letresidual vectors rl=yl where l = 1,2 , … , L, the support set Λold=⌀ and the number of iterationt=0. Iteration: (1) Let the temporary support setΛtemp=Λold∪max_indices(∑l=1Labs(σl_omp),α) whereσl_omp=OMPrl,A,α is the recovered signal by applying the OMP algorithm on thel-th residual vector. (2) Let Λ n e w = Λ t e m p - m i n _ i n d i c e s ( ∑ l = 1 L a b s ( A Λ t e m p † y l ) , β ). (3) Update the residual vector as r l = y l - A Λ n e w A Λ n e w † y l where l=1,2,…,L. (4) t=t+1. If ∑l=1Lrl2<ε or t>T, terminate the iteration and let the final support setΛ=Λold; otherwise, letΛold=Λnew and go back to step (1). Output:Compute the estimated signals xl^ by (5).The index selection stage of HOFBP utilizes SOMP. Because OMP is a simplified version of SOMP for one signal, then HOFBP also seems closely related to OMP. HMP combines SOMP and SSP, and HOFBP combines SOMP and FBP. Because SSP and FBP are similar to each other, then, HOFBP also seems closely related to HMP. Although HOFPB seems closely related to OMP and HMP, because SOMP and FBP are basic algorithms of HOFPB, we still say that HOFBP is an algorithm combining SOMP and FBP.The Halting Conditions. The terminate threshold ε should be set as a very small positive number. In the noiseless case, we set ε=10-3. In noisy cases, it is set as ε=10-SNR/20, where SNR=10log10YF2/NoiseF2 is the signal-to-noise ratio (SNR). Moreover, in case of a failure, in order to avoid too many iterations, T is used to limit the size of the support set. In order to ensure all correct atoms are found before the support set size reaches T, it should be set larger than the real sparsity K. Because K is unknown, T can be set large enough or simply as M. If the sum of all residual vectors’ l2-norms is smaller than the terminate threshold ε or the number of iterations exceeds the maximum allowed number of iterations, the iteration process is terminated.The Choice of Forward and Backward Step Size. The forward step size α and backward step size β are two important parameters for HOFBP. The choices are critical to the performances of the algorithm. α should be set larger than 1 and smaller than the measurement number M. According to the research result of the paper [18], α∈[0.2M,0.3M] leads to the best performance. According to the definition in HOFBP, β should be set smaller than α. The smaller the value of α-β, the more carefully the indices are selected. When β=α-1, HOFBP has the best performance.


## 4.3. Relations to Other GP Algorithms

We can see that HOFBP has more strengths over other greedy pursuit algorithms. Firstly, compared with SOMP, it has a mechanism to revaluate the selected indices. Secondly, compared with SSP, it can distinguish adjoining atoms well by individually applying OMP on local residual vectors. Moreover, it does not need sparsity as a prior information. Thirdly, compared with HMP, the support set expands iteratively, which makes it does not need sparsity as a prior information. Lastly, compared with FBP, in the forward step, it selects atoms one by one, then it can distinguish adjoining atoms well.We noted the existence of two improved versions of FBP: Iterative Forward-Backward Pursuit (IFBP) [40] and Adaptive Forward-Backward Orthogonal Matching Pursuit (AFBOMP) [41]. However, our algorithm is very different from them. IFBP is an algorithm embedding FBP into the iterative framework. The performance improvement is mainly due to the iterative framework. In AFBOMP, an adaptive method of selecting α and β is proposed. By using this method, FBP’s performance is improved. Different from them, the main innovation of our algorithm is the atom selection strategy. The innovations of these three algorithms come from different aspects.


## 5. Computational Complexity Analysis

SOMP and SSP are extended versions of OMP and SP. According to [30], both OMP and SP have the computational complexity of O(KMN). Because there are L measurement vectors in the MMV problem, both SOMP and SSP have computational complexity of O(KLMN). According to [17], the computational complexity of HMP mainly comes from the computations of applying OMP on L local residual vectors. Thus, each iteration of HMP has the computational complexity of O(KLMN). Since HMP performs O(K) iterations, it has the total computational complexity of O(K2LMN).In each iteration of HOFBP, the computational complexity is mainly composed of two parts. Part 1 mainly comes from applying OMP onL residual vectors. Because α is taken as sparsity, the operation of applying OMP on each residual vector has the complexity of O(αMN). Therefore, the operation of applying OMP on L residual vectors has the computational complexity of O(αLMN). Part 2 mainly comes from the projections of L measurement vectors on the space corresponding to the support set in step 2 of the iteration. In step 3 of the t-th(t>0) iteration, the size of support set is t(α-β), then, the computations of projecting L measurement vectors have the complexity of O(t2(α-β)2ML). Generally speaking, t(α-β)<α≪N, in comparison with the complexity of part 1, and the complexities of part 2 are marginal. Therefore, each iteration of HOFBP has the computational complexity of O(αLMN). HOFBP has a maximum of T/(α-β) iterations, then the total computational complexity of HOFBP is limited by O(TαLMN/(α-β)).


## 6. Simulation Results and Analysis

In this part, we conduct some numerical simulations to evaluate the performance and efficiency of HOFBP. The performance and efficiency of HOFBP are compared with those of HMP, SSP, SOMP, IFBP, and AFBOMP. SOMP, SSP, and HMP are traditional GP algorithms. IFBP and AFBOMP are improved versions of FBP. For IFBP, AFBOMP, and HOFBP, the forward step sizeα is set as 0.25M. Noisy and noiseless cases are both considered.We use the following hypothesis in the simulations.(I)
The measurement matrixA is constructed by randomly selecting M rows from a discrete Fourier matrix which has the size of N×N. Every column of A is normalized to unit l2 norm.(II)
Randomly generateL signals. They have same nonzero elements positions. K indices of nonzero elements positions are uniformly chosen from the set S={1,2,…,N}. The nonzero elements of each signal are independently chosen from the standard normal distribution.(III)
All curves in the following figures are obtained by averaging the results over200 independent trials. In each trial, the signals and measurement matrix A are generated independently.(IV)
Perfect recovery percentages are used to evaluate all algorithms’ performances. In thei-th trial, if X(i)-X^(i)2<0.05X(i)2, where X(i) and X^(i), respectively, denote the original joint signal and the recovered joint signal in the i-th trial, we announce the i-th trial is a perfect recovery. Supposing there is SE perfect recovery in 200 trials, the perfect recovery percentage is SE/200.(V)
Average run times are used to evaluate all algorithms’ efficiencies. They are computed by averaging the run times of200 independent trials.(VI)
The parameters of algorithms are set as follows: the limit of the size of the estimated support set in HOFBP, IFBP, and AFBOMP is set asT=M; the termination threshold of HOFBP, SOMP, IFBP, and AFBOMP is set as ε=10-3; the stop conditions of SSP and HMP are the real sparsity K.(VII)
The simulations are performed in Matlab R2014a running on a computer with 3.7 GHz, Intel Xeon E3 processor, with 16.0 GB of RAM, and Windows 7 system.


### 6.1. Robustness against the Step Size

Firstly, we analyze the influence of the value of step sizeα-β on the performances and efficiencies of HOFBP, IFBP, and AFBOMP. Figure 1 shows the perfect recovery percentages and average run times of these three algorithms against the sparsity K changing from 20 to 38 in the noiseless case. We fix other parameters as L=6, N=200, and M=40. α-β varies in 1,3,5,7.Figure 1
Perfect recovery percentages and average run times of HOFBP, IFBP, and AFBOMP with differentα-β against K with M=40, N=200, L=6.
(a) (b)For these three algorithms, we can see thatα-β=1 leads to the best performances. With the increase of the value of α-β, the performances of these three algorithms decline. That is because the smaller the value of α-β, the more carefully the indices are selected. For the efficiencies, we can see that the smaller the value of α-β is, the more time these three algorithms consume. The reason is that the smaller the value of α-β is, the more iterations are needed for the algorithms to reach convergence.We can see that, for the sameα-β, HOFBP has the highest recovery percentage. For efficiency, HOFBP consumes more time than AFBOMP, however, much less time than IFBP.


### 6.2. Robustness against Sparsity

Figure2 shows the perfect recovery percentages and average run times of all algorithms against the sparsity K changing from 20 to 38 in the noiseless case. We fix other parameters as L=6, N=200, and M=40. HMP performs better than SOMP at the beginning. With the increase of sparsity, HMP begins to perform worse than SOMP at the sparsity 24. SSP performs nearly consistently worse than SOMP, which is consistent with our previous analysis. HOFBP performs better than the other algorithms. IFBP has similar performance with HOFBP; however, it consumes much more time than HOFBP. For efficiency, HOFBP consumes similar time to HMP and AFBOMP but performs better than them.Figure 2
Perfect recovery percentages and average run times of different algorithms againstK with M=40, N=200, L=6.
(a) (b)


### 6.3. Robustness against the Measurement Number

Figure3 shows the perfect recovery percentages and average run times of all algorithms against the measurement number M changing from 34 to 60 in the noiseless case. We fix other parameters as L=6, N=200, and K=30. With the increase of the measurement number, the perfect recovery percentages of all algorithms increase. HMP and SSP perform consistently worse than SOMP. HOFBP has the best performance among all algorithms. IFBP performs similarly to HOFBP; however, it consumes much more time. For efficiency, HOFBP is similar to AFBOMP and HMP, but with better performance. Moreover, HOFBP needs the least measurement number to reach perfect recovery.Figure 3
Perfect recovery percentages and average run times of different algorithms againstM with K=30, N=200, L=6.
(a) (b)


### 6.4. Robustness against the Number of Signals

Figure4 shows the perfect recovery percentages and average run times of all algorithms against the number of signals L changing from 1 to 25 in the noiseless case. We fix other parameters as M=45, N=200, K=30. With the increase of L, the perfect recovery percentages of all algorithms increase. This is because the more measurement vectors are, the more correlations the algorithms can utilize. We can see that HMP performs worse than SOMP when L<11. This illustrates that HMP has less strengths than SOMP in the case of relative less signals. HOFBP performs better than other algorithms. Moreover, it needs less number of signals to reach perfect recovery. For efficiency, we have similar results to those of Sections 4.2 and 4.3.Figure 4
Perfect recovery percentages and average run times of different algorithms againstL with M=45, K=30, N=200.
(a) (b)


### 6.5. Recovery from Noisy Observation

In this part, we test algorithms’ performances against noise. We assume the noise is i.i.d Gaussian noise. The signal-to-noise ratio (SNR) is used to measure the relation between the measurement signal and the noise, which is defined as(8)SNR=10log10YF2NoiseF2.What should be pointed out is that, in the noisy case, the terminate parameters of SOMP, HOFBP, IFBP, and AFBOMP are set asε=10-SNR/20. The terminate parameters of SSP and HMP are still set as the real sparsity.Figure5 shows the perfect recovery percentages and average run times of all algorithms against the change of noise. Other parameters are fixed as M=45, K=30, N=200, and L=6. With the increase of SNR, the perfect recovery percentages of all algorithms increase. SSP and HMP perform consistently worse than SOMP. AFBOMP performs worst, which means it is not suitable for the noisy situation. HOFBP performs consistently better than other algorithms. For efficiency, we have similar result to those of Sections 4.2 and 4.3.Figure 5
Perfect recovery percentages and average run times of different algorithms againstSNR with M=45, K=30, N=200, and L=6.
(a) (b)


## 6.1. Robustness against the Step Size

Firstly, we analyze the influence of the value of step sizeα-β on the performances and efficiencies of HOFBP, IFBP, and AFBOMP. Figure 1 shows the perfect recovery percentages and average run times of these three algorithms against the sparsity K changing from 20 to 38 in the noiseless case. We fix other parameters as L=6, N=200, and M=40. α-β varies in 1,3,5,7.Figure 1
Perfect recovery percentages and average run times of HOFBP, IFBP, and AFBOMP with differentα-β against K with M=40, N=200, L=6.
(a) (b)For these three algorithms, we can see thatα-β=1 leads to the best performances. With the increase of the value of α-β, the performances of these three algorithms decline. That is because the smaller the value of α-β, the more carefully the indices are selected. For the efficiencies, we can see that the smaller the value of α-β is, the more time these three algorithms consume. The reason is that the smaller the value of α-β is, the more iterations are needed for the algorithms to reach convergence.We can see that, for the sameα-β, HOFBP has the highest recovery percentage. For efficiency, HOFBP consumes more time than AFBOMP, however, much less time than IFBP.


## 6.2. Robustness against Sparsity

Figure2 shows the perfect recovery percentages and average run times of all algorithms against the sparsity K changing from 20 to 38 in the noiseless case. We fix other parameters as L=6, N=200, and M=40. HMP performs better than SOMP at the beginning. With the increase of sparsity, HMP begins to perform worse than SOMP at the sparsity 24. SSP performs nearly consistently worse than SOMP, which is consistent with our previous analysis. HOFBP performs better than the other algorithms. IFBP has similar performance with HOFBP; however, it consumes much more time than HOFBP. For efficiency, HOFBP consumes similar time to HMP and AFBOMP but performs better than them.Figure 2
Perfect recovery percentages and average run times of different algorithms againstK with M=40, N=200, L=6.
(a) (b)


## 6.3. Robustness against the Measurement Number

Figure3 shows the perfect recovery percentages and average run times of all algorithms against the measurement number M changing from 34 to 60 in the noiseless case. We fix other parameters as L=6, N=200, and K=30. With the increase of the measurement number, the perfect recovery percentages of all algorithms increase. HMP and SSP perform consistently worse than SOMP. HOFBP has the best performance among all algorithms. IFBP performs similarly to HOFBP; however, it consumes much more time. For efficiency, HOFBP is similar to AFBOMP and HMP, but with better performance. Moreover, HOFBP needs the least measurement number to reach perfect recovery.Figure 3
Perfect recovery percentages and average run times of different algorithms againstM with K=30, N=200, L=6.
(a) (b)


## 6.4. Robustness against the Number of Signals

Figure4 shows the perfect recovery percentages and average run times of all algorithms against the number of signals L changing from 1 to 25 in the noiseless case. We fix other parameters as M=45, N=200, K=30. With the increase of L, the perfect recovery percentages of all algorithms increase. This is because the more measurement vectors are, the more correlations the algorithms can utilize. We can see that HMP performs worse than SOMP when L<11. This illustrates that HMP has less strengths than SOMP in the case of relative less signals. HOFBP performs better than other algorithms. Moreover, it needs less number of signals to reach perfect recovery. For efficiency, we have similar results to those of Sections 4.2 and 4.3.Figure 4
Perfect recovery percentages and average run times of different algorithms againstL with M=45, K=30, N=200.
(a) (b)


## 6.5. Recovery from Noisy Observation

In this part, we test algorithms’ performances against noise. We assume the noise is i.i.d Gaussian noise. The signal-to-noise ratio (SNR) is used to measure the relation between the measurement signal and the noise, which is defined as(8)SNR=10log10YF2NoiseF2.What should be pointed out is that, in the noisy case, the terminate parameters of SOMP, HOFBP, IFBP, and AFBOMP are set asε=10-SNR/20. The terminate parameters of SSP and HMP are still set as the real sparsity.Figure5 shows the perfect recovery percentages and average run times of all algorithms against the change of noise. Other parameters are fixed as M=45, K=30, N=200, and L=6. With the increase of SNR, the perfect recovery percentages of all algorithms increase. SSP and HMP perform consistently worse than SOMP. AFBOMP performs worst, which means it is not suitable for the noisy situation. HOFBP performs consistently better than other algorithms. For efficiency, we have similar result to those of Sections 4.2 and 4.3.Figure 5
Perfect recovery percentages and average run times of different algorithms againstSNR with M=45, K=30, N=200, and L=6.
(a) (b)


## 7. Conclusion

The strengths and drawbacks of previous GP algorithms are analyzed. On this basis, we combine SOMP and FBP to propose the HOFBP algorithm. The choices of the parameters are critical to the performance of this algorithm. However, the simulation results show that, by using proper parameters, HOFBP has better performance than other GP algorithms at the expense of more computing time in some cases. Moreover, HOFBP does not need sparsity as a prior information. It is important to give the exact expressions for the mean square error in the construction using our algorithm. This will be researched in our future work. Moreover, in my opinion, AFBOMP, IFBP, and HOFBP can be combined to improve the performance further. This also can be researched in the future work.



---

*Source: 5965020-2018-07-17.xml*