21.4. Compressive Sampling Matching Pursuit#

In this section, we look at an algorithm named CoSaMP (Compressive Sampling Matching Pursuit) developed in [61]. This algorithm follows in the tradition of orthogonal matching pursuit while bringing in other ideas from many recent developments in the field leading to an algorithm which is much more robust, fast and provides much stronger guarantees.

This algorithm has many impressive features:

  • It can work with a variety of sensing matrices.

  • It works with minimal number of measurements (c.f. basis pursuit).

  • It is robust against presence of measurement noise (OMP does not have such a claim).

  • It provides optimal error guarantees for every signal (sparse signals, compressible signals, completely arbitrary signals).

  • The algorithm is quite efficient in terms of resource (memory, CPU) requirements.

As we move along in this section, we will understand the algorithm and validate all of these claims. Hopefully through this process we will have a very good understanding of how to evaluate the quality of signal recovery algorithm.

21.4.1. Algorithm#

The algorithm itself is presented in Algorithm 21.4.

Algorithm 21.4 (CoSaMP for iterative sparse signal recovery)

Inputs:

  • Sensing matrix ΦCM×N

  • Measurement vector: yCM where y=Φx+e

  • Sparsity level: K

Outputs:

  • z: a K-sparse approximation of the signal:xCN

Initialization:

  1. z00 # Initial approximation

  2. r0y # Residual yΦz

  3. k0 # Iteration counter

Algorithm:

  1. If halting criteria is satisfied: break.

  2. kk+1.

  3. pΦHrk1 # Form signal proxy

  4. Ω=supp(p|2K) # Identify 2K large components

  5. TΩsupp(zk1) # Merge supports

  6. bTΦTy # Estimation by least-squares

  7. bTc0

  8. zkb|K # Prune to obtain next approximation

  9. rkyΦzk # Update residual

  10. Go to step 1.

Let us set out the notation before proceeding further. Most of the things are as usual, with few minor updates.

  • xCN represents the signal which is to be estimated through the algorithm. x is unknown to us within the algorithm.

  • N is the dimension of ambient signal space, KN is the sparsity level of of the approximation of x that we are estimating and M is the dimension of measurement space (number of measurements K<MN).

  • We note that x itself may not be K-sparse.

  • Our algorithm is designed to estimate a K-sparse approximation of x.

  • ΦCM×N represents the sensing matrix. It is known to the algorithm.

  • y=Φx+e represents the measurement vector belonging to CM. This is known to the algorithm.

  • eCM represents the measurement noise which is unknown to us within the algorithm.

  • k represents the iteration (or step) counter within the algorithm.

  • zCN represents our estimate of x. z is updated iteratively.

  • zk represents the estimate of x at the end of k-th iteration.

  • We start with z0=0 and update it in each iteration.

  • pCN represents a proxy for xzk1. We will explain it shortly.

  • T, Ω etc. represent index sets (subsets of {1,2,,N}).

  • For any vCN and any index set T{1,2,,N}, vT can mean either of the two things:

    • A vector in C|T| consisting of only those entries in v which are indexed by T.

    • A vector in CN whose entries indexed by T are same as that of v while entries indexed by {1,2,,N}T are set all to zero.

      (21.19)#vT(i)={v(i)if iT;0otherwise.
  • For any vCN and any integer 1nN, v|n means a vector in CN which consists of the n largest (in magnitude) entries of v (at the corresponding indices) while rest of the entries in v|n are 0.

  • With an index set T, ΦT means an M×|T| matrix consisting of selected columns of Φ indexed by T.

  • rCM represents the difference between the actual measurement vector y and estimated measurement vector Φz.

  • Ideal estimate of r would be e itself. But that won’t be possible to achieve in general.

  • x|K is the best K-sparse approximation of x which can be estimated by our algorithm (Definition 18.13).

Example 21.1 (Clarifying the notation in CoSaMP)

  1. Let us consider

    x=(1,5,8,0,0,3,0,0,0,0)
  2. Then

    x|2=(0,5,8,0,0,0,0,0,0,0)
  3. Also

    x|4=x

    since x happens to be 4-sparse.

  4. For Λ={1,2,3,4}, we have

    x{1,2,3,4}=(1,5,8,0,0,0,0,0,0,0)

    or

    x{1,2,3,4}=(1,5,8,0)

    in different contexts.

