# Tree-Based Backtracking Orthogonal Matching Pursuit for Sparse Signal Reconstruction

Authors: Yigang Cen; Fangfei Wang; Ruizhen Zhao; Lihong Cui; Lihui Cen; Zhenjiang Miao; Yanming Cen
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Compressed sensing (CS) is a theory which exploits the sparsity characteristic of the original signal in signal sampling and coding. By solving an optimization problem, the original sparse signal can be reconstructed accurately. In this paper, a new Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) algorithm is presented with the idea of the tree model in wavelet domain. The algorithm can convert the wavelet tree structure to the corresponding relations of candidate atoms without any prior information of signal sparsity. Thus, the atom selection process will be more structural and the search space can be narrowed. Moreover, according to the backtracking process, the previous chosen atoms’ reliability can be detected and the unreliable atoms can be deleted at each iteration, which leads to an accurate reconstruction of the signal ultimately. Compared with other compressed sensing algorithms, simulation results show the proposed algorithm’s superior performance to that of several other OMP-type algorithms.


---

## Body

## 1. Introduction

Compressive sensing (CS) [1, 2] aims to recover sparse or compressible signal with low amount of information and high probability. It breaks the traditional rule of Nyquist sampling theorem, which states that a signal’s information is preserved if it is uniformly sampled at a rate at least two times faster than its Fourier bandwidth. By this state-of-the-art signal compression and processing theory, the signal sampling frequency, the cost of processing time, data storage, and transmission can be greatly reduced.For a given orthogonal basisΨ={ψ1,…,ψN}, the signal x∈RN×1 can be represented in terms of the coefficient vector α as
(1)x=∑k=1Nψkαk=Ψα.The corresponding inverse transformation isα=ΨHx, ΨΨH=ΨHΨ=I, and Ψ∈CN×N. Here, I is the identity matrix. We say that x is K-sparse under the orthogonal basis Ψ if only K≪N coefficients αk of x are nonzero.Usually, the signal is not sparse but its coefficient can be considered to be sparse or compressible after some transformations, such as the wavelet transformation.Suppose that a matrixΦ represents the M×N measurement matrix. Then α is accomplished by collecting a measurement vector y of dimension M with M≪N. y can be expressed as y=Φα. Then, (1) becomes
(2)y=Φα=ΦΨHx.Φ is called as the CS measurement matrix and its columns are called atoms. The matrix Φ is rank deficient and hence loses information in general. The CS reconstruction problem wishes to recover the coefficient vector α from the set of M linear measurements y. Since M<N, the reconstruction of α from y is generally ill-posed.The two major algorithmic approaches to sparse recovery are methods based on (l1) minimization and iterative methods (matching pursuits). We now briefly describe these methods, as follows.


### 1.1.  (l1) Minimization

The sparse recovery of this approach can be stated as the problem of finding the sparsest signalα=ΨHx with the given measurements y:
(3)(l0):min∥ΨHx∥l0s.t.y=ΦΨHx.Donoho and his associates advocated the principle that for some measurement matricesΦ, the highly nonconvex combinatorial optimization problem (l0) should be equivalent to its convex relaxation:
(4)(l1):min∥ΨHx∥l1s.t.y=ΦΨHx.
Reference [3] showed that if the measurement matrix satisfies the restricted isometry property (RIP), then a K-sparse signal can be recovered exactly; that is,
(5)(1-δK)∥x∥22≤∥Φx∥22≤(1+δK)∥x∥22.δK is called as the Restricted Isometry Constant of Φ. It has been shown that (l1) minimization can recover a sparse signal exactly under various conditions on restricted isometry constants, see [4, 5]. Then, the convex problem (l1) can be solved using method of convex and even linear programming.


### 1.2. Orthogonal Matching Pursuit (OMP)

