PaperlayerView exact revision on arXiv

Source-linked AI summary

CoSaMP: Iterative signal recovery from incomplete and inaccurate samples

D. Needell, J. A. Tropp

arXiv:0803.2392v2math.NAcs.IT

TL;DR

Compressive sampling needs recovery methods that handle varied, noisy samples while using few samples, offering optimal error guarantees, and controlling computational resources. CoSaMP is an iterative reconstruction algorithm that achieves these goals with rigorous runtime and storage bounds, relying on matrix–vector multiplies and often running in O(N*log^2(N)).

  • Problem

    Existing recovery approaches do not simultaneously support varied sampling schemes, few samples, noise robustness, optimal error guarantees, and provably efficient resource use.

  • Method

    CoSaMP iteratively forms residual proxies to identify large components, estimates their values by least squares, and updates the approximation.

  • Results

    CoSaMP provides recovery guarantees for noisy measurements under a restricted-isometry condition, with O(L · log(∥x∥2 /η)) runtime and O(N) working storage.

  • Takeaways & Limitations

    The algorithm is designed for practical efficiency because it requires only matrix–vector multiplies with the sampling matrix and often runs in O(N*log^2(N)).

  • Takeaways & Limitations

    The appropriate precision for solving CoSaMP’s least-squares subproblems depends on relative computational costs and is ultimately best settled empirically.

Abstract

from arXiv · show

Compressive sampling offers a new paradigm for acquiring signals that are compressible with respect to an orthonormal basis. The major algorithmic challenge in compressive sampling is to approximate a compressible signal from noisy samples. This paper describes a new iterative recovery algorithm called CoSaMP that delivers the same guarantees as the best optimization-based approaches. Moreover, this algorithm offers rigorous bounds on computational cost and storage. It is likely to be extremely efficient for practical problems because it requires only matrix-vector multiplies with the sampling matrix. For many cases of interest, the running time is just O(N*log^2(N)), where N is the length of the signal.

1. Introduction

Compressive sampling seeks to recover information-rich signals from a small number of nonadaptive measurements rather than complete representations. The paper introduces CoSaMP, an iterative reconstruction algorithm designed to combine robust, accurate recovery with provably efficient computation and storage.

  • Motivation: Many signals are compressible, so complete acquisition wastes measurements and retains information that is later discarded.Compressive sampling instead uses a small number of nonadaptive samples that can approximate suitable signals well.
  • Motivation: A practical reconstruction method should support varied sampling schemes, use few samples, tolerate noise, provide optimal error guarantees, and use resources efficiently.The paper notes that no prior approach simultaneously achieves all five goals.
  • CoSaMP: CoSaMP is a novel greedy pursuit that combines orthogonal matching pursuit with combinatorial ideas to accelerate recovery and establish guarantees unavailable for OMP.Its analysis is also inspired by ROMP and convex-relaxation methods.
  • Guarantees: 2s-sparse approximation: under δ2s ≤ c, CoSaMP recovers an arbitrary signal from noisy samples using a sampling matrix with the stated restricted-isometry condition.The theorem allows arbitrary additive noise and accepts samples of the form u = Φx + e.
  • Efficiency: O(L · log(∥x∥2 /η)): the running time depends on matrix–vector multiplication cost, while working storage is O(N).For partial Fourier sampling, the total runtime is O(N log N · |R-SNR|), which is nearly linear in signal length for most signals of interest.

2. The CoSaMP Algorithm

CoSaMP iteratively identifies large residual components using signal proxies and provides explicit guarantees on recovery, iteration count, computational cost, and storage. Under exact arithmetic, it reaches an s-sparse approximation within at most 6(s + 1) iterations, while each iteration costs O(L) time and total storage is O(N).

  • Algorithm: CoSaMP forms a proxy from the sampling matrix and its adjoint to identify large signal or residual components, then uses the resulting support for the next approximation.The proxy relies on the restricted isometry property, and samples are updated to reflect the current residual at each iteration.
  • Performance guarantees: Each iteration makes significant progress when approximation error exceeds the unrecoverable energy, with reconstruction SNR improving by about 3 decibels per iteration until the noise floor.The stated bound is R-SNR > 3 −min{3k, SNR −13}.
  • Iteration count: 6(s + 1) iterations suffice under exact arithmetic to produce an s-sparse approximation meeting the theorem’s error guarantee.For signals with comparable largest entries, the iteration count may be as small as O(log s), whereas exponentially decaying entries can require Ω(s) iterations.
  • Resource requirements: O(L) time is required per iteration, where L bounds a multiplication with Φ or Φ∗, and the algorithm uses O(N) storage.When both Φ and Φ∗ have fast multiplies, a typical value is L = O(N log N), including for partial Fourier matrices.
  • Resource requirements: The total runtime is O(L log(∥x∥2 /η)) for the stated precision parameter η, while generic p-compressible signals require O(log s) iterations.The constant in the O(log s) iteration bound depends on p.