21.4.1.1. The Signal Proxy#

As we have learnt by now that the most challenging part in a signal recovery algorithm is to identify the support of the K-sparse approximation. OMP identifies one index in the support at each iteration and hopes that it never makes any mistakes. It ends up taking K iterations and solving K least squares problems. If there could be a simple way which could identify the support quickly (even if roughly), that can help in tremendously accelerating the algorithm. The fundamental innovation in CoSaMP is to identify the support through the signal proxy.

  1. If x is K-sparse and Φ satisfies RIP (see Restricted Isometry Property) of order K with the restricted isometry constant δK1 then it can be argued that p=ΦHΦx can serve as a proxy for the signal x.

  2. In particular the largest K entries of p point towards the largest K entries of x.

  3. Although we don’t know x inside the algorithm, yet we have access to y=Φx (assuming error to be 0), the proxy can be easily estimated by computing p=ΦHy.

This is the key idea in CoSaMP.

  1. Rather than identifying just one new index in support of x, it tries to identify whole of support by picking up the largest 2K entries of p.

  2. It then solves a least squares problem around the columns of Φ indexed by this support set.

  3. It keeps only the K largest entries from the least squares solution as an estimate of x.

  4. Of course there is no guarantee that in a single attempt the support of x will be recovered completely.

  5. Hence the residual between actual measurement vector y and estimated measurement vector Φz is computed and it is used to identify other indices of supp(x) iteratively.

21.4.1.2. Core of Algorithm#

There are two things which are estimated in each iteration of algorithm

  1. K-sparse estimate of x at k-th step: zk.

  2. Residual at k-th step: rk.

Residual and proxy:

  1. We start with a trivial estimate z0=0 and improve it in each iteration.

  2. r0 is nothing but the measurement vector y.

  3. As explained before rk is the difference between actual measurement vector y and the estimated measurement vector Φzk.

  4. This rk is used for computing the signal proxy at k-th step.

  5. Concretely assuming e=0

    p=ΦHrk=ΦH(yΦzk)=ΦHΦ(xzk).
  6. Thus p is the proxy of the difference between the original signal and the estimated signal.

During each iteration, the algorithm performs following tasks

  1. [Identification] The algorithm forms a proxy of the residual rk1 and locates the 2K largest entries of the proxy.

  2. [Support merger] The set of newly identified indices is merged with the set of indices that appear in the current approximation zk1 (i.e. supp(zk1)).

  3. [Estimation] The algorithm solves a least squares problem to approximate the signal x on the merged set of indices.

  4. [Pruning] The algorithm produces a new approximation zk by retaining only the largest K entries in the least squares solution.

  5. [Residual update] Finally the new residual rk between original measurement vector y and estimated measurement vector Φzk is computed.

The steps are repeated until the halting criteria is reached. We will discuss the halting criteria in detail later. A possible halting criteria is when the norm of residual rk reaches a very small value.

21.4.2. CoSaMP Analysis Sparse Case#

In this subsection, we will carry out a detailed theoretical analysis of CoSaMP algorithm for sparse signal recovery.

We will make following assumptions in the analysis:

  • The sparsity level K is fixed and known in advance.

  • The signal x is K-sparse (i.e. xΣKCN).

  • The sensing matrix Φ satisfies RIP of order 4K with δ4K0.1.

  • Measurement error eCM is arbitrary.

For each iteration, we need to define one more quantity: recovery error

dk=xzk.

The quantity dk captures the difference between actual signal x and estimated signal zk at the end of k-th iteration.

  1. Under ideal recovery, dk should become 0 as k increases.

  2. Since we don’t know x within the algorithm, we cannot see dk directly.

  3. We do have following equation

    rk=yΦzk=Φx+eΦzk=Φ(xzk)+e=Φdk+e.
  4. We will show that in each iteration, CoSaMP reduces the recovery error by a constant factor while adding a small multiple of the measurement noise e2.

  5. As a result, when recovery error is large compared to measurement noise, the algorithm makes substantial progress in each step.

  6. The algorithm stops making progress when recovery error is of the order of measurement noise.