An alternative approach to sparse recovery is via iterative algorithms, which find the support of theK-sparse signal α progressively. Once S=supp(α) is found correctly, it is easy to compute the signal α from its measurements y as α=(ΦS)-1y, where ΦS denotes the measurement matrix Φ restricted to columns indexed by S.A basic iterative algorithm is Orthogonal Matching Pursuit (OMP) [6]. OMP recovers the support of α, one index at a time, in n steps. Under a hypothetical assumption that  is an isometry, that is, the columns of Φ are orthogonal, the signal α can be exactly recovered from its measurements y as α=Φ*y.The problem is that theM×N matrix Φ is never an isometry in the interesting range where the number of measurements M is smaller than the ambient dimension N. Even though the matrix is not an isometry, one can still use the notion of coherence in recovery of sparse signals. In that setting, greedy algorithms are used with incoherent dictionaries to recover such signals, see [7, 8]. In our setting, for the commonly used random matrices, one expects the columns to be approximately orthogonal, and the observation vector α=Φ*y to be a good approximation to the original signal α.Tropp and Gilbert [6] analyzed the performance of OMP for Gaussian measurement matrices Φ; a similar result holds for general sub-gaussian matrices. They proved that, for every fixed K-sparse N-dimensional signal α and a random Gaussian measurement matrix Φ, OMP recovers (the support of) α from the measurements y correctly with high probability, provided the number of measurements is M~KlogN.The (l1)-minimization method has the strongest known guarantees of sparse recovery. Once the measurement matrix Φ satisfies the Restricted Isometry Condition, this method works correctly for all sparse signals α. (l1)-minimization is based on linear programming, which has its advantages and disadvantages. One thinks of linear programming as a black box and any development of fast solvers will reduce the running time of the sparse recovery method. On the other hand, it is not very clear what this running time is, as there is no strongly polynomial time algorithm in linear programming yet. All known solvers take time polynomial not only in the dimension of the program N but also on certain condition numbers of the program. While for some classes of random matrices the expected running time of linear programming solvers can be bounded, estimating condition numbers is hard for specific matrices. For example, there is no result yet showing that the Restricted Isometry Condition implies that the condition numbers of the corresponding linear program is polynomial in N.OMP is quite fast, both theoretically and experimentally. It makesn iterations, where each iteration amounts to a multiplication by a N×M matrix Φ* (computing the observation vector α) and solving a least squares problem in dimensions at most M×n. This yields strongly polynomial running time. In practice, OMP is observed to perform faster and is easier to implement than (l1)-minimization. For more details, see [6]. OMP is quite transparent; at each iteration, it selects a new coordinate from the support of the signal α in a very specific and natural way. In contrast, the known (l1)-minimization solvers, such as the simplex method and interior point methods, compute a path toward the solution. However, the geometry of (l1) is clear, whereas the analysis of greedy algorithms can be difficult simply because they are iterative.On the other hand, OMP has weaker guarantees of exact recovery. Unlike (l1)-minimization, the guarantees of OMP are nonuniform: for each fixed sparse signal α and not for all signals, the algorithm performs correctly with high probability. Rauhut has shown that uniform guarantees for OMP are impossible for natural random measurement matrices [9].Moreover, OMP’s condition on measurement matrices given in [6] is more restrictive than the Restricted Isometry Condition. In particular, it is not known whether OMP succeeds in the important class of partial Fourier measurement matrices.These open problems about OMP, first stated in [6] and often reverberated in the Compressed Sensing community, motivated the recent works on the modified OMP algorithms, such as the model-based Compressive Sensing [10], Tree-Based Orthogonal Matching Pursuit [11], Compressive Sampling Matching Pursuit (CoSaMP) [12], Regularized Orthogonal Matching Pursuit (ROMP) [13], and Backtracking-Based Matching Pursuit (BAOMP) [14]. ROMP and CoSaMP require the sparsity level as an input parameter. However, in the most practical applications, this information may not be known before reconstruction. Although the sparsity level is not required for the OMP and BAOMP, they do not use the characteristics of the sparse representation, such as the tree structure of wavelet transform. In this paper, a new Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) algorithm is presented based on the tree model in wavelet domain. Our algorithm converts the wavelet tree structure to the corresponding relations of candidate atoms without the prior information of signal sparsity level. Also, combing with the backtracking algorithm, the unreliable atoms can be deleted. Compared with OMP, ROMP, and BAOMP algorithms, the atom selection process will be more traceable, normalizable, and structural.


## 1.1.  (l1) Minimization

The sparse recovery of this approach can be stated as the problem of finding the sparsest signalα=ΨHx with the given measurements y:
(3)(l0):min∥ΨHx∥l0s.t.y=ΦΨHx.Donoho and his associates advocated the principle that for some measurement matricesΦ, the highly nonconvex combinatorial optimization problem (l0) should be equivalent to its convex relaxation:
(4)(l1):min∥ΨHx∥l1s.t.y=ΦΨHx.
Reference [3] showed that if the measurement matrix satisfies the restricted isometry property (RIP), then a K-sparse signal can be recovered exactly; that is,
(5)(1-δK)∥x∥22≤∥Φx∥22≤(1+δK)∥x∥22.δK is called as the Restricted Isometry Constant of Φ. It has been shown that (l1) minimization can recover a sparse signal exactly under various conditions on restricted isometry constants, see [4, 5]. Then, the convex problem (l1) can be solved using method of convex and even linear programming.


## 1.2. Orthogonal Matching Pursuit (OMP)