3. Restricted Isometry Consequences

The restricted isometry property yields norm bounds on submatrices, approximate orthogonality between disjoint column sets, control of higher restricted isometry constants, and an energy bound for nonsparse signals.

  • Norm bounds: For any index set T of at most r elements, the restricted isometry inequalities bound the relevant norms by 1 ± δr.These bounds follow because the singular values of Φ_T lie between √(1 − δr) and √(1 + δr).
  • Approximate orthogonality: Disjoint column sets S and T with combined cardinality at most r span nearly orthogonal subspaces, with deviation controlled by δr.The spectral-norm estimate uses the eigenvalue bounds implied by the restricted isometry property.
  • Higher restricted isometry constants: For positive integers c and r, the higher restricted isometry constants satisfy δcr ≤ c · δ2r.The proof partitions the associated matrix into r × r blocks and applies block Gershgorin bounds, using δr for diagonal blocks and δ2r for off-diagonal blocks.
  • Energy bound: The energy bound extends the sampling-matrix control from sparse signals to every signal, enabling sparse-signal results to transfer to the general case.The proposition assumes the upper restricted isometry inequality, namely 1 + δr for vectors with ∥x∥0 ≤ r.

4. The Iteration Invariant: Sparse Case

For s-sparse signals, each CoSaMP iteration preserves s-sparsity and reduces approximation error by a constant factor while adding a small multiple of the noise. The proof proceeds through identification, support merger, least-squares estimation, and pruning, then solves the resulting error recursion.

  • Iteration invariant: Each iteration reduces approximation error by a constant factor while adding a small multiple of the noise.When approximation error is large relative to noise, the algorithm makes substantial progress identifying the unknown signal.
  • Identification: The identification phase selects at most 2s residual components with the largest proxy magnitudes.The selected set contains components where the residual signal has substantial energy.
  • Support merger: Merging the identified set with the current approximation support produces a support set containing at most 3s indices.The support-merger bound shows that signal components outside this set have very little energy.
  • Estimation: The estimation phase computes coefficient values on the merged support by solving a least-squares problem.The stated bound assumes infinite-precision least squares; iterative solvers add an extra error term, analyzed in Section 5.
  • Error recursion: 15 ∥e∥2 bounds the accumulated noise contribution after solving the error recursion.The recursion uses the geometric factor (1 + 0.5 + 0.25 + . . . ) · 7.5 ∥e∥2.

5. Analysis of Iterative Least-squares

The analysis shows that CoSaMP can solve its least-squares subproblems with iterative methods whose performance follows from restricted isometry. Initializing with the current approximation and using three iterations preserves the sparse iteration invariant.

  • Conditions and convergence: The condition δ4s ≤0.1 ensures a well-conditioned least-squares problem because the merged support contains at most 3s components.This condition underlies the performance of Richardson’s iteration and conjugate gradient.
  • Conditions and convergence: Richardson’s iteration computes the least-squares solution using only matrix–vector multiplies with A and A∗ and converges linearly.For A = ΦT with |T| ≤3s, δ3s ≤δ4s ≤0.1 implies rapid convergence.
  • Initialization: CoSaMP naturally initializes iterative least squares with the current signal approximation, which becomes an increasingly good starting point as iterations progress.The initial iterate’s error is controlled by the current approximation error.
  • Iteration count: After at most three iterations, Richardson’s iteration and conjugate gradient produce a sufficiently accurate signal estimate when initialized with z0 = ak−1.The conjugate-gradient bound is slightly better.
  • Sparse iteration invariant: Three least-squares iterations preserve Theorem 4.1’s sparse iteration invariant for either Richardson’s iteration or conjugate gradient.The theorem uses the current approximation ak−1 as the initializer.

6. Extension to General Signals

The section extends CoSaMP’s guarantees from sparse to arbitrary signals by absorbing the signal tail into an effective noise vector. Decomposing the signal enables the sparse-case analysis to yield a final error bound for general signals.

  • Reduction to Sparse Case: A general signal’s noisy samples can be rewritten as samples of its best sparse approximation plus an effective noise vector incorporating the approximation tail.This reduction follows by decomposing x = x_s + (x−x_s), so e = Φ(x−x_s)+e.
  • Reduction to Sparse Case: √1 + δ_s ≤1.05 follows from δ_s ≤δ_4s ≤0.1 in the reduction’s norm estimate.This numerical simplification is used when controlling the effective noise term.
  • Proof of Theorem 2.1: The proof completes Theorem 2.1 by applying the sparse iteration invariant or its iterative least-squares analogue to the reduced representation.The argument then uses lower and upper triangle inequalities before substituting the effective-noise estimate.
  • Final Bound: 2 + 9.375 ∥x −x_s∥2 + 7.875 √s ∥x −x_s∥1 + 7.5 ∥e∥2 is the resulting error expression.The bound combines the sparse-approximation tail in ℓ2 and ℓ1 with the original noise magnitude.