We will consider some iteration k1 and analyze the behavior of each of the steps: identification, support merger, estimation, pruning and residual update one by one.

  1. At the beginning of the iteration we have

    rk1=Φ(xzk1)+e=Φdk1+e.
  2. The quantities of interest are zk1, rk1 and dk1 out of which dk1 is unknowable.

  3. In order to simplify the equations below we will write

    r=rk1,z=zk1 and d=dk1.
Table 21.1 Terms used in analysis#

Term

Description

Details

x

K sparse signal

Φ

Sensing matrix

e

Measurement noise

y

Measurement vector

y=Φx+e

z

A K sparse estimate of x at the beginning of the iteration

z=zk1=b|K, bT=ΦTy

r

Residual at the beginning of the iteration

r=rk1=yΦzk1

d

Recovery error at the beginning of the iteration

d=dk1=xzk1

p

The proxy for d

pΦHr=ΦH(Φdk1+e)

Ω

The index set of largest 2K entries in p

21.4.2.1. Identification#

We start our analysis with the identification step in the main loop of CoSaMP (Algorithm 21.4).

  1. We compute p=ΦHr and consider Ω=supp(p|2K) as the index set of largest 2K indices in p.

  2. We will show that most of the energy in d (the recovery error) is concentrated in the entries indexed by Ω.

  3. In other words, we will be looking for a bound of the form

    dΩc2d2.
  4. Since x is K-sparse and z is K-sparse hence d is 2K-sparse.

    1. In first iteration where z0=0, support of d0 is same as support of x.

    2. In later iterations it may include more indices.

  5. Let

    Γ=supp(d).
  6. We remind that Ω is known to us while Γ is unknown to us.

  7. Ideal case would be when ΓΩ.

  8. Then we would have recovered whole of support of x; recovering x would therefore be easier.

  9. Let us examine what are the differences between the two sets.

  10. We have |Ω|2K and |Γ|2K.

  11. Moreover, Ω indexes 2K largest entries in p.

  12. Thus we have (due to Theorem 18.18)

    (21.20)#pΓ2pΩ2

    where pΓ,pΩCN are obtained from p as per Definition 18.11.

  13. Squaring (21.20) on both sides and expanding we get

    γΓ|pγ|2ωΩ|pω|2.
  14. Now we do hope that some of the entries are common in both sides. Those entries are indexed by ΓΩ.

  15. The remaining entries on L.H.S. are indexed by ΓΩ and on the R.H.S. are indexed by ΩΓ.

  16. Hence we have

    (21.21)#pΓΩ22pΩΓ22pΓΩ2pΩΓ2.
  17. For both ΓΩ and ΩΓ we have

    |ΓΩ|2K and |ΩΓ|2K.
  18. The worst case happens when both sets are totally disjoint.

  19. Let us first examine the R.H.S and find an upper bound for it.

    pΩΓ2=ΦΩΓHr2=ΦΩΓH(Φd+e)2ΦΩΓHΦd2+ΦΩΓHe2.

    At this juncture, its worthwhile to scan various results in Restricted Isometry Property since we are going to need many of them in the following steps.

  20. Let us look at the term ΦΩΓHΦd2.

  21. We have |ΩΓ|2K.

  22. Further d is 2K sparse and sits over Γ which is disjoint with ΩΓ.

  23. Together |ΓΩ|4K.

  24. A straightforward application of Corollary 18.5 gives us

    ΦΩΓHΦd2δ4Kd2.
  25. Similarly using Theorem 18.57 we get

    ΦΩΓHe21+δ2Ke2.
  26. Combining the two we get

    (21.22)#pΩΓ2δ4Kd2+1+δ2Ke2.
  27. We now look at L.H.S. and find a lower bound for it.

    pΓΩ2=ΦΓΩHr2=ΦΓΩH(Φd+e)2.
  28. We will split d as

    d=dΓΩ+dΩ.
  29. Further we will use a form of triangular inequality as

    a+bab.
  30. We expand L.H.S.

    ΦΓΩH(Φd+e)2=ΦΓΩH(Φ(dΓΩ+dΩ)+e)2ΦΓΩHΦdΓΩ2ΦΓΩHΦdΩ2ΦΓΩHe2.
  31. As before

    ΦΓΩHΦdΓΩ2=ΦΓΩHΦΓΩdΓΩ2.
  32. We can use the lower bound from Theorem 18.59 to give us

    ΦΓΩHΦΓΩdΓΩ2(1δ2K)dΓΩ2

    since |ΓΩ|2K.

  33. For the other two terms we need to find their upper bounds since they appear in negative sign.

  34. Applying Corollary 18.5 we get

    ΦΓΩHΦdΩ2δ2Kd2.
  35. Again using Theorem 18.57 we get

    ΦΓΩHe21+δ2Ke2.
  36. Putting the three bounds together, we get a lower bound for L.H.S.

    (21.23)#pΓΩ2(1δ2K)dΓΩ2δ2Kd21+δ2Ke2.
  37. Finally plugging the upper bound from (21.22) and lower bound from (21.23) into (21.21) we get

    (1δ2K)dΓΩ2δ2Kd21+δ2Ke2δ4Kd2+1+δ2Ke2(1δ2K)dΓΩ2(δ2K+δ4K)d2+21+δ2Ke2dΓΩ2(δ2K+δ4K)d2+21+δ2Ke21δ2K.
  38. This result is nice but there is a small problem. We would like to get rid of Γ from L.H.S. and get dΩc instead.

  39. Recall that

    Ωc=(ΓΩc)(ΓcΩc)

    where ΓΩc and ΓcΩc are disjoint.

  40. Thus

    dΩc=dΓΩc+dΓcΩc.
  41. Since d is 0 on Γc (recall that Γ=supp(d)), hence

    dΓcΩc=0.
  42. As a result

    dΩc=dΓΩc=dΓΩ.
  43. This gives us

    (21.24)#dΩc2(δ2K+δ4K)d2+21+δ2Ke21δ2K.
  44. Let us simplify this equation by using δ2Kδ4K0.1 assumed at the beginning of the analysis.

  45. This gives us

    1δ2K0.9,δ2K+δ4K0.2,21+δ2K21.1=2.098.
  46. Thus we get

    dΩc20.2d2+2.098e2.9.
  47. Simplifying

    (21.25)#dΩc20.223d2+2.331e2.
  48. This result tells us that in the absence of measurement noise, at least 78% of energy in the recovery error is indeed concentrated in the entries indexed by Ω.

  49. Moreover, if the measurement noise is small compared to the size of recovery error, then this concentration continues to hold.

  50. Thus even if we don’t know Γ directly, Ω is a close approximation for Γ.