An alternative approach to sparse recovery is via iterative algorithms, which find the support of theK-sparse signal α progressively. Once S=supp(α) is found correctly, it is easy to compute the signal α from its measurements y as α=(ΦS)-1y, where ΦS denotes the measurement matrix Φ restricted to columns indexed by S.A basic iterative algorithm is Orthogonal Matching Pursuit (OMP) [6]. OMP recovers the support of α, one index at a time, in n steps. Under a hypothetical assumption that  is an isometry, that is, the columns of Φ are orthogonal, the signal α can be exactly recovered from its measurements y as α=Φ*y.The problem is that theM×N matrix Φ is never an isometry in the interesting range where the number of measurements M is smaller than the ambient dimension N. Even though the matrix is not an isometry, one can still use the notion of coherence in recovery of sparse signals. In that setting, greedy algorithms are used with incoherent dictionaries to recover such signals, see [7, 8]. In our setting, for the commonly used random matrices, one expects the columns to be approximately orthogonal, and the observation vector α=Φ*y to be a good approximation to the original signal α.Tropp and Gilbert [6] analyzed the performance of OMP for Gaussian measurement matrices Φ; a similar result holds for general sub-gaussian matrices. They proved that, for every fixed K-sparse N-dimensional signal α and a random Gaussian measurement matrix Φ, OMP recovers (the support of) α from the measurements y correctly with high probability, provided the number of measurements is M~KlogN.The (l1)-minimization method has the strongest known guarantees of sparse recovery. Once the measurement matrix Φ satisfies the Restricted Isometry Condition, this method works correctly for all sparse signals α. (l1)-minimization is based on linear programming, which has its advantages and disadvantages. One thinks of linear programming as a black box and any development of fast solvers will reduce the running time of the sparse recovery method. On the other hand, it is not very clear what this running time is, as there is no strongly polynomial time algorithm in linear programming yet. All known solvers take time polynomial not only in the dimension of the program N but also on certain condition numbers of the program. While for some classes of random matrices the expected running time of linear programming solvers can be bounded, estimating condition numbers is hard for specific matrices. For example, there is no result yet showing that the Restricted Isometry Condition implies that the condition numbers of the corresponding linear program is polynomial in N.OMP is quite fast, both theoretically and experimentally. It makesn iterations, where each iteration amounts to a multiplication by a N×M matrix Φ* (computing the observation vector α) and solving a least squares problem in dimensions at most M×n. This yields strongly polynomial running time. In practice, OMP is observed to perform faster and is easier to implement than (l1)-minimization. For more details, see [6]. OMP is quite transparent; at each iteration, it selects a new coordinate from the support of the signal α in a very specific and natural way. In contrast, the known (l1)-minimization solvers, such as the simplex method and interior point methods, compute a path toward the solution. However, the geometry of (l1) is clear, whereas the analysis of greedy algorithms can be difficult simply because they are iterative.On the other hand, OMP has weaker guarantees of exact recovery. Unlike (l1)-minimization, the guarantees of OMP are nonuniform: for each fixed sparse signal α and not for all signals, the algorithm performs correctly with high probability. Rauhut has shown that uniform guarantees for OMP are impossible for natural random measurement matrices [9].Moreover, OMP’s condition on measurement matrices given in [6] is more restrictive than the Restricted Isometry Condition. In particular, it is not known whether OMP succeeds in the important class of partial Fourier measurement matrices.These open problems about OMP, first stated in [6] and often reverberated in the Compressed Sensing community, motivated the recent works on the modified OMP algorithms, such as the model-based Compressive Sensing [10], Tree-Based Orthogonal Matching Pursuit [11], Compressive Sampling Matching Pursuit (CoSaMP) [12], Regularized Orthogonal Matching Pursuit (ROMP) [13], and Backtracking-Based Matching Pursuit (BAOMP) [14]. ROMP and CoSaMP require the sparsity level as an input parameter. However, in the most practical applications, this information may not be known before reconstruction. Although the sparsity level is not required for the OMP and BAOMP, they do not use the characteristics of the sparse representation, such as the tree structure of wavelet transform. In this paper, a new Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) algorithm is presented based on the tree model in wavelet domain. Our algorithm converts the wavelet tree structure to the corresponding relations of candidate atoms without the prior information of signal sparsity level. Also, combing with the backtracking algorithm, the unreliable atoms can be deleted. Compared with OMP, ROMP, and BAOMP algorithms, the atom selection process will be more traceable, normalizable, and structural.


## 2. Tree-Based Backtracking Orthogonal Matching Pursuit (TBOMP) Algorithm

In this section, we will first review the wavelet tree structure. Second, the proposed TBOMP algorithm will be presented in detail.


### 2.1. Wavelet Tree Structure

