Sparse Signal Recovery Using gOMP Assisted ... - Semantic Scholar
1
Sparse Signal Recovery Using gOMP Assisted mOLS
arXiv:1511.08575v1 [cs.IT] 27 Nov 2015
Samrat Mukhopadhyay1, Siddhartha Satpathi2 and Mrityunjoy Chakraborty3 Department of Electronics and Electrical Communication Engineering, Indian Institute of Technology, Kharagpur, INDIA 1 Email: [email protected],[email protected], 3 [email protected] Abstract—Because of fast convergence in finite number of steps and low computational complexity, signal recovery from compressed measurements using greedy algorithms have generated a large amount of interest in recent years. Among these greedy algorithms OMP is well studied and recently its generalization, gOMP, have also drawn attention. On the other hand OLS and its generalization mOLS have been studied in recent years because of their potential advantages in signal recovery guarantee compared to OMP or gOMP. But OLS and mOLS have the shortcomings of high computational complexity. In this paper we propose a new algorithm which uses gOMP to preselect a N length set out of which mOLS selects its L best coefficients. We have shown that our new algorithm, named gOMPamOLS, is guaranteed to reconstruct a K sparse signals perfectly from compressed measurements if the restricted √ L isometry constant satisfies δLK+N < √ . Moreover √ L+ K +L experimental results indicate that gOMPamOLS is robust under sensing matrix with correlated dictionary in terms of signal reconstruction and is much faster than OLS. Index Terms—Compressive Sensing, gOMP, mOLS, restricted isometry property
I. I NTRODUCTION There has been considerable interest for some time to recover signals from their partial measurements. Mathematically the problem can be expressed as finding solution for the following set of linear equations y = Φx where Φ ∈ Rm×n and m < n. Since this is an undetermined set of linear equations, the number of solutions is infinite, without further information about structure of x. A natural problem that then follows is to find solutions for these set of equations which has the lowest sparsity i.e. lowest number of nonzero elements, i.e the constarined l0 norm minimization problem min kxk0
x∈Rn
s.t. y = Φx However, this problem requires combinatorial search which is computationally intractable. The way out to this problem was found by Candes and Tao [1] in a seminal paper where it was shown that the l0 norm can be solved exactly using l1 norm optimization if the number of measurements is O(K ln n) where K is the sparsity of the vector x. The problem that
is required to solve under this condition is min kxk1
x∈Rn
s.t. y = Φx This is called the basis pursuit(BP) [2]. This can be solved by standard Linear programming approach. A. Greedy approaches Linear programming approach has cubic complexity on the number of columns n, which can be too slow for large n (∼ 1000 is common in many applications). In that case greedy approaches play the role of useful alternatives. A greedy algorithm utilizes some greedy rule to reconstruct the vector at each step of the algorithm. They are often faster than LP. Orthogonal Matching Pursuit(OMP) [3][4], Orthogonal Least Squares(OLS) [5][6] are among the earliest of these greedy algorithms. These algorithms greedily select indices at each step of the algorithm and append them to the already constructed support create a succesively increasing support over which they take projections to get reconstructed signal. The only difference between these two algorithms is in the method of selecting an index in the identification step. At an iteration k OMP chooses an index i which maximizes the absolute
correlation rk−1 , φi where rk−1 is the residual error vector at step k − 1 and φi is a column of the measurement matrix Φ. In OLS the identification is performed by choosing an index i which would minimize the resulting residual vector norm. Wang et al introduced gOMP algorithm [7] and mOLS algorithm [8] as generalizations of OMP and OLS algorithms respectively. In gOMP and mOLS in the identification step, the algorithms choose multiple instead of single index. See Table. I for details. The motivation for introducing both gOMP as well as mOLS is that the time complexity of the algorithms can be markedly improved if at each iteration the algorithms choose multiple indices instead of just one. This increases the probability of choosing a correct index at an iteration and hence, for a given sparse signal with sparsity K, the algorithms converge faster than OMP or OLS. Wang et al [7] has shown that the computational complexity for OMP is about 2Kmn + 3K 2 m floating point operations per second(flops) where m × n is the size of the measurement matrix and the computational complexity of gOMP is about 2smn + (2N 2 + N )s2 m flops where N is the number of indices that gOMP chooses at a time at each iteration and s is the number of iterations required for the algorithm to stop. By extensive simulation they have shown that gOMP indeed takes a small number of
2
TABLE I: gOMP
AND
mOLS A LGORITHM
(a) gOMP A LGORITHM
(b) mOLS A LGORITHM
Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅ While (krk k2 ≥ ǫ and k < K) k = k+1 X k−1 2 Identify: hk = arg max kφT k2 i r S:|S|=N i∈S
Augment: T k = T k−1 ∪ hk Estimate: xk = arg min Update: rk =
Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅ While (krk k2 ≥ ǫ and k < K) k =k+1 X Identify: hk = arg min kP⊥ yk22 T k−1 ∪{i} S:|S|=L i∈S
Augment: T k = T k−1 ∪ hk Estimate: xk = arg min
ky − Φuk2
u:u∈Rn , supp(u)⊂T k y − Φ T k xk
End While
Update: rk =
ky − Φuk2
u:u∈Rn , supp(u)⊂T k y − Φ T k xk
End While
Output: estimated support set Tˆ = arg min kxk − xkS k2 and K-sparse S:|S|=K
Output: estimated support set Tˆ = arg min kxk − xkS k2 and K-sparse S:|S|=K
ˆ satisfying x ˆ Tˆ = xkˆ , x ˆ H\Tˆ = 0 signal x
ˆ satisfying x ˆ Tˆ = xkˆ , x ˆ H\Tˆ = 0 signal x
T
simulation s to recover the K sparse signal which effectively reduces its computational complexity. Similar conclusions are drawn for mOLS compared to OLS by Wang et al [8] through extensive simulations. Other algorithms of the same flavour as gOMP or mOLS are stagewise OMP(StOMP) [9], regularized OMP(ROMP) [10], Compressed Sampling Matching Pursuit(CoSaMP) [11] and Subspace Pursuit (SP) [12]. OMP has been the representative of these algorithms because of its simplicity and recovery performance. Tropp and Gilbert [4] have shown that OMP can recover the original sparse vector from a few measurements with exceedingly high probability when the measurement matrix has entries i.i.d Gaussian. But recently Soussen et al [13] have shown that OLS performs considerably better than OMP in recovering sparse signals from few measurements when the measurement matrix has correlated entries. They have also shown that there exists measurement dictionaries for which OMP fails to recover the original signal whereas OLS is unaffected. B. Our contributions in this paper • In this paper we propose a new algorithm, referred to as gOMPamOLS (short for gOMP assisted mOLS), that exploits the capability of gOMP to select “good” indices at a step as well as the robustness of mOLS in maintaining high recovery performance in the presence measurement matrices with correlated dictionaries. • We find a criteria under which our proposed algorithm can exactly recover the unknown K-sparse signal which is based upon the restricted isometry property of the measurement matrix. A matrix Φ ∈ Rm×n is said to have restricted isometry property (RIP) of order K if ∃ δ > 0 such that ∀x ∈ Rn with sparsity K, the matrix Φ satisfies the following inequalities (1 − δ)kxk22 ≤ kΦxk22 ≤ (1 + δ)kxk22 The smallest δ for which this property is satisfied is called the restricted isometry constant of order K of the matrix Φ and is denoted by δK . We show that gOMPamOLS with gOMP selecting N indices and mOLS selecting
T
L ≤ N indices (See Table. II for details of the algorithm) is guaranteed to recover a K sparse vector exactly from compressed measurements if the matrix Φ obeys RIP with constant √ L √ δLK+N < √ K+L+ L •
We compare simulation results for recovery performances of OMP, OLS and gOMPamOLS and show that gOMPamOLS performs much better than OMP in recovering sparse vector for measurement matrices with correlated dictionary. Also, gOMPmaOLS is much faster than OLS and thus becomes a suitable upgradation of both OMP and OLS. II. P RELIMINARIES
A. Notations The following notations will be used throughout the paper. We will denote by k · kp the lp norm of a vector v, i.e., Pn 1/p kvkp = ( i=1 |vi |p ) . Φ ∈ Rm×n will denote the measurement matrix with m rows and n columns with m < n. The i th column of Φ will be denoted by φi ; i = 1, 2, · · · , n. We assume that all the columns of Φ have unit l2 norm. K will denote the sparsity level estimated beforehand. xS will denote the vector x restricted to a subset S. H will denote the set of all the indices {1, 2, · · · , n}. T will refer to the support set of the original vector x. Similarly, ΦS will denote the submatrix of Φ formed with the columns restricted to index set S. If ΦS has column rank |S| (|S| < m), then Φ†S = (ΦTS ΦS )−1 ΦTS will denote the Moore-Penrose pseudoinverse of ΦS . PS = ΦS Φ†S will denote the projection operator on span(ΦS ) and P⊥ S = I− PS will denote the projection operator on the orthogonal complement of span(ΦS ). B. Lemmas The following lemmas will be useful for analysis of our algorithm.
3
Lemma 2.1 (monotonicity of RIC, Lemma 1 of [12]). If a measurement matrix satisfies RIP of orders K1 , K2 and K1 ≤ K2 , then δK1 ≤ δK2 . Lemma 2.2 (Lemma 1 of [12]). If x ∈ Rn is a vector with support S1 , and S1 ∩ S2 = ∅, then, if δ|S1 |+|S2 | < 1, kΦTS1 Φxk2 ≤ δ|S1 |+|S2 | kxk2
IV. S IGNAL RECOVERY
USING G OMPAM OLS ALGORITHM
We seek the conditions under which gOMPamOLS algorithm is guaranteed to recover a K-sparse signal within K iterations. We say that gOMPamOLS algorithm makes a success at an iteration if it selects at least one correct index at that iteration. Success at the first iteration: The contraction step gives the set S 1 which consists of the indices corresponding to the N largest coordinates of ΦT y. Now, k = 0 =⇒ kP⊥ T k φi k2 = kφi k2 = 1. Hence
Lemma 2.3 (Lemma 3 of [14]). For any Λ ⊂ H, define AΛ := P⊥ Λ Φ. If I1 , I2 ⊂ H such that I1 ∩ I2 = ∅, and δ|I1 |+|I2 | < 1 then, ∀u ∈ Rn such that supp(u) ⊂ I2 , δ|I1 |+|I2 | 1− kuk22 ≤ kAI1 uk22 ≤ (1 + δ|I1 |+|I2 | )kuk22 1 − δ|I1 |+|I2 | and, ! δ|I1 |+|I2 | 2 1− kΦuk22 ≤ kAI1 uk22 ≤ (1 + δ|I1 |+|I2 | )kΦuk22 1 − δ|I1 |+|I2 |
T 1 = H1
X | φi , r0 |2 = arg max ⊥ 2 S⊂S 1 :|S|=L i∈S kPT k φi k2 X | hφi , yi |2 = arg max S⊂S 1 :|S|=L i∈S
= arg max kΦTS yk22
(1)
S⊂S 1 :|S|=L
III. P ROPOSED A LGORITHM
So, kΦTT 1 yk =
TABLE II: gOMPamOLS A LGORITHM Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅, set selected by gOMP S 0 = ∅, While (krk k2 ≥ ǫ and k < K) k =k+1 X k−1 2 Contract: hk = arg max kφT k2 i r S:|S|=N i∈S yk22 Identify: hk = arg min kP⊥ T k−1 ∪{i} i∈S k Augment: T k = T k−1 ∪ hk k
Estimate: x
=
arg min
This results in the following
(2)
1
(3)
Thus, kΦTT 1 yk2
ky −
Φuk2 Update: rk = y − ΦT k xk
r
L kΦTT yk2 K r L = kΦTT ΦT xT k2 K r L ≥ (1 − δK )kxk2 K ≥
(4)
Also,
Output: estimated support set Tˆ = arg min kxk − xkS k2 and ˆ satisfying K-sparse signal x
kΦTS yk2
1 1 1 √ kΦTT 1 yk2 ≥ √ kΦTS 1 yk2 ≥ √ kΦTT yk2 L N K
u:u∈Rn , supp(u)⊂T k
End While
max
S⊂S 1 :|S|=L
S:|S|=K ˆ Tˆ = xkˆ , x T
ˆ H\Tˆ = 0 x
Our proposed algorithm is described in Table. II. It is different from the conventional mOLS algorithm (described in Table. Ib) in the identification step as well as a new step, that we call Contraction. A conventional mOLS algorithm searches for the new index in H. This search procedure is conducted by going to each index and measuring the norm of residue if this index were used in the augmentation step. The L indices which give the minimum residue norms arranged in ascending order, are chosen. For large n, this step itself is computationally expensive. We modify this step by breaking it into a contraction step and an identification step. The contraction step uses gOMP (descried in Table. Ia) to contract the set of potential candidates to a set of cardinality N (< n), choosing the indices corresponding to the N largest coordinates of ΦT rk−1 , norm wise. The resulting set is then considered for the mOLS identification step.
kΦTS 1 yk2
r
N kΦTT yk2 K r N = kΦTT ΦT xT k2 K r N (1 − δK )kxk2 ≥ K
≥
(5)
Now if no correct index were selected in S 1 , then the first iteration is bound to fail. Then kΦTS 1 yk2 = kΦTS 1 ΦT xT k2 ≤ δK+N kxk2
(6)
So, S 1 is guaranteed to contain at least one correct index if r N δK+N < (1 − δK ) √ √ √K =⇒ KδK+N + N δK < N (7) 1 This uses the fact that the average energy of the elements of set T 1 is greater than the average energy of the elements in the set S 1 which in turn is greater than the average energy of the elements in the set T
4
Even if S 1 consists of at least one correct index, no correct index may not be selected at the first iteration. Then kΦTT 1 yk2 = kΦTT 1 ΦT xT k2 ≤ δK+L kxk2
(8)
Thus sufficient condition for success in the first iteration, given the condition in Eq. (7) is guaranteed by r L δK+L < (1 − δK ) √ √K √ (9) =⇒ KδK+L + LδK < L Thus, success at the first iteration is guaranteed under the conditions in Eq. (7) and Eq. (9). Now, if we have √ L √ δN +K < √ (10) L+ K then,√ since L ≤ N , by monotonicity δL+K ≤ δN +K < L √ √ which satisfies Eq. (7). Also, since L ≤ N , under L+ K √ √ L N √ ≤√ √ the condition in Eq. (10), δN +K < √ L+ K N+ K and hence condition in Eq. (9) is satisfied. Thus under the condition in Eq. (10) the first iteration is guaranteed to succeed. Success at the (k + 1)th iteration (k > 0): We assume that the algorithm selects at least one index at each of the previous k(≥ 1) iterations. Let |T ∩ T k | = ck . Then by the induction hypothesis, ck ≥ k. Also, let that not all the correct indices are recovered, i.e. ck < K. Also, we define mk := |Sk ∩ T |, ∀k ≥ 1. Then if the ith iteration succeeded, then mi ≥ 1. The necessary conditions for the success of the (k + 1)th iteration are the following: S k+1 ∩ T 6= ∅ H k+1 ∩ S k+1 ∩ T 6= φ
(11) (12)
We will find the sufficient conditions that would ensure these. Conditions to ensure Eq. (11): We introduce some notations necessary in the analysis henceforth. The notations are drawn from [7] and [14]. We denote 2
β1k = max φi , rk (13) i∈T
αkN = min | φi , rk | where arg max
kΦTS rk k22
we find rk =y − ΦT k Φ†T k y =ΦT xT − PT k ΦT xT
=(ΦT ∩T k xT ∩T k + ΦT \T k xT \T k ) − PT k (ΦT ∩T k xT ∩T k + ΦT \T k xT \T k )
=(ΦT ∩T k xT ∩T k + ΦT \T k xT \T k )
− ΦT ∩T k xT ∩T k − PT k ΦT \T k xT \T k
=ΦT \T k xT \T k − ΦT k uT k =ΦT k ∪T x′T k ∪T
where x′T k ∪T
x k = T \T −uT k
Then αkN = min
i∈W k+1
2 φi , rk
=
kΦTW k+1 rk k22 N kΦTW k+1 \T k rk k22
=
kΦTW k+1 \T k ΦT ∪T k x′T ∪T k k22
≤
N
(∵ ΦTT k rk = 0)
N 1 2 ≤ δLk+K−ck +N kx′T ∪T k k2 N 1 2 kx′T ∪T k k2 (15) ≤ δLK+N N The last step uses the fact that ck ≥ k, k ≤ K along with monotonicity of RIC. To find a lower bound of β1k , we proceed as follows
(14)
≥
Then at (k + 1)th iteration, we can see that the set S k+1 is guaranteed to consist of at least one index from T if we have β1k > αkN . We can find a condition for this along the lines of [14]. As noted in [14], it can be seen that rk ∈ Span(ΦT ∪T k ). This is because, looking at the relationship between T, T k in Fig. 1 and keeping in mind the property of the projection operator,
=
S⊂H\T :|S|=N
β1k = kΦTT rk k2∞
i∈W k+1
W k+1 :=
Fig. 1: Venn diagram for S k+1 , T, T k
= ≥ ≥
kΦTT rk k22 K 1 k[ΦT ΦT k \T ]T rk k22 K 1 kΦTT ∪T k ΦT ∪T k x′T ∪T k k22 K 1 (1 − δLk+K−ck )2 kx′T ∪T k k22 K 1 (1 − δLK+N )2 kx′T ∪T k k22 K
(16)
5
Thus from Eq. (15) and Eq. (16) we see that √ N √ δLK+N < √ K+ N
On the other hand (17)
will assure that atleast one correct index from T is chosen in S k+1 . Condition to ensure Eq. (12): Eq. (17) guarantees mk+1 = |T ∩ S k+1 | ≥ 1. Now, the identification step in the gOMPamOLS algorithm can be expressed in a way convenient to our analysis by the virtue of the following lemma Lemma 4.1. At the (k + 1)th iteration, the identification step chooses the set
X | φi , rk | k+1 h = arg max ⊥ S:S⊂S k+1 ,|S|=L i∈S kPT k φi k2
X | φi , rk |2 = arg max ⊥ 2 S:S⊂S k+1 ,|S|=L i∈S kPT k φi k2 Proof. This lemma is a direct consequence of Proposition 1 of [8]. Nevertheless we provide a different proof in Appendix A that utilizes properties of a vector space. Now, along the lines of [8] let us define
| φi , rk | u1 = max i∈S k+1 ∩T kP⊥k φi k2 T
| φi , rk | vL = min where i∈V k+1 kP⊥k φi k2 T
X | φi , rk |2 k+1 V = arg max ⊥ 2 S⊂S k+1 \T :|S|=L i∈S kPT k φi k2
| φi , rk | vL = min i∈V k+1 kP⊥k φi k2 T
mini∈V k+1 | φi , rk | ≤ mini∈V k+1 kP⊥ φk Tk i 2
(18)
(19)
Then if we can ensure that u1 > vL , we can ensure the condition in Eq. (12). Now,
| φi , rk | u1 = max i∈S k+1 ∩T kP⊥k φi k2
T ≥ max | φi , rk | i∈S k+1 ∩T
kΦTS k+1 ∩T rk k2 √ mk T kΦT \T k rk k2 ≥ √ K − ck
≥
Now, ∀i ∈ S k+1 \ T , we have kP⊥ T k φi k2 = kAT k ei k2 2 ! δLk+1 ≥ 1− kφi k22 1 − δLk+1 2 ! δLk+1 = 1− 1 − δLk+1 and
1 min | φi , rk | ≤ √ kΦTV k+1 rk k2 k+1 i∈V L 1 ≤ √ kΦTS k+1 \T rk k2 L 1 = √ kΦTS k+1 \(T ∪T k ) rk k2 (∵ ΦTT k rk = 0) L 1 ≤ √ kΦTS k+1 \(T ∪T k ) ΦT ∪T k x′T ∪T k k2 L 1 ≤ √ δN −mk +Lk+K−ck kx′T ∪T k k2 L 1 ≤ √ δN +LK−1 kx′T ∪T k k2 (23) L The last two steps follow from the observations |(S k+1 \ (T ∪ T k )) ∪ (T ∪ T k )| = |S k+1 \ (T ∪ T k )| + |T ∪ T k | ≤ |S k+1 \ T | + |T ∪ T k | and mk ≥ 1, ck ≥ k, k ≤ K. Thus 1 vL ≤ s δ kx′ k k 2 N +LK−1 T ∪T 2 δLk+1 L 1 − 1−δ Lk+1
1 δ kx′ k k (24) <s 2 N +LK−1 T ∪T 2 δLK+N −1 L 1 − 1−δ LK+N −1
Then u1 > vL will be satisfied if we have 1 s δ kx′ k k 2 N +LK−1 T ∪T 2 δLK+N −1 L 1 − 1−δ LK+N −1
The last step results because energy of each element in S k+1 ∩ T is greater than individual energy of rest of the elements of T \ T k . Now, kΦTT \T k rk k2 = k[ΦT \T k
ΦT k ]T rk k2
=⇒
= kΦTT ∪T k ΦT ∪T k x′T ∪T k k2
≥ (1 − δLK+N −1 )kx′T ∪T k k2 (From 16) (20)
Hence 1 u1 ≥ √ (1 − δLK+N )kx′T ∪T k k2 K
(22)
r
K L
1 < √ (1 − δLK+N −1 )kx′T ∪T k k2 K √ δLK+N −1 1 − x2 x := (25) ≤ x 1 − δLK+N −1
Thus (25) is satisfied under δLK+N −1
(21)
√ L √
Sparse Signal Recovery Using gOMP Assisted mOLS
arXiv:1511.08575v1 [cs.IT] 27 Nov 2015
Samrat Mukhopadhyay1, Siddhartha Satpathi2 and Mrityunjoy Chakraborty3 Department of Electronics and Electrical Communication Engineering, Indian Institute of Technology, Kharagpur, INDIA 1 Email: [email protected],[email protected], 3 [email protected] Abstract—Because of fast convergence in finite number of steps and low computational complexity, signal recovery from compressed measurements using greedy algorithms have generated a large amount of interest in recent years. Among these greedy algorithms OMP is well studied and recently its generalization, gOMP, have also drawn attention. On the other hand OLS and its generalization mOLS have been studied in recent years because of their potential advantages in signal recovery guarantee compared to OMP or gOMP. But OLS and mOLS have the shortcomings of high computational complexity. In this paper we propose a new algorithm which uses gOMP to preselect a N length set out of which mOLS selects its L best coefficients. We have shown that our new algorithm, named gOMPamOLS, is guaranteed to reconstruct a K sparse signals perfectly from compressed measurements if the restricted √ L isometry constant satisfies δLK+N < √ . Moreover √ L+ K +L experimental results indicate that gOMPamOLS is robust under sensing matrix with correlated dictionary in terms of signal reconstruction and is much faster than OLS. Index Terms—Compressive Sensing, gOMP, mOLS, restricted isometry property
I. I NTRODUCTION There has been considerable interest for some time to recover signals from their partial measurements. Mathematically the problem can be expressed as finding solution for the following set of linear equations y = Φx where Φ ∈ Rm×n and m < n. Since this is an undetermined set of linear equations, the number of solutions is infinite, without further information about structure of x. A natural problem that then follows is to find solutions for these set of equations which has the lowest sparsity i.e. lowest number of nonzero elements, i.e the constarined l0 norm minimization problem min kxk0
x∈Rn
s.t. y = Φx However, this problem requires combinatorial search which is computationally intractable. The way out to this problem was found by Candes and Tao [1] in a seminal paper where it was shown that the l0 norm can be solved exactly using l1 norm optimization if the number of measurements is O(K ln n) where K is the sparsity of the vector x. The problem that
is required to solve under this condition is min kxk1
x∈Rn
s.t. y = Φx This is called the basis pursuit(BP) [2]. This can be solved by standard Linear programming approach. A. Greedy approaches Linear programming approach has cubic complexity on the number of columns n, which can be too slow for large n (∼ 1000 is common in many applications). In that case greedy approaches play the role of useful alternatives. A greedy algorithm utilizes some greedy rule to reconstruct the vector at each step of the algorithm. They are often faster than LP. Orthogonal Matching Pursuit(OMP) [3][4], Orthogonal Least Squares(OLS) [5][6] are among the earliest of these greedy algorithms. These algorithms greedily select indices at each step of the algorithm and append them to the already constructed support create a succesively increasing support over which they take projections to get reconstructed signal. The only difference between these two algorithms is in the method of selecting an index in the identification step. At an iteration k OMP chooses an index i which maximizes the absolute
correlation rk−1 , φi where rk−1 is the residual error vector at step k − 1 and φi is a column of the measurement matrix Φ. In OLS the identification is performed by choosing an index i which would minimize the resulting residual vector norm. Wang et al introduced gOMP algorithm [7] and mOLS algorithm [8] as generalizations of OMP and OLS algorithms respectively. In gOMP and mOLS in the identification step, the algorithms choose multiple instead of single index. See Table. I for details. The motivation for introducing both gOMP as well as mOLS is that the time complexity of the algorithms can be markedly improved if at each iteration the algorithms choose multiple indices instead of just one. This increases the probability of choosing a correct index at an iteration and hence, for a given sparse signal with sparsity K, the algorithms converge faster than OMP or OLS. Wang et al [7] has shown that the computational complexity for OMP is about 2Kmn + 3K 2 m floating point operations per second(flops) where m × n is the size of the measurement matrix and the computational complexity of gOMP is about 2smn + (2N 2 + N )s2 m flops where N is the number of indices that gOMP chooses at a time at each iteration and s is the number of iterations required for the algorithm to stop. By extensive simulation they have shown that gOMP indeed takes a small number of
2
TABLE I: gOMP
AND
mOLS A LGORITHM
(a) gOMP A LGORITHM
(b) mOLS A LGORITHM
Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅ While (krk k2 ≥ ǫ and k < K) k = k+1 X k−1 2 Identify: hk = arg max kφT k2 i r S:|S|=N i∈S
Augment: T k = T k−1 ∪ hk Estimate: xk = arg min Update: rk =
Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅ While (krk k2 ≥ ǫ and k < K) k =k+1 X Identify: hk = arg min kP⊥ yk22 T k−1 ∪{i} S:|S|=L i∈S
Augment: T k = T k−1 ∪ hk Estimate: xk = arg min
ky − Φuk2
u:u∈Rn , supp(u)⊂T k y − Φ T k xk
End While
Update: rk =
ky − Φuk2
u:u∈Rn , supp(u)⊂T k y − Φ T k xk
End While
Output: estimated support set Tˆ = arg min kxk − xkS k2 and K-sparse S:|S|=K
Output: estimated support set Tˆ = arg min kxk − xkS k2 and K-sparse S:|S|=K
ˆ satisfying x ˆ Tˆ = xkˆ , x ˆ H\Tˆ = 0 signal x
ˆ satisfying x ˆ Tˆ = xkˆ , x ˆ H\Tˆ = 0 signal x
T
simulation s to recover the K sparse signal which effectively reduces its computational complexity. Similar conclusions are drawn for mOLS compared to OLS by Wang et al [8] through extensive simulations. Other algorithms of the same flavour as gOMP or mOLS are stagewise OMP(StOMP) [9], regularized OMP(ROMP) [10], Compressed Sampling Matching Pursuit(CoSaMP) [11] and Subspace Pursuit (SP) [12]. OMP has been the representative of these algorithms because of its simplicity and recovery performance. Tropp and Gilbert [4] have shown that OMP can recover the original sparse vector from a few measurements with exceedingly high probability when the measurement matrix has entries i.i.d Gaussian. But recently Soussen et al [13] have shown that OLS performs considerably better than OMP in recovering sparse signals from few measurements when the measurement matrix has correlated entries. They have also shown that there exists measurement dictionaries for which OMP fails to recover the original signal whereas OLS is unaffected. B. Our contributions in this paper • In this paper we propose a new algorithm, referred to as gOMPamOLS (short for gOMP assisted mOLS), that exploits the capability of gOMP to select “good” indices at a step as well as the robustness of mOLS in maintaining high recovery performance in the presence measurement matrices with correlated dictionaries. • We find a criteria under which our proposed algorithm can exactly recover the unknown K-sparse signal which is based upon the restricted isometry property of the measurement matrix. A matrix Φ ∈ Rm×n is said to have restricted isometry property (RIP) of order K if ∃ δ > 0 such that ∀x ∈ Rn with sparsity K, the matrix Φ satisfies the following inequalities (1 − δ)kxk22 ≤ kΦxk22 ≤ (1 + δ)kxk22 The smallest δ for which this property is satisfied is called the restricted isometry constant of order K of the matrix Φ and is denoted by δK . We show that gOMPamOLS with gOMP selecting N indices and mOLS selecting
T
L ≤ N indices (See Table. II for details of the algorithm) is guaranteed to recover a K sparse vector exactly from compressed measurements if the matrix Φ obeys RIP with constant √ L √ δLK+N < √ K+L+ L •
We compare simulation results for recovery performances of OMP, OLS and gOMPamOLS and show that gOMPamOLS performs much better than OMP in recovering sparse vector for measurement matrices with correlated dictionary. Also, gOMPmaOLS is much faster than OLS and thus becomes a suitable upgradation of both OMP and OLS. II. P RELIMINARIES
A. Notations The following notations will be used throughout the paper. We will denote by k · kp the lp norm of a vector v, i.e., Pn 1/p kvkp = ( i=1 |vi |p ) . Φ ∈ Rm×n will denote the measurement matrix with m rows and n columns with m < n. The i th column of Φ will be denoted by φi ; i = 1, 2, · · · , n. We assume that all the columns of Φ have unit l2 norm. K will denote the sparsity level estimated beforehand. xS will denote the vector x restricted to a subset S. H will denote the set of all the indices {1, 2, · · · , n}. T will refer to the support set of the original vector x. Similarly, ΦS will denote the submatrix of Φ formed with the columns restricted to index set S. If ΦS has column rank |S| (|S| < m), then Φ†S = (ΦTS ΦS )−1 ΦTS will denote the Moore-Penrose pseudoinverse of ΦS . PS = ΦS Φ†S will denote the projection operator on span(ΦS ) and P⊥ S = I− PS will denote the projection operator on the orthogonal complement of span(ΦS ). B. Lemmas The following lemmas will be useful for analysis of our algorithm.
3
Lemma 2.1 (monotonicity of RIC, Lemma 1 of [12]). If a measurement matrix satisfies RIP of orders K1 , K2 and K1 ≤ K2 , then δK1 ≤ δK2 . Lemma 2.2 (Lemma 1 of [12]). If x ∈ Rn is a vector with support S1 , and S1 ∩ S2 = ∅, then, if δ|S1 |+|S2 | < 1, kΦTS1 Φxk2 ≤ δ|S1 |+|S2 | kxk2
IV. S IGNAL RECOVERY
USING G OMPAM OLS ALGORITHM
We seek the conditions under which gOMPamOLS algorithm is guaranteed to recover a K-sparse signal within K iterations. We say that gOMPamOLS algorithm makes a success at an iteration if it selects at least one correct index at that iteration. Success at the first iteration: The contraction step gives the set S 1 which consists of the indices corresponding to the N largest coordinates of ΦT y. Now, k = 0 =⇒ kP⊥ T k φi k2 = kφi k2 = 1. Hence
Lemma 2.3 (Lemma 3 of [14]). For any Λ ⊂ H, define AΛ := P⊥ Λ Φ. If I1 , I2 ⊂ H such that I1 ∩ I2 = ∅, and δ|I1 |+|I2 | < 1 then, ∀u ∈ Rn such that supp(u) ⊂ I2 , δ|I1 |+|I2 | 1− kuk22 ≤ kAI1 uk22 ≤ (1 + δ|I1 |+|I2 | )kuk22 1 − δ|I1 |+|I2 | and, ! δ|I1 |+|I2 | 2 1− kΦuk22 ≤ kAI1 uk22 ≤ (1 + δ|I1 |+|I2 | )kΦuk22 1 − δ|I1 |+|I2 |
T 1 = H1
X | φi , r0 |2 = arg max ⊥ 2 S⊂S 1 :|S|=L i∈S kPT k φi k2 X | hφi , yi |2 = arg max S⊂S 1 :|S|=L i∈S
= arg max kΦTS yk22
(1)
S⊂S 1 :|S|=L
III. P ROPOSED A LGORITHM
So, kΦTT 1 yk =
TABLE II: gOMPamOLS A LGORITHM Input: measurement vector y ∈ Rm , sensing matrix Φ ∈ Rm×n , sparsity level K, true support set T = supp(x), total set H = {1, 2, · · · , n} Initialize: counter k = 0, residue r0 = y, estimated support set, T 0 = ∅, set selected by gOMP S 0 = ∅, While (krk k2 ≥ ǫ and k < K) k =k+1 X k−1 2 Contract: hk = arg max kφT k2 i r S:|S|=N i∈S yk22 Identify: hk = arg min kP⊥ T k−1 ∪{i} i∈S k Augment: T k = T k−1 ∪ hk k
Estimate: x
=
arg min
This results in the following
(2)
1
(3)
Thus, kΦTT 1 yk2
ky −
Φuk2 Update: rk = y − ΦT k xk
r
L kΦTT yk2 K r L = kΦTT ΦT xT k2 K r L ≥ (1 − δK )kxk2 K ≥
(4)
Also,
Output: estimated support set Tˆ = arg min kxk − xkS k2 and ˆ satisfying K-sparse signal x
kΦTS yk2
1 1 1 √ kΦTT 1 yk2 ≥ √ kΦTS 1 yk2 ≥ √ kΦTT yk2 L N K
u:u∈Rn , supp(u)⊂T k
End While
max
S⊂S 1 :|S|=L
S:|S|=K ˆ Tˆ = xkˆ , x T
ˆ H\Tˆ = 0 x
Our proposed algorithm is described in Table. II. It is different from the conventional mOLS algorithm (described in Table. Ib) in the identification step as well as a new step, that we call Contraction. A conventional mOLS algorithm searches for the new index in H. This search procedure is conducted by going to each index and measuring the norm of residue if this index were used in the augmentation step. The L indices which give the minimum residue norms arranged in ascending order, are chosen. For large n, this step itself is computationally expensive. We modify this step by breaking it into a contraction step and an identification step. The contraction step uses gOMP (descried in Table. Ia) to contract the set of potential candidates to a set of cardinality N (< n), choosing the indices corresponding to the N largest coordinates of ΦT rk−1 , norm wise. The resulting set is then considered for the mOLS identification step.
kΦTS 1 yk2
r
N kΦTT yk2 K r N = kΦTT ΦT xT k2 K r N (1 − δK )kxk2 ≥ K
≥
(5)
Now if no correct index were selected in S 1 , then the first iteration is bound to fail. Then kΦTS 1 yk2 = kΦTS 1 ΦT xT k2 ≤ δK+N kxk2
(6)
So, S 1 is guaranteed to contain at least one correct index if r N δK+N < (1 − δK ) √ √ √K =⇒ KδK+N + N δK < N (7) 1 This uses the fact that the average energy of the elements of set T 1 is greater than the average energy of the elements in the set S 1 which in turn is greater than the average energy of the elements in the set T
4
Even if S 1 consists of at least one correct index, no correct index may not be selected at the first iteration. Then kΦTT 1 yk2 = kΦTT 1 ΦT xT k2 ≤ δK+L kxk2
(8)
Thus sufficient condition for success in the first iteration, given the condition in Eq. (7) is guaranteed by r L δK+L < (1 − δK ) √ √K √ (9) =⇒ KδK+L + LδK < L Thus, success at the first iteration is guaranteed under the conditions in Eq. (7) and Eq. (9). Now, if we have √ L √ δN +K < √ (10) L+ K then,√ since L ≤ N , by monotonicity δL+K ≤ δN +K < L √ √ which satisfies Eq. (7). Also, since L ≤ N , under L+ K √ √ L N √ ≤√ √ the condition in Eq. (10), δN +K < √ L+ K N+ K and hence condition in Eq. (9) is satisfied. Thus under the condition in Eq. (10) the first iteration is guaranteed to succeed. Success at the (k + 1)th iteration (k > 0): We assume that the algorithm selects at least one index at each of the previous k(≥ 1) iterations. Let |T ∩ T k | = ck . Then by the induction hypothesis, ck ≥ k. Also, let that not all the correct indices are recovered, i.e. ck < K. Also, we define mk := |Sk ∩ T |, ∀k ≥ 1. Then if the ith iteration succeeded, then mi ≥ 1. The necessary conditions for the success of the (k + 1)th iteration are the following: S k+1 ∩ T 6= ∅ H k+1 ∩ S k+1 ∩ T 6= φ
(11) (12)
We will find the sufficient conditions that would ensure these. Conditions to ensure Eq. (11): We introduce some notations necessary in the analysis henceforth. The notations are drawn from [7] and [14]. We denote 2
β1k = max φi , rk (13) i∈T
αkN = min | φi , rk | where arg max
kΦTS rk k22
we find rk =y − ΦT k Φ†T k y =ΦT xT − PT k ΦT xT
=(ΦT ∩T k xT ∩T k + ΦT \T k xT \T k ) − PT k (ΦT ∩T k xT ∩T k + ΦT \T k xT \T k )
=(ΦT ∩T k xT ∩T k + ΦT \T k xT \T k )
− ΦT ∩T k xT ∩T k − PT k ΦT \T k xT \T k
=ΦT \T k xT \T k − ΦT k uT k =ΦT k ∪T x′T k ∪T
where x′T k ∪T
x k = T \T −uT k
Then αkN = min
i∈W k+1
2 φi , rk
=
kΦTW k+1 rk k22 N kΦTW k+1 \T k rk k22
=
kΦTW k+1 \T k ΦT ∪T k x′T ∪T k k22
≤
N
(∵ ΦTT k rk = 0)
N 1 2 ≤ δLk+K−ck +N kx′T ∪T k k2 N 1 2 kx′T ∪T k k2 (15) ≤ δLK+N N The last step uses the fact that ck ≥ k, k ≤ K along with monotonicity of RIC. To find a lower bound of β1k , we proceed as follows
(14)
≥
Then at (k + 1)th iteration, we can see that the set S k+1 is guaranteed to consist of at least one index from T if we have β1k > αkN . We can find a condition for this along the lines of [14]. As noted in [14], it can be seen that rk ∈ Span(ΦT ∪T k ). This is because, looking at the relationship between T, T k in Fig. 1 and keeping in mind the property of the projection operator,
=
S⊂H\T :|S|=N
β1k = kΦTT rk k2∞
i∈W k+1
W k+1 :=
Fig. 1: Venn diagram for S k+1 , T, T k
= ≥ ≥
kΦTT rk k22 K 1 k[ΦT ΦT k \T ]T rk k22 K 1 kΦTT ∪T k ΦT ∪T k x′T ∪T k k22 K 1 (1 − δLk+K−ck )2 kx′T ∪T k k22 K 1 (1 − δLK+N )2 kx′T ∪T k k22 K
(16)
5
Thus from Eq. (15) and Eq. (16) we see that √ N √ δLK+N < √ K+ N
On the other hand (17)
will assure that atleast one correct index from T is chosen in S k+1 . Condition to ensure Eq. (12): Eq. (17) guarantees mk+1 = |T ∩ S k+1 | ≥ 1. Now, the identification step in the gOMPamOLS algorithm can be expressed in a way convenient to our analysis by the virtue of the following lemma Lemma 4.1. At the (k + 1)th iteration, the identification step chooses the set
X | φi , rk | k+1 h = arg max ⊥ S:S⊂S k+1 ,|S|=L i∈S kPT k φi k2
X | φi , rk |2 = arg max ⊥ 2 S:S⊂S k+1 ,|S|=L i∈S kPT k φi k2 Proof. This lemma is a direct consequence of Proposition 1 of [8]. Nevertheless we provide a different proof in Appendix A that utilizes properties of a vector space. Now, along the lines of [8] let us define
| φi , rk | u1 = max i∈S k+1 ∩T kP⊥k φi k2 T
| φi , rk | vL = min where i∈V k+1 kP⊥k φi k2 T
X | φi , rk |2 k+1 V = arg max ⊥ 2 S⊂S k+1 \T :|S|=L i∈S kPT k φi k2
| φi , rk | vL = min i∈V k+1 kP⊥k φi k2 T
mini∈V k+1 | φi , rk | ≤ mini∈V k+1 kP⊥ φk Tk i 2
(18)
(19)
Then if we can ensure that u1 > vL , we can ensure the condition in Eq. (12). Now,
| φi , rk | u1 = max i∈S k+1 ∩T kP⊥k φi k2
T ≥ max | φi , rk | i∈S k+1 ∩T
kΦTS k+1 ∩T rk k2 √ mk T kΦT \T k rk k2 ≥ √ K − ck
≥
Now, ∀i ∈ S k+1 \ T , we have kP⊥ T k φi k2 = kAT k ei k2 2 ! δLk+1 ≥ 1− kφi k22 1 − δLk+1 2 ! δLk+1 = 1− 1 − δLk+1 and
1 min | φi , rk | ≤ √ kΦTV k+1 rk k2 k+1 i∈V L 1 ≤ √ kΦTS k+1 \T rk k2 L 1 = √ kΦTS k+1 \(T ∪T k ) rk k2 (∵ ΦTT k rk = 0) L 1 ≤ √ kΦTS k+1 \(T ∪T k ) ΦT ∪T k x′T ∪T k k2 L 1 ≤ √ δN −mk +Lk+K−ck kx′T ∪T k k2 L 1 ≤ √ δN +LK−1 kx′T ∪T k k2 (23) L The last two steps follow from the observations |(S k+1 \ (T ∪ T k )) ∪ (T ∪ T k )| = |S k+1 \ (T ∪ T k )| + |T ∪ T k | ≤ |S k+1 \ T | + |T ∪ T k | and mk ≥ 1, ck ≥ k, k ≤ K. Thus 1 vL ≤ s δ kx′ k k 2 N +LK−1 T ∪T 2 δLk+1 L 1 − 1−δ Lk+1
1 δ kx′ k k (24) <s 2 N +LK−1 T ∪T 2 δLK+N −1 L 1 − 1−δ LK+N −1
Then u1 > vL will be satisfied if we have 1 s δ kx′ k k 2 N +LK−1 T ∪T 2 δLK+N −1 L 1 − 1−δ LK+N −1
The last step results because energy of each element in S k+1 ∩ T is greater than individual energy of rest of the elements of T \ T k . Now, kΦTT \T k rk k2 = k[ΦT \T k
ΦT k ]T rk k2
=⇒
= kΦTT ∪T k ΦT ∪T k x′T ∪T k k2
≥ (1 − δLK+N −1 )kx′T ∪T k k2 (From 16) (20)
Hence 1 u1 ≥ √ (1 − δLK+N )kx′T ∪T k k2 K
(22)
r
K L
1 < √ (1 − δLK+N −1 )kx′T ∪T k k2 K √ δLK+N −1 1 − x2 x := (25) ≤ x 1 − δLK+N −1
Thus (25) is satisfied under δLK+N −1
(21)
√ L √