7. Discussion and Related Work

CoSaMP is compared with convex-relaxation, greedy, and sublinear recovery methods, combining comparable recovery guarantees with favorable computational and sampling properties. The discussion also situates CoSaMP’s multiple-component selection strategy within earlier work.

  • Convex relaxation and greedy methods: CoSaMP’s error bound is equivalent to convex-relaxation guarantees, modulo the exact value of the constants.Convex relaxation requires restricted-isometry hypotheses, while CoSaMP avoids the extra logarithmic factors identified for ROMP and is described as essentially optimal.
  • Convex relaxation and greedy methods: OMP reconstructs a fixed s-sparse signal after s iterations with m = Cs log N measurements, except with probability N−1.The result uses an m × N matrix with independent, zero-mean subgaussian entries, but OMP lacks uniform guarantees and fails for some sparse signals.
  • Sublinear algorithms: Sublinear methods require more measurements by logarithmic factors and highly structured samples, which may limit their practical usefulness.HHS pursuit achieves the error bound (7.2), while chaining pursuit has the weaker bound ∥x −a∥1 ≤C log N ∥x −xs∥1.
  • Algorithm comparison: CoSaMP achieves the best performance among linear and superlinear algorithms across the comparison criteria in Table 2.The criteria include general samples, optimal sample count, uniformity, stability, and running time; interior-point methods cost O(m2N1.5).
  • Algorithmic origins: CoSaMP selects many components per iteration, a strategy incorporated in several earlier algorithms and essential to Fourier sampling.The discussion notes that the first appearance of this idea is unclear, while identifying earlier work by Gilbert et al. and later compressive-sampling methods.

Appendix A. Algorithmic Variations

The appendix discusses algorithmic variations, focusing on three natural halting criteria and their consequences. These include fixed iteration counts and norm-based evidence from current samples or the proxy.

  • Halting criteria: The appendix identifies three natural approaches to halting the algorithm.They are fixed iterations, current-sample norms, and proxy-entry magnitudes.
  • Halting criteria: The first criterion stops the algorithm after a fixed number of iterations.This is the approach discussed in the paper’s body.
  • Halting criteria: A second criterion uses the current sample norm ∥v∥2 as evidence about the residual norm ∥r∥2.The criterion relates the norm of the current samples to the norm of the residual.
  • Halting criteria: A third criterion uses the proxy-entry magnitude ∥y∥∞ to bound the residual-entry magnitude ∥r∥∞.This approach compares the maximum proxy-entry magnitude with the maximum residual-entry magnitude.

A.1. Halting Rules.

The section develops halting rules for sparse-signal recovery, first relating an a priori residual threshold to final approximation quality and then showing when the criteria are triggered. The guarantees rely on Proposition 3.1 and bounds such as δ2s ≤ δ4s ≤ 0.1.

  • Setup: For sparse signals, the residual framework defines r = x − a, noisy residual samples v = Φr + e, and proxy y = Φ∗v.The general case is reduced to sampling a sparse signal, with x s-sparse and a an s-sparse approximation.
  • Halting II: Halting II states that ∥v∥2 ≤ ε is triggered once the residual satisfies the theorem’s specified condition.The argument uses Proposition 3.1 and δ2s ≤ 0.1.

A.2. Other Variations.

The section presents three CoSaMP variations: one updates the residual, another adds a final least-squares refinement, and a third prunes the merged support before solving. The residual-based version has a theorem analogous to Theorem 2.1, but no verified analog of Theorem 2.2.

  • (1) Residual-based update: The first variation approximates the current residual at each iteration rather than the entire signal, similarly to HHS Pursuit.Its estimation step solves least squares using the current samples instead of the original samples.
  • (1) Residual-based update: The residual approximation is merged with the previous signal approximation, then pruned to an s-sparse signal.The merger sets c = a_{k−1} + b, and pruning sets a_k = c_s.
  • (1) Residual-based update: This residual-based algorithm satisfies a result similar to Theorem 2.1, but the authors could not verify an analog of Theorem 2.2.The authors nevertheless describe this version as promising.
  • (2) Final least-squares refinement: A second variation solves another least-squares problem on the final support and outputs an s-sparse approximation, but noise means it is not guaranteed to improve the result.The output should never be much worse than the end-of-loop approximation.
  • (3) Pre-solve support pruning: A third variation prunes the merged support to s entries before least squares; its analysis was not attempted.Proxy values may serve as surrogates for unknown new-approximation values, and the identification step then selects new components each iteration.