Consider a signalx of length N=2L, after L-level wavelet transformations, the set of K-tree sparse signals is defined as
(6)Γk={x=υLν+∑i=L1∑j=12L-iωi,jψi,j:ω|ΩC=0,|Ω|=K},
where ν is the scaling function and ψi,j is the wavelet function at scale i and offset j. The wavelet transform consists of the scale coefficient υL and wavelet coefficients ωi,j at scale i, 1≤i≤L, and position j, with 1≤j≤2L-i.Suppose thatα=[υL,ωL,0,ωL-1,0,ωL-1,1,ωL-2,0,…]T is the vector of the scaling and wavelet coefficients of x with the maximum decomposition level L. Also, it is a set of wavelet coefficients Ω forms a connected subtree [10]. The set Ω defines a subspace of signals whose support is contained in Ω, which means that all wavelet coefficients outside Ω are approximately zero. The nested structure of wavelet coefficients creates a parent/child relationship between wavelet coefficients at different scales. We say that ωi+1,⌊j/2⌋ (⌊·⌋ denotes rounded down) is the parent of ωi,j. Also, ωi-1,2j and ωi-1,2j+1 are the children of ωi,j. These relations can be expressed graphically by the wavelet coefficient tree in Figure 2(a). The relationship between the parent and child nodes is that the index value of the parent node in a level is 1/2 times the index of the child node.A kind of tree structure (greedy tree) was proposed in [15]. For the greedy tree, if a coefficient is significant then its child and all of its grandchildren are likely significant [11]. Figure 1 depicts two cases of greedy tree approximation. The number of each node is the wavelet coefficient modulus. Nodes not labeled depict zeros. In the first case, the wavelet coefficients decay monotonically along the tree branches toward the leaves. Suppose that the wavelet tree Ω containing P wavelet coefficients; that is, |Ω|=P. The P-term greedy tree approximation (here, we assume that P=4) can be proceeded in three steps: (1) find the p, p≤4 largest wavelet coefficient terms; (2) form the smallest connected rooted subtree that contains all of these p coefficients; and then (3) increase p until |Ω|=4.Greedy tree search.
(a)
(b)Wavelet tree structure.
(a)
Wavelet tree structure
(b)
Process of tree nodes selection in the TBOMPInitializingp=2, two coefficients 10 and 8 will be found and will form a minimum, connected subtree Ω. Gradually increase p until p=4, the greedy tree approximation forms the connected rooted subtree Ω, 10-8-4-3, with 4 nodes that maximize the sum of the wavelet coefficients in the subtree. This process was shown in Figure 1(a), the error is small. Another case was shown in Figure 1(b), when the wavelet coefficients do not decay monotonically along the tree branches toward the leaves, an isolated significant coefficient away from the root will be selected, either of its all ancestor coefficients. These ancestor coefficients may be very small, which will increase the approximation error. For example, initializing p=2, then two coefficients 10 and 8 will be found and the resulted subtree is 10-0-0-8 with p=|Ω|=4. Obviously, the error is large.We can see that the process of greedy tree approximation is simple, but when the tree includes isolated large coefficients far from the tree root, the approximation error will be increased. Thus, backtracking is imposed to deleting the wrong nodes selected by the greedy tree. This will be illustrated in the Section2.2.


### 2.2. Tree-Based Backtracking Orthogonal Matching Pursuit (TBOMP) Algorithm