We summarize the analysis for the identification step in the following lemma.

Lemma 21.2

Let Φ satisfy RIP of order 4K with δ4K0.1. At every iteration k in CoSaMP algorithm, let dk1=(xzk1) be the recovery error and p=ΦHrk1 be the signal proxy. The set Ω=supp(p|2K) contains at most 2K indices and

(21.26)#dΩck120.223dk12+2.331e2.

In other words, most of the energy in dk1 is concentrated in the entries indexed by Ω when measurement noise is small.

21.4.2.2. Support Merger#

The next step in a CoSaMP iteration is support merger. Recalling from Algorithm 21.4, the step involves computing

T=Ωsupp(z)

i.e. we merge the support of current signal estimate z with the newly identified set of indices Ω.

In previous lemma we were able to show how well the recovery error is concentrated over the index set Ω. But we didn’t establish anything concrete about how x is concentrated. In the support merger step, we will be able to show that x is highly concentrated over the index set T. For this we will need to find an upper bound on xTc2 i.e. the energy of entries in x which are not covered by T.

  1. We recall that |Ω|2K and since z is a K-sparse estimate of x hence |supp(z)|K.

  2. Thus

    |T|3K.
  3. Clearly

    Tc=Ωcsupp(z)cΩc.
  4. Further since supp(z)T hence zTc=0.

  5. Thus

    xTc=(xz)Tc=dTc.
  6. Since TcΩc hence

    dTc2dΩc2.
  7. Combining these facts we can write

    xTc2=dTc2dΩc2.
  8. We summarize this analysis in following lemma.