Appendix B. Iteration Count · B.1. Additional Notation.

Appendix B shows that CoSaMP’s iteration count depends on signal structure, quantified by the signal profile, and gives a sparse-case bound. It introduces auxiliary band-tail signals and an iteration-by-iteration alternative used in the proof.

  • Appendix B. Iteration Count: CoSaMP’s iteration count depends heavily on the signal structure, with rapidly decaying entries potentially requiring an iteration or more per component.A large residual component contaminates every component of the signal proxy, delaying progress on smaller entries.
  • Appendix B. Iteration Count: Comparable signal entries can be located simultaneously, allowing CoSaMP to find all s components after about log s iterations.This follows by reducing the residual norm below the smallest-entry magnitude.
  • Appendix B. Iteration Count: The signal profile counts the nonempty component bands, representing the number of coefficient orders of magnitude and satisfying profile(x) ≤ s for an s-sparse signal.Component bands group signal components with comparable magnitudes.
  • Appendix B. Iteration Count: For an s-sparse signal with p = profile(x), Theorem B.1 bounds the iterations needed for CoSaMP to produce an approximation satisfying the stated guarantee.The supplied passage introduces the theorem and its approximation guarantee, but does not include the displayed formulas.
  • Appendix B. Iteration Count: 6(s + 1) is an absolute upper bound on the iteration count, since the fixed-s maximum occurs at p = s and log4/3 5.6 < 6.The theorem’s iteration bound therefore never exceeds 6(s + 1).
  • Appendix B. Iteration Count: The proof uses ak for the kth approximation, Sk = supp(ak) for its support, and rk = x −ak for the residual, whose norm represents approximation error.The notation is introduced for k = 0, 1, 2, . . . .
  • B.1. Additional Notation.: For nonnegative integer j, the auxiliary signal yj collects the portion of x in bands Bj, Bj+1, Bj+2, . . . .These auxiliary signals play a key role in the analysis.
  • B.1. Additional Notation.: Each iteration’s proof establishes an alternative: either the approximation error is small, or it is dominated by unidentified signal energy and the next error is a constant factor smaller.The first case means the algorithm is effectively finished.

B.2. Proof of Theorem B.1.

The proof establishes a two-alternative iteration invariant, then shows that one alternative identifies the full signal support within a bounded number of iterations. It finally uses this support result and the iteration invariant to establish the theorem’s approximation guarantee.

  • B.2. Proof of Theorem B.1.: Lemma B.2 shows that every iteration satisfies at least one of two alternatives governing residual reduction and approximation behavior.The proof derives the alternatives from the estimation, pruning, and sparse iteration results.
  • B.2. Proof of Theorem B.1.: The second alternative implies complete support recovery: once a band enters the support, residual decay makes it persist in all subsequent approximations.The argument proceeds band by band, using the residual threshold condition and the monotonic residual decline from (B.5).
  • B.2. Proof of Theorem B.1.: β−1 = 0.75 reduces the error by that factor per iteration, yielding an upper bound on the iterations needed to identify each successive nonempty band.The proof sums the per-band iteration requirements to bound the total support-identification count.
  • B.2. Proof of Theorem B.1.: For each iteration k ≥⌊k⋆⌋, the approximation support equals the signal support, and k⋆ is bounded in terms of the signal profile p.The proof uses p = |J|, the geometric mean–arithmetic mean inequality, Jensen’s inequality, and the fact that all bands contain s elements in total.
  • B.2. Proof of Theorem B.1.: If the first alternative never holds during the first K iterations, Lemma B.3 gives exact support recovery by iteration K, and the resulting error estimate contradicts that alternative’s failure.Therefore, some iteration k ≤K satisfies (B.3), after which repeated applications of Theorem 4.1 complete the argument.

B.3. Proof of Theorem 2.2.

The proof bounds CoSaMP’s iteration count for general signals by invoking the sparse-case result and optimizing the profile parameter. This yields a signal-independent upper bound of 6(s + 1) iterations.

  • Iteration-count bound: Theorem B.4 states that CoSaMP’s iteration count depends on the signal profile p = profile(x_s).The theorem considers an arbitrary signal x and guarantees an approximation after a profile-dependent number of iterations.
  • Proof of Theorem 2.2: The proof represents noisy samples as u = Φx_s + e and applies the sparse iteration-count theorem.The sparse result supplies the iteration-dependent approximation guarantee, which is then combined with the noise estimate.
  • Proof of Theorem 2.2: 6(s + 1) iterations suffice as a signal-independent upper bound for CoSaMP.The bound follows because the iteration estimate is maximized at p = s.
Loading 0803.2392v2…