Our proposed Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) is as follows.Algorithm 1 (TBOMP).
Symbol Description
ω—wavelet high frequency coefficient vector;
ω
^—reconstruction wavelet high frequency coefficient vector;
A—measurement matrix, y=Aω;
a
i—the ith column vector of A, 1≤i≤N;
μ
1
,
μ
2—parameters of thresholds, μ1,μ2∈[0,1];
Λ
n—index set, Λ denotes the index set of all columns {ai} of matrix A;
n
max—number of maximum iterations allowed;
Γ
n—atom-deleting set in the nth iteration;
C
n—candidate set of the root atoms in the nth iteration;
F
n—family set that consists of the subtrees corresponding to the root nodes in Cn.
Initialization.
r
0
=
y (initial residual), Λ0=∅, Γ0=∅, and C0=∅.
Loop
(1)
Initial selection: select the candidate setCn with absolute values of correlations satisfying:
(7)|〈rn-1,aCn〉|≥μ1·maxi∈Λn|〈rn-1,ai〉|,Λn=Λ∖Λn-1.
(2)
According to the 2-times relationship of wavelet tree node indices, find the wavelet tree rooted at each node inCn. Then the family set Fn consists of the atoms indexed by Cn and all of their families can be found.
For example, assume thatCn={cn1,cn2,…,cnQ}, then the wavelet subtrees rooted at cn1,cn2,…,cnQ will be found, respectively, in this step. The index sets of these Q trees are denoted as Fn1,Fn2,…,FnQ.
(3)
Computeω^Fnqn=(AFnqHAFnq)-1AFnqHy, 1≤q≤Q.
(4)
FindFnq~ such that ωFnq~n minimizing the residual as follows:
(8)q~=argmin1≤q≤Q∥y-AFnqω^Fnqn∥2.
(5)
Select atom deleting index setΓn satisfying
(9)|ω^Λn-1∪Fnq~n|≤μ2·max|ω^Fnq~n|.
(6)
SetΛn={Λn-1∪Fnq~}∖Γn, a{i:i∈Λn}=∅, and update the residual as follows:
(10)rn=y-AΛnω^Λnn.
(7)
If∥rn∥2<ε or if n=nmax, quit the iteration; otherwise, set n=n+1, go to step 1.
End Loop.Output. the estimated support set Λn and the nonzero values ω^Λn=(AΛnHAΛn)-1AΛnHy.As seen in the above algorithm, we combined the characteristics of tree structure and the BAOMP algorithm. In the first step, TBOMP selects candidate setCn whose correlations between the columns of ΦΛn and the residual rn-1 are not smaller than μ1·maxi∈Λn|〈rn-1,ai〉|, Λn=Λ∖Λn-1. Here, the constant μ1 is used to adaptively decide how many atoms are chosen at each time. Then the atoms corresponding to the elements of Cn are set as the root nodes of subtrees. As we mentioned in Section 2.1, due to the 2-times relationship between the indices of parent and child nodes, the subtree of each atom corresponding to an index in Cn can be found to form the family set Fnq, which consists of the indices of the family atoms in the qth subtree. In the third step, least square method is applied to obtain the reconstruction wavelet high frequency coefficients ω^Fnqn corresponding to the atoms indexed by Fnq. Then the optimal subtree indexed by Fnq~ will be selected according to step (4). In this step, there may exist insignificant atoms in aFnq~. This is because that we only simply applied the 2-times relationship discipline in the searching processing of subtrees. Thus, the backtracking deleting method is introduced in the algorithm to identify the true support set of Fnq~. The backtracking deleting set Γn consists of the indices corresponding to all the reconstructed coefficients satisfying (9). Then, the index set is updated by Λn={Λn-1∪Fnq~}∖Γn at this iteration. According to the atoms corresponding to the indices in the set Λn, the reconstruction coefficients ω^Λnn can be computed. Finally, update the residual by (10) and go to the next iteration. If ∥rn∥2<ε or n=nmax, quit the iteration.In the TBOMP, the process of tree nodes selection was shown in Figure2; the first step of the algorithm is to select candidate set Cn by (7). For example, suppose that C1={ωL,0,ωL-2,3} was chosen at the first iteration. The nodes of subtree ② rooted at ωL,0 and the family nodes rooted at ωL-2,3 are the significant coefficients needed to be found. According to the 2-times relationship of wavelet tree node indices and Figure 2(b), ωL-1,0 and ωL-2,0 are the child and grandchild nodes of ωL,0. Thus, subtree ① rooted at the node ωL,0 will be found in the first iteration.Now we assume that the subtree ① is the optimal tree corresponding toωL,0. At the end of this iteration, the backtracking algorithm will remove the node ωL-2,0 according to step 5 of the TBOMP algorithm described above. In the remaining iteration, node ωL-2,1 will be choosen as the child node of ωL-1,0. Ultimately, subtree ② will be found accurately. Analogously, the searching process of the subtree rooted at the node ωL-2,3 is the same, and it can be proceeded simultaneously.These characteristics of tree structure provide a new way for the study of reconstruction algorithm. Thanks to the tree structure of wavelet coefficients, when the signal is sparsely represented by the wavelet transform, it also provides a clew for the selection of atoms in the reconstruction algorithm. This will greatly improve the reliability of the atom selection.The coefficients of wavelet decomposition include low-frequency coefficients and high-frequency coefficients (scaling coefficients and wavelet coefficients inα). The more levels of wavelet decomposition, the less low-frequency coefficients, and more important information is reserved in the high-frequency coefficients. Compared with the high-frequency coefficients, the number of low-frequency coefficients are much less if the decomposition level is big enough. Since the low-frequency coefficients play an important role in the wavelet reconstruction, in our proposed algorithm, only the high-frequency coefficients are measured by measurement matrix. For the reconstruction, we combine the reconstructed high-frequency coefficients ω^ and the unprocessed low-frequency coefficients. Then the inverse wavelet transform is applied to obtain a reconstructed x^ of the original signal x.


## 2.1. Wavelet Tree Structure