Lemma 21.3

Let Ω be a set of at most 2K entries. The set T=Ωsupp(zk1) contains at most 3K entries, and

(21.27)#xTc2dΩc2.

Note that this inequality says nothing about how Ω is chosen. But we can use an upper bound on dΩc2 from Lemma 21.2 to get an upper bound on xTc2.

21.4.2.3. Estimation#

The next step is a least square estimate of x over the columns indexed by T. Recalling from Algorithm 21.4 we compute

bT=ΦTy=ΦT(Φx+e).
  1. We also set bTc=0.

  2. We have

    b=bT+bTc=bT.
  3. Since |T|3K, b is a 3K-sparse approximation of x.

  4. What we need here is a bound over the approximation error xb2.

  5. We have already obtained a bound on xTc2.

  6. If an upper bound on xb2 depends on xTc2 and measurement noise e2 that would be quite reasonable.

  7. We start with splitting x over T and Tc.

    x=xT+xTc.
  8. Since supp(b)T hence

    xb=xT+xTcbT=(xTbT)+xTc.
  9. This gives us

    xb2xTbT2+xTc2.
  10. We will now expand the term xTbT2.

    xTbT2=xTΦT(Φx+e)2=xTΦT(ΦxT+ΦxTc+e)2
  11. Now since ΦT is full column rank (since Φ satisfies RIP of order 3K) hence ΦTΦT=I.

  12. Also ΦxT=ΦTxT (since column columns indexed by T are involved).

  13. This helps us cancel xTΦTΦxT.

  14. Thus

    xTbT2=ΦT(ΦxTc+e)2(ΦTHΦT)1ΦTHΦxTc2+ΦTe2.
  15. Let us look at the terms on R.H.S. one by one.

  16. Let v=ΦTHΦxTc.

  17. Then

    (ΦTHΦT)1ΦTHΦxTc2=(ΦTHΦT)1v2.
  18. This perfectly fits Theorem 18.59 with |T|3K giving us

    (ΦTHΦT)1ΦTΦxTc211δ3KΦTHΦxTc2.
  19. Further, applying Corollary 18.5 we get

    ΦTHΦxTc2δ4KxTc2

    since |Tsupp(x)|4K.

  20. For the measurement noise term, applying Theorem 18.58 we get

    ΦTe211δ3Ke2.
  21. Combining the above inequalities we get

    xb2[1+δ4K1δ3K]xTc2+11δ3Ke2.
  22. Recalling our assumption that δ3Kδ4K0.1 we can simplify the constants to get

    xb21.112xTc2+1.0541e2.

We summarize the analysis in this step in the following lemma.

Lemma 21.4

Let T be a set of at most 3K indices, and define the least squares signal estimate b by the formula

bT=ΦTy=ΦT(Φx+e) and bTc=0.

If Φ satisfies RIP of order 4K with δ4K0.1 then the following holds

xb21.112xTc2+1.0541e2.

Note that the lemma doesn’t make any assumptions about how T is arrived at except that |T|3K. Also, this analysis assumes that the least squares solution is of infinite precision given by ΦTy. The approximation error in an iterative least squares solution needs to be separately analyzed to identify the number of required steps so that the least squares error is negligible.

21.4.2.4. Pruning#

The last step in the main loop of CoSaMP (Algorithm 21.4) is pruning.

We compute

zk=b|K

as the next estimate of x by picking the K largest entries in b.

  1. We note that both x and b|K can be regarded as K term approximations of b.

  2. As established in Theorem 18.18, b|K is the best K-term approximation of b.

  3. Thus

    bb|K2bx2.
  4. Now

    xzk2=xb+bb|K2xb2+bb|K22xb2.
  5. This helps us establish that although the 3K-sparse approximation b is closer to x compared to b|K, but b|K is also not too bad approximation of x while using only K entries at most.