Consider a signalx of length N=2L, after L-level wavelet transformations, the set of K-tree sparse signals is defined as
(6)Γk={x=υLν+∑i=L1∑j=12L-iωi,jψi,j:ω|ΩC=0,|Ω|=K},
where ν is the scaling function and ψi,j is the wavelet function at scale i and offset j. The wavelet transform consists of the scale coefficient υL and wavelet coefficients ωi,j at scale i, 1≤i≤L, and position j, with 1≤j≤2L-i.Suppose thatα=[υL,ωL,0,ωL-1,0,ωL-1,1,ωL-2,0,…]T is the vector of the scaling and wavelet coefficients of x with the maximum decomposition level L. Also, it is a set of wavelet coefficients Ω forms a connected subtree [10]. The set Ω defines a subspace of signals whose support is contained in Ω, which means that all wavelet coefficients outside Ω are approximately zero. The nested structure of wavelet coefficients creates a parent/child relationship between wavelet coefficients at different scales. We say that ωi+1,⌊j/2⌋ (⌊·⌋ denotes rounded down) is the parent of ωi,j. Also, ωi-1,2j and ωi-1,2j+1 are the children of ωi,j. These relations can be expressed graphically by the wavelet coefficient tree in Figure 2(a). The relationship between the parent and child nodes is that the index value of the parent node in a level is 1/2 times the index of the child node.A kind of tree structure (greedy tree) was proposed in [15]. For the greedy tree, if a coefficient is significant then its child and all of its grandchildren are likely significant [11]. Figure 1 depicts two cases of greedy tree approximation. The number of each node is the wavelet coefficient modulus. Nodes not labeled depict zeros. In the first case, the wavelet coefficients decay monotonically along the tree branches toward the leaves. Suppose that the wavelet tree Ω containing P wavelet coefficients; that is, |Ω|=P. The P-term greedy tree approximation (here, we assume that P=4) can be proceeded in three steps: (1) find the p, p≤4 largest wavelet coefficient terms; (2) form the smallest connected rooted subtree that contains all of these p coefficients; and then (3) increase p until |Ω|=4.Greedy tree search.
(a)
(b)Wavelet tree structure.
(a)
Wavelet tree structure
(b)
Process of tree nodes selection in the TBOMPInitializingp=2, two coefficients 10 and 8 will be found and will form a minimum, connected subtree Ω. Gradually increase p until p=4, the greedy tree approximation forms the connected rooted subtree Ω, 10-8-4-3, with 4 nodes that maximize the sum of the wavelet coefficients in the subtree. This process was shown in Figure 1(a), the error is small. Another case was shown in Figure 1(b), when the wavelet coefficients do not decay monotonically along the tree branches toward the leaves, an isolated significant coefficient away from the root will be selected, either of its all ancestor coefficients. These ancestor coefficients may be very small, which will increase the approximation error. For example, initializing p=2, then two coefficients 10 and 8 will be found and the resulted subtree is 10-0-0-8 with p=|Ω|=4. Obviously, the error is large.We can see that the process of greedy tree approximation is simple, but when the tree includes isolated large coefficients far from the tree root, the approximation error will be increased. Thus, backtracking is imposed to deleting the wrong nodes selected by the greedy tree. This will be illustrated in the Section2.2.


## 2.2. Tree-Based Backtracking Orthogonal Matching Pursuit (TBOMP) Algorithm

Our proposed Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) is as follows.Algorithm 1 (TBOMP).
Symbol Description
ω—wavelet high frequency coefficient vector;
ω
^—reconstruction wavelet high frequency coefficient vector;
A—measurement matrix, y=Aω;
a
i—the ith column vector of A, 1≤i≤N;
μ
1
,
μ
2—parameters of thresholds, μ1,μ2∈[0,1];
Λ
n—index set, Λ denotes the index set of all columns {ai} of matrix A;
n
max—number of maximum iterations allowed;
Γ
n—atom-deleting set in the nth iteration;
C
n—candidate set of the root atoms in the nth iteration;
F
n—family set that consists of the subtrees corresponding to the root nodes in Cn.
Initialization.
r
0
=
y (initial residual), Λ0=∅, Γ0=∅, and C0=∅.
Loop
(1)
Initial selection: select the candidate setCn with absolute values of correlations satisfying:
(7)|〈rn-1,aCn〉|≥μ1·maxi∈Λn|〈rn-1,ai〉|,Λn=Λ∖Λn-1.
(2)
According to the 2-times relationship of wavelet tree node indices, find the wavelet tree rooted at each node inCn. Then the family set Fn consists of the atoms indexed by Cn and all of their families can be found.
For example, assume thatCn={cn1,cn2,…,cnQ}, then the wavelet subtrees rooted at cn1,cn2,…,cnQ will be found, respectively, in this step. The index sets of these Q trees are denoted as Fn1,Fn2,…,FnQ.
(3)
Computeω^Fnqn=(AFnqHAFnq)-1AFnqHy, 1≤q≤Q.
(4)
FindFnq~ such that ωFnq~n minimizing the residual as follows:
(8)q~=argmin1≤q≤Q∥y-AFnqω^Fnqn∥2.
(5)
Select atom deleting index setΓn satisfying
(9)|ω^Λn-1∪Fnq~n|≤μ2·max|ω^Fnq~n|.
(6)
SetΛn={Λn-1∪Fnq~}∖Γn, a{i:i∈Λn}=∅, and update the residual as follows:
(10)rn=y-AΛnω^Λnn.
(7)
If∥rn∥2<ε or if n=nmax, quit the iteration; otherwise, set n=n+1, go to step 1.
End Loop.Output. the estimated support set Λn and the nonzero values ω^Λn=(AΛnHAΛn)-1AΛnHy.As seen in the above algorithm, we combined the characteristics of tree structure and the BAOMP algorithm. In the first step, TBOMP selects candidate setCn whose correlations between the columns of ΦΛn and the residual rn-1 are not smaller than μ1·maxi∈Λn|〈rn-1,ai〉|, Λn=Λ∖Λn-1. Here, the constant μ1 is used to adaptively decide how many atoms are chosen at each time. Then the atoms corresponding to the elements of Cn are set as the root nodes of subtrees. As we mentioned in Section 2.1, due to the 2-times relationship between the indices of parent and child nodes, the subtree of each atom corresponding to an index in Cn can be found to form the family set Fnq, which consists of the indices of the family atoms in the qth subtree. In the third step, least square method is applied to obtain the reconstruction wavelet high frequency coefficients ω^Fnqn corresponding to the atoms indexed by Fnq. Then the optimal subtree indexed by Fnq~ will be selected according to step (4). In this step, there may exist insignificant atoms in aFnq~. This is because that we only simply applied the 2-times relationship discipline in the searching processing of subtrees. Thus, the backtracking deleting method is introduced in the algorithm to identify the true support set of Fnq~. The backtracking deleting set Γn consists of the indices corresponding to all the reconstructed coefficients satisfying (9). Then, the index set is updated by Λn={Λn-1∪Fnq~}∖Γn at this iteration. According to the atoms corresponding to the indices in the set Λn, the reconstruction coefficients ω^Λnn can be computed. Finally, update the residual by (10) and go to the next iteration. If ∥rn∥2<ε or n=nmax, quit the iteration.In the TBOMP, the process of tree nodes selection was shown in Figure2; the first step of the algorithm is to select candidate set Cn by (7). For example, suppose that C1={ωL,0,ωL-2,3} was chosen at the first iteration. The nodes of subtree ② rooted at ωL,0 and the family nodes rooted at ωL-2,3 are the significant coefficients needed to be found. According to the 2-times relationship of wavelet tree node indices and Figure 2(b), ωL-1,0 and ωL-2,0 are the child and grandchild nodes of ωL,0. Thus, subtree ① rooted at the node ωL,0 will be found in the first iteration.Now we assume that the subtree ① is the optimal tree corresponding toωL,0. At the end of this iteration, the backtracking algorithm will remove the node ωL-2,0 according to step 5 of the TBOMP algorithm described above. In the remaining iteration, node ωL-2,1 will be choosen as the child node of ωL-1,0. Ultimately, subtree ② will be found accurately. Analogously, the searching process of the subtree rooted at the node ωL-2,3 is the same, and it can be proceeded simultaneously.These characteristics of tree structure provide a new way for the study of reconstruction algorithm. Thanks to the tree structure of wavelet coefficients, when the signal is sparsely represented by the wavelet transform, it also provides a clew for the selection of atoms in the reconstruction algorithm. This will greatly improve the reliability of the atom selection.The coefficients of wavelet decomposition include low-frequency coefficients and high-frequency coefficients (scaling coefficients and wavelet coefficients inα). The more levels of wavelet decomposition, the less low-frequency coefficients, and more important information is reserved in the high-frequency coefficients. Compared with the high-frequency coefficients, the number of low-frequency coefficients are much less if the decomposition level is big enough. Since the low-frequency coefficients play an important role in the wavelet reconstruction, in our proposed algorithm, only the high-frequency coefficients are measured by measurement matrix. For the reconstruction, we combine the reconstructed high-frequency coefficients ω^ and the unprocessed low-frequency coefficients. Then the inverse wavelet transform is applied to obtain a reconstructed x^ of the original signal x.


## 3. Simulation Results