We summarize it in following lemma.

Lemma 21.5

The pruned estimate zk=b|K satisfies

(21.28)#xzk2=xb|K22xb2.

21.4.2.5. The CoSaMP Iteration Invariant#

Having analyzed each of the steps in the main loop of CoSaMP algorithm, it is time for us to combine the analysis. Concretely, we wish to establish how much progress CoSaMP makes in each iteration. Following theorem provides an upper bound on how the recovery error changes from one iteration to next iteration. This theorem is known as the iteration invariant for the sparse case.

Theorem 21.4 (CoSaMP iteration invariant for exact sparse recovery)

Assume that x is K-sparse. Assume that Φ satisfies RIP of order 4K with δ4K0.1. For each iteration number k1, the signal estimate zk is K-sparse and satisfies

xzk212xzk12+7.5e2.

In particular

xzk22kx2+15e2.

This theorem helps us establish that if measurement noise is small, then the algorithm makes substantial progress in each iteration. The proof makes use of the lemmas developed above.

Proof. We run the proof in backtracking mode. We start from zk and go back step by step through pruning, estimation, support merger, and identification to connect it with zk1.

  1. From Lemma 21.5 we have

    xzk2=xb|K22xb2.
  2. Applying Lemma 21.4 for the least squares estimation step gives us

    2xb22(1.112xTc2+1.0541e2)=2.224xTc2+2.1082e2.
  3. Lemma 21.3 (Support merger) tells us that

    xTc2dΩck12.
  4. This gives us

    xzk22.224dΩck12+2.1082e2.
  5. From identification step we have Lemma 21.2

    dΩck120.223dk12+2.331e2.
  6. This gives us

    xzk22.224(0.223dk12+2.331e2)+2.1082e20.5dk12+7.5e2=12xzk12+7.5e2.

    The constants have been simplified to make them look better.

  7. For the 2nd result in this theorem, we add up the error at each stage as

    (1+21+22++2(k1))7.5e227.5e2=15e2.
  8. At k=1 we have z0=0.

  9. This gives us the result

    xzk22kx2+15e2.

21.4.3. CoSaMP Analysis General Case#

Having completed the analysis for the sparse case (where the signal x is K-sparse) it is time for us to generalize the analysis for the case where x is assumed to be an arbitrary signal in CN. Although it may look hard at first sight but there is a simple way to transform the problem into the problem of CoSaMP for the sparse case.

  1. We decompose x into its K-sparse approximation and the approximation error.

  2. Further, we absorb the approximation error term into the measurement error term.

  3. Since sparse case analysis is applicable for an arbitrary measurement error, this approach gives us an upper bound on the performance of CoSaMP over arbitrary signals.


  1. We start with writing

    x=xx|K+x|K

    where x|K is the best K-term approximation of x (Theorem 18.18).

  2. Thus we have

    y=Φx+e=Φ(xx|K+x|K)+e=Φx|K+Φ(xx|K)+e.
  3. We define

    e^=Φ(xx|K)+e.
  4. This lets us write

    y=Φx|K+e^.
  5. In this formulation, the problem is equivalent to recovering the K-sparse signal x|K from the measurement vector y.

  6. The results of CoSaMP Analysis Sparse Case and in particular the iteration invariant Theorem 21.4 apply directly.

  7. The remaining problem is to estimate the norm of modified error e^.

  8. We have

    e^2=Φ(xx|K)+e2Φ(xx|K)2+e2.
  9. Another result for RIP on the energy bound of embedding of arbitrary signals from Theorem 18.68 gives us

    Φ(xx|K)21+δK[xx|K2+1Kxx|K1].
  10. Thus we have an upper bound on e^2 given by

    e^21+δK[xx|K2+1Kxx|K1]+e2.
  11. Since we have assumed throughout that δ4K0.1, it gives us

    e^21.05[xx|K2+1Kxx|K1]+e2.
  12. This inequality is able to combine measurement error and approximation error into a single expression.

  13. To fix ideas further, we define the notion of unrecoverable energy in CoSaMP.

Definition 21.2 (Unrecoverable energy)

The unrecoverable energy in CoSaMP algorithm is defined as

ν=[xx|K2+1Kxx|K1]+e2.

This quantity measures the baseline error in the CoSaMP recovery consisting of measurement error and approximation error.

It is obvious that

e^21.05ν.

We summarize this analysis in following lemma.

Lemma 21.6

Let xCN be an arbitrary signal. Let x|K be its best K-term approximation. The measurement vector y=Φx+e can be expressed as y=Φx|K+e^ where

e^21.05[xx|K2+1Kxx|K1]+e21.05ν.

We now move ahead with the development of the iteration invariant for the general case.

  1. Invoking Theorem 21.4 for the iteration invariant for recovery of a sparse signal gives us

    x|Kzk212x|Kzk12+7.5e^2.
  2. What remains is to generalize this inequality for the arbitrary signal x itself.

  3. We can write

    x|Kx+xzk212x|Kx+xzk12+7.5e^2.
  4. We simplify L.H.S. as

    x|Kx+xzk2=(xzk)(xx|K)2xzk2xx|K2.
  5. On the R.H.S. we expand as

    x|Kx+xzk12xzk12+xx|K2.
  6. Combining the two we get

    xzk20.5xzk12+1.5xx|K2+7.5e^2.
  7. Putting the estimate of e^2 from Lemma 21.6 we get

    xzk20.5xzk12+9.375xx|K2+7.875Kxx|K1+7.5e2.
  8. Now

    9.375xx|K2+7.875Kxx|K1+7.5e210([xx|K2+1Kxx|K1]+e2).
  9. Thus we write a simplified expression

    xzk20.5xzk12+10ν

    where ν is the unrecoverable energy (Definition 21.2).

We can summarize the analysis for the general case in the following theorem.

Theorem 21.5 (CoSaMP iteration invariant for general signal recovery)

Assume that Φ satisfies RIP of order 4K with δ4K0.1. For each iteration number k1, the signal estimate zk is K-sparse and satisfies

xzk212xzk12+10ν.

In particular

(21.29)#xzk22kx2+20ν.

Proof. 1st result was developed in the development before the theorem. For the 2nd result in this theorem, we add up the error at each stage as

(1+21+22++2(k1))10ν210ν=20ν.

At k=1 we have z0=0. This gives us the result.

21.4.3.1. SNR Analysis#

  1. A sensible though unusual definition of signal-to-noise ratio (as proposed in [61]) is as follows

    (21.30)#SNR=10log(x2ν).
  2. The whole of unrecoverable energy is treated as noise.

  3. The signal 2 norm rather than its square is being treated as the measure of its energy. This is the unusual part.

  4. Yet the way ν has been developed, this definition is quite sensible.

  5. Further we define the reconstruction SNR or recovery SNR as the ratio between signal energy and recovery error energy.

    (21.31)#R-SNR=10log(x2xz2).
  6. Both SNR and R-SNR are expressed in dB.

  7. Certainly we have

    R-SNRSNR.
  8. Let us look closely at the iteration invariant

    xzk22kx2+20ν.
  9. In the initial iterations, 2kx2 term dominates in the R.H.S.

  10. Assuming

    2kx220ν

    we can write

    xzk222kx2.
  11. This gives us

    xzk22k+1x2x2xz22k1R-SNR10(k1)log23k3.
  12. In the later iterations, the 20ν term dominates in the R.H.S.

  13. Assuming

    2kx220ν

    we can write

    xzk2220ν=40ν.
  14. This gives us

    xzk240νx2xz2140x2νR-SNRSNR10log40SNR16=SNR133.
  15. We combine these two results into the following

    (21.32)#R-SNRmin{3k,SNR13}3.
  16. This result tells us that in the initial iterations the reconstruction SNR keeps improving by 3dB per iteration till it hits the noise floor given by SNR16 dB.

  17. Thus roughly the number of iterations required for converging to the noise floor is given by

    kSNR133.