In this section, several experiments will be given for the TBOMP algorithm. In the first experiment, the original signalx is a one-dimensional blocks signal with length N=256. It was recovered from M=64 measurements by using the Gaussian random measurement matrix. The wavelet decomposition level is 4 and the wavelet function is Db1. Figure 3 shows the reconstruction result of 7th iterations by using the TBOMP algorithm.Reconstruction signal by TBOMP algorithm.
(a)
Original signal
(b)
Wavelet coefficients of 4-level decompositions. The coefficients are arrayed as{a4,d4,d3,d2,d1}
(c)
Wavelet coefficients recovery of the TBOMP algorithmIn the first iteration of the TBOMP algorithm, according to the parent-child relations of wavelet tree, some unreliable atoms will be chosen, which leads to a wrong reconstruction result. As marked by the cycles in Figure4(a). Then according to the backtracking deleting method, the wrong selected atoms can be deleted. After the second and the third iterations, some atoms are still not found. After the 7th iteration, the reconstruction result (Figure 3(c)) with TBOMP algorithm is exactly same as the original wavelet coefficients shown in Figure 3(b).Reconstruction of TBOMP algorithm.
(a)
Reconstructed wavelet coefficients after the first selection of the wavelet nodes
(b)
Backtracking deleting of the wrong nodes after the first iteration of TBOMP
(c)
Reconstructed wavelet coefficients after the third selection of the wavelet nodes
(d)
Backtracking deleting of the wrong nodes after the third iteration of TBOMPSimilar results can be obtained for other signals. Reconstruction results of Doppler and Heavysine signals by using our TBOMP algorithm are shown in Figure5. Here, we compared our reconstruction results with the classical OMP algorithm, M/N=1/4.Reconstruction results of Doppler and HeavySine signal by TBOMP and OMP algorithms.
(a)
Original and Reconstruction signals of Doppler signal by TBOMP
(b)
A zoom-in view of Figure5(a)
(c)
Original and Reconstruction signals of Doppler signal by OMP
(d)
A zoom-in view of Figure5(c)
(e)
Original and Reconstruction signals of Heavysine signal by TBOMP
(f)
A zoom-in view of Figure5(e)
(g)
Original and Reconstruction signals of Heavysine signal by OMP
(h)
A zoom-in view of Figure5(g)In the next experiment, we will compare the TBOMP with some popular algorithms such as OMP, ROMP, and BAOMP. Here, only the high frequency coefficients are measured; the low-frequency coefficients will not be processed [16]. The wavelet function is choosen as the “sym8” in MATLAB. The decomposition level is 5 for these four algorithms. Define SNR=20log10(std(x)/std(x^-x)), where std denotes the standard deviation. Because of the randomness of the sensing matrix, numerical result at each time is different. Hereafter, we use the same sensing matrix in one experiment for these four algorithms.We use the Bumps signal of lengthN=2048 and change the values of M simultaneously in order to guarantee the same experiment condition. After 5 layers of wavelet decompositions, there are 64 low-pass coefficients in the 5th decomposition layer and total 1984 high-pass coefficients in the 5 decomposition layers. In order to obtain a fair comparison, in the Figure 6, the measurements number used in these four algorithms is 500-64=436. For sake of simplicity, when we mention that M measurements in the TBOMP, we means that M is the sum of the low-pass coefficient number and the measurement number of the high-pass coefficients.Comparison signal of TBOMP and BAOMP in time domain.
(a)
TBOMP algorithm,M=500, SNR=35.3231 (dB)
(b)
BAOMP algorithm,M=500, SNR=33.5404 (dB)WhenM=500, the compression ratio is about 1⁄4. The reconstruction results of Bumps signal of TBOMP and BAOMP are shown Figure 6. The SNR of TBOMP is about 1.8 dB higher than the BAOMP.Since ROMP requires the sparsity levelK to be known for exact recovery, in the experiments, the best sparsity value K of the wavelet coefficients can be estimated according to repeated experiments and then used in the simulations. Figure 7 shows the SNR comparison results for different values of M. The values of M are selected as 200, 500, 800, 1100, and 1400, respectively. For each M, we conduct the experiment 10 independent trials and calculate the average SNR. It is obviously that the reconstruction result of TBOMP algorithm is superior to others.Figure 7
SNR comparison for different values ofM.


## 4. Conclusion

Sparse reconstruction algorithm is one of the three core problems (signal sparse representation, measurement matrix design, and reconstruction algorithm design) of CS. The existed sparse reconstruction algorithms such as ROMP and CoSaMP algorithms employ the sparsityK as the prior knowledge for exact recovery, which has many limitations for the realistic applications. However, although the sparsity level are not required for OMP and BAOMP algorithms, they do not use the characteristics of special sparse basis to improve the performance of the algorithms. In this paper, a new Tree-based Backtracking Orthogonal Matching Pursuit (TBOMP) algorithm was proposed based on the tree model in wavelet domain. Our algorithm can convert the wavelet tree structures to the corresponding relations of candidate atoms without any prior information of signal sparsity level. Moreover, the unreliable atoms can be deleted according to the backtracking algorithm. Compared with other compressive sensing algorithms (OMP, ROMP, and BAOMP), the signal reconstruction results of TBOMP outperform the above mentioned CS algorithms.



---

*Source: 864132-2013-11-10.xml*