PaperlayerView exact revision on arXiv

Source-linked AI summary

Message Passing Algorithms for Compressed Sensing

David L. Donoho, Arian Maleki, Andrea Montanari

arXiv:0907.3574v1cs.ITcond-mat.dis-nnstat.CO

TL;DR

The paper addresses the expensive computation of LP-based compressed-sensing reconstruction and the poorer tradeoffs of existing fast iterative thresholding. It introduces a belief-propagation-inspired modification and finds LP-matching sparsity–undersampling performance, with state evolution accurately describing AMP and agreeing with prior geometric results.

  • Problem

    LP reconstruction offers strong sparsity–undersampling tradeoffs but is expensive for large applications, while existing fast iterative thresholding methods fall short of those tradeoffs.

  • Method

    The paper modifies iterative thresholding with a message-passing algorithm inspired by belief propagation and analyzes it using state evolution.

  • Results

    AMP matches the theoretical sparsity–undersampling tradeoff proved for LP-based reconstruction, while state evolution accurately predicts its behavior and failure region.

  • Takeaways & Limitations

    A fast iterative algorithm can achieve reconstruction performance equivalent to corresponding LP methods on random problems, supported by simulations and theory.

Abstract

from arXiv · show

Compressed sensing aims to undersample certain high-dimensional signals, yet accurately reconstruct them by exploiting signal characteristics. Accurate reconstruction is possible when the object to be recovered is sufficiently sparse in a known basis. Currently, the best known sparsity-undersampling tradeoff is achieved when reconstructing by convex optimization -- which is expensive in important large-scale applications. Fast iterative thresholding algorithms have been intensively studied as alternatives to convex optimization for large-scale problems. Unfortunately known fast algorithms offer substantially worse sparsity-undersampling tradeoffs than convex optimization. We introduce a simple costless modification to iterative thresholding making the sparsity-undersampling tradeoff of the new algorithms equivalent to that of the corresponding convex optimization procedures. The new iterative-thresholding algorithms are inspired by belief propagation in graphical models. Our empirical measurements of the sparsity-undersampling tradeoff for the new algorithms agree with theoretical calculations. We show that a state evolution formalism correctly derives the true sparsity-undersampling tradeoff. There is a surprising agreement between earlier calculations based on random convex polytopes and this new, apparently very different theoretical formalism.

I. INTRODUCTION AND OVERVIEW

The paper develops a fast message-passing modification of iterative thresholding that matches LP reconstruction’s sparsity–undersampling tradeoff, supported by simulations and state-evolution analysis.

  • Motivation: Compressed sensing recovers high-dimensional signals from fewer measurements by exploiting sparse structure, but LP-based reconstruction is expensive at large scale.The motivating example involves an image with 10^6 pixels and an LP with tens of thousands of constraints and millions of variables.
  • Method: The proposed message-passing algorithm iteratively updates the current estimate and residual using scalar threshold functions and a transpose measurement operator.The algorithm is inspired by belief propagation in graphical models and adds a term absent from standard iterative thresholding.
  • Motivation: Fast iterative thresholding has low per-iteration cost and storage requirements but previously offered a worse sparsity–undersampling tradeoff than LP reconstruction.The paper targets this gap without abandoning the scalability advantages of iterative methods.
  • Results: AMP achieves a sparsity–undersampling tradeoff matching the theoretical tradeoff proved for LP-based reconstruction, with both methods succeeding or failing in the same parameter-space regions.The comparison uses sparsity and undersampling axes in the large-dimensional limit.
  • Results: State evolution models reconstruction mean squared error through a scalar MSE map whose fixed points predict AMP failure regions.The analytically evaluated region coincides within numerical precision with the LP failure region, while Monte Carlo results agree within statistical precision.
  • Conclusion: The paper reports a fast iterative algorithm with LP-equivalent reconstruction performance on random problems, supported by simulations and a theoretical formalism.It also connects the LP success/failure phases to a simple phase-boundary formula derived from a different theoretical principle.

B. Phase Transitions

The paper characterizes compressed-sensing phase transitions using sparsity and undersampling ratios, and relates iterative reconstruction behavior to signal structure and interference suppression.

  • Phase-transition framework: In the large-system limit, (δ, ρ) partitions into success and failure phases, with δ measuring undersampling and ρ measuring sparsity.The phase boundary describes where exact reconstruction typically succeeds or fails.
  • Signal structure: Knowing that a signal is sparse and positive improves the phase transition over knowing only that it is sparse.The ordering ρCG(δ; +) > ρCG(δ; ±) expresses this advantage.
  • Signal structure: For sparse nonnegative and signed vectors, ρCG(δ; +, ±) asymptotically behaves as (2 log(1/δ))^-1 as δ → 0.This indicates that substantial undersampling is possible when sparsity is sufficient.
  • Signal structure: Vectors bounded in [−1, 1] have ρCG(δ; □) = 0 for δ < 1/2, so these bounds provide no phase-transition benefit in that regime.The bounded setting becomes useful only with a limited amount of undersampling.
  • Analytical phase boundaries: The phase-boundary formula agrees numerically with earlier combinatorial-geometry expressions while using methods unrelated to combinatorial geometry.The paper presents this agreement as a principal result.
  • Iterative reconstruction: The iteration uses thresholding controlled by the current MSE and can operate using applications of A and A∗ without solving linear systems.This supports rapid application to operators such as Fourier transforms.

E. State Evolution

State Evolution models AMP’s mean squared error through a scalar recursive MSE map and uses its fixed points to identify reconstruction phase transitions.

  • MSE recursion: State Evolution recursively tracks formal MSE through a scalar MSE map Ψ that represents iteration-to-iteration error change.The map is defined as a one-dimensional dynamical system with fixed parameters.
  • MSE recursion: The parameter space divides according to whether Ψ decreases MSE toward zero or fails to reach zero.Region (I) has Ψ(σ2) < σ2 throughout the relevant positive-MSE range; Region (II) does not evolve to zero.
  • Phase threshold: The State Evolution sparsity threshold ρSE is the largest ρ for which the parameter setting lies in the convergence region.This threshold depends on the signal setting χ and undersampling, with λ and signal-distribution notation introduced in the formalism.
  • Phase threshold: For the three canonical problems, ρSE(δ; χ, λ, FX) is independent of the signal distribution FX when sparsity and bounded second moment are prescribed.This permits writing the threshold without the distribution argument.
  • Main correspondence: High-precision evaluations yield the paper’s Finding 1: formal MSE reaches zero over the same (δ, ρ) region as corresponding convex optimization.The result connects State Evolution’s phase boundary with the convex-optimization phase diagram.

F. Failure of standard iterative algorithms

Standard iterative thresholding does not realize the favorable formal phase diagram, but a message-passing modification makes State Evolution accurate and achieves LP-like asymptotic performance.

  • Standard iteration failure: Formal MSE calculations do not describe standard iterative thresholding, whose algorithms fail to achieve the optimal phase diagram even when optimally tuned.Earlier empirical studies found transitions at substantially smaller sparsity values than LP-based methods.
  • Standard iteration failure: The failure arises because standard iterations create strong dependencies among estimate entries after the first iteration, invalidating State Evolution’s moment predictions.Numerical simulations directly show that the MSE map does not track the actual MSE.
  • Message passing: The paper modifies iterative thresholding using message-passing ideas from graphical-model inference and error-correcting codes.The modification associates variables with messages on a complete bipartite graph.
  • Message passing: Message passing requires updating Nn messages rather than N estimates, increasing complexity, while neglecting weak index dependence recovers iterative-thresholding equations.This motivates the approximate message-passing formulation.
  • State Evolution validation: Properly tuned MP-based algorithms are asymptotically as powerful as LP reconstruction, despite standard iterative thresholding’s weaker phase diagram.The algorithm’s tuning parameter can be selected using State Evolution.
  • State Evolution validation: For AMP, State Evolution accurately predicts MSE, false alarms, missed detections, and convergence behavior in large dimensions.Empirical quantities agree with the formal predictions, including convergence below and nonconvergence above ρSE.

II. DETAILS ABOUT THE MSE MAPPING

The section derives the state-evolution threshold from the MSE map, using local stability and concavity to establish convergence and distribution-independence. It also obtains explicit optimal tuning parameters.

  • MSE threshold and local stability: The state-evolution threshold is defined by requiring the MSE map Ψ(σ2) to remain below the σ2 line for all σ2 > 0.Because Ψ(0) = 0, local stability requires the derivative at σ2 = 0 to be at most 1.
  • Distribution-independence: For χ = +, the local stability threshold depends on the signal distribution only through sparsity and is independent of FX.The same independence is established for ρSE(δ; χ).
  • Convergence: Concavity of the MSE map makes the local stability condition sufficient for convergence to zero error.The resulting convergence to the correct solution is exponential.
  • Concavity proof: The proof of concavity reduces the analysis to deterministic nonzero signals because convex combinations of concave functions remain concave.The contribution from X = 0 vanishes.
  • Optimal tuning: The argument also yields explicit optimal tuning parameters by maximizing the local stability threshold.Practical use of the formula is accompanied by a notice in the Supplemental Information.

III. DISCUSSION

The discussion characterizes when state evolution is distribution-insensitive and examines the scope of possible improvements, matrix ensembles, and large-operator implementations.

  • State-evolution scope: For the threshold functions studied here, monotonicity and concavity make the optimally tuned AMP phase transition independent of the empirical distribution of x0.The discussion cautions that state evolution may be inaccurate without these properties.
  • Limitations: The improvement limit is conditional on accurate state evolution, and the comparison is stated asymptotically as δ → 0.The bound is ˜ρSE(δ; χ) ≤ ρSE(δ; χ)(1 + o(1)).
  • Possible improvements: In the high-undersampling limit, alternative nonlinearities offer essentially unimprovable state-evolution phase transitions.Reconstruction experiments also suggest little improvement from other nonlinearities.
  • Signal distributions: Theoretical phase transitions are insensitive to coefficient-amplitude distributions, and empirical MP results show similar insensitivity.This supports robustness across signal-amplitude distributions within the studied setting.
  • Measurement matrices: Comparable phase transitions occur for random matrix ensembles beyond Gaussian matrices.The evidence comes from Supporting Information experiments.
  • Large-scale applications: Partial Fourier operators produce observed MP phase transitions comparable to those for random A.This matters for applications where very large matrices are applied as operators rather than explicitly represented.

APPENDIX

The appendix supplies thresholding formulas, derivative-based local-stability calculations, concavity arguments, and the resulting thresholds for the +, ±, and □ settings.

  • Parameter conventions: The appendix notes that its threshold parameter λ is implicitly rescaled relative to traditional standard-deviation units by a δ-dependent factor.The conversion is introduced through λsd and the appropriately rescaled λSE.
  • χ = +: For χ = +, the threshold function sets η(x; λσ) to x − λσ above threshold and zero otherwise.Its derivative is then used to obtain the local stability threshold.
  • χ = ±: For χ = ±, soft thresholding retains x − λσ or x + λσ beyond the positive or negative threshold and returns zero otherwise.The derivatives are expressed using indicators of |x| ≥ λσ.
  • χ = □: For χ = □, the proposed nonlinearity clips inputs to the interval [−1, +1] and has no threshold parameter.The local stability threshold is ρLS(δ; □) = (2 − δ−1)+.

D. Relation to Minimax Thresholding ǫ

The minimax-thresholding analysis chooses thresholds using sparsity and effective noise rather than the unknown signal distribution, yielding guaranteed state-evolution convergence below a matching phase transition.

  • Minimax setup: The minimax threshold is defined over sparse distributions supported on nonnegative or unrestricted real values, with η chosen according to χ.The corresponding minimax threshold MSE uses expectation under each distribution class.
  • Threshold selection: At each iteration, the minimax threshold is scaled by the effective noise level and does not depend on the unknown distribution F.The resulting formal MSE guarantee is the minimax threshold risk scaled by the squared effective noise level.
  • MSE reduction: State evolution with minimax thresholding guarantees a fractional MSE reduction at each iteration when the minimax-risk bound is below the current σ2.The formal evolution is consequently bounded by the stated risk-based recursion.
  • Convergence guarantee: Below ρMM(δ; χ), state evolution with minimax thresholding converges to zero formal MSE for every marginal distribution F in the relevant class.The theorem also bounds the dynamics by the preceding risk recursions.
  • Relation to optimal tuning: The minimax threshold requires δ and ρ, whereas the main-text optimal threshold depends only on δ; setting ρ at the minimax transition removes that requirement.The minimax and state-evolution transitions coincide through M∗(ρδ; χ) = δ if and only if ρ = ρSE(δ; χ).
  • Derivation: The derivation uses scale invariance of η and concavity of the MSE map to relate the minimax threshold risk to the state-evolution threshold.The signal distribution is rescaled through the operator SaF.

E. Convergence Rate of State Evolution

The section characterizes AMP convergence through state evolution and relates its phase transition to combinatorial-geometry calculations. It also contrasts optimal and minimax thresholding through their convergence behavior.

  • Threshold policies: Optimal thresholding can differ from minimax thresholding because minimax adapts to the actual sparsity ratio, while the optimal policy assumes the largest plausible ratio.Below the state-evolution phase transition, both methods converge, but the ρ-adaptive minimax threshold generally provides better speed guarantees.
  • Convergence rate: The AMP convergence rate is expressed through the state-evolution factor ωSE(δ, ρ; χ), with formal MSE bounded by ωSE(δ, ρ; χ)^t EX2.The equivalent exponential rate parameter is b(δ, ρ) = −log(ωSE(δ, ρ)).
  • Phase transitions: State evolution and combinatorial geometry produce phase transitions ρSE and ρCG that are rigorously equivalent in the highly undersampled limit δ → 0.The paper notes that their agreement was observed numerically in the main text and established rigorously in this limit.
  • Asymptotics: The extremal distribution attaining the supremum is a two-point mixture for χ = + and a three-point symmetric mixture for χ = ±.These distributions determine the worst-case thresholding-risk calculations used in the asymptotic analysis.
  • Asymptotics: For signed and nonnegative sparse signals, the asymptotic minimax quantity satisfies M∗(ϵ; χ) = 2 log(ϵ)(1 + o(1)) as ϵ → 0.The section states that corresponding precise asymptotic expressions for µχ(ϵ) are available.

G. Rigorous Asymptotic Optimality of Soft Thresholding

This section analyzes whether alternative nonlinearities can improve soft thresholding under state evolution. It concludes that soft thresholding is asymptotically minimax-optimal in the extreme-undersampling regime.

  • General nonlinearities: Under state evolution, choosing a minimax threshold for a general nonlinearity is guaranteed to reduce MSE when the stated condition holds.When the condition is satisfied, the formal MSE evolves fully to zero.
  • General nonlinearities: The generalized minimax phase transition eρSE(δ; χ) is defined by the largest ρ for which state evolution drives the formal MSE to zero.For every signal distribution F and ρ below eρSE(δ), the formal MSE evolves to zero.
  • Optimality result: The analysis asks whether a specially chosen nonlinearity can substantially outperform soft thresholding in minimax phase transitions.The comparison is made for χ ∈ {+, ±} under the state-evolution formalism.
  • Optimality result: No other nonlinearity can outperform optimally tuned soft thresholding in minimax phase transitions for χ ∈ {+, ±} in the limit of extreme undersampling.The result is interpreted through the range of sparsity ratios guaranteeing the stable-zero-fixed-point region of the MSE map.
  • Proof strategy: The minimax-risk comparison permits optimization over all measurable scalar nonlinearities, making the optimality claim broader than a comparison among a few named thresholding rules.The proof uses the fact that any specific nonlinearity cannot do better than the minimax risk.

J. Estimating Phase Transitions

The paper estimates finite-dimensional phase transitions from Monte Carlo success data and tests state-evolution predictions against AMP observations. Its experiments compare multiple algorithms and track interference and reconstruction observables across iterations.

  • Estimating phase transitions: Finite-N phase transitions are defined as the sparsity ratio ρ where the reconstruction success probability reaches 50%.The estimate depends implicitly on undersampling δ, problem size N, and the success tolerance.
  • Estimating phase transitions: The experiments model success counts as binomial observations whose probability varies with ρ according to a logistic generalized linear model.For fixed δ, the data are represented by triples (k, M, S(k, n, N)) with ρ = k/n.
  • Estimating phase transitions: Phase-transition estimates are obtained for each algorithm-parameter and problem-suite combination, then optimized over parameters or across suites for minimax performance.The resulting observed minimax transitions are presented for several algorithms.
  • Empirical phase transitions: The standard-suite experiments use a uniform spherical measurement ensemble, sparse constant-amplitude coefficient vectors, and M = 20 independent instances per grid point.For AMP, the study uses a focused search around ρCG(δ) and runs T = 1000 iterations before reporting final-iteration mean squared error.
  • Empirical phase transitions: Figure 4 compares observed phase transitions for AMP, IST, IHT, TST, OMP, and ℓ1, with the ℓ1 curve coinciding with the theoretical state-evolution transition ρSE.The other curves represent empirical results.
  • Testing state evolution: QQ plots test whether mutual access interference is approximately Gaussian during AMP iterations and whether its standard deviation decreases over time.For the USE ensemble, approximate Gaussianity persists through 90 iterations; a partial-Fourier example is also examined.
  • Testing state evolution: State evolution predicts the iteration-by-iteration behavior of MSE, missed-detection rate, false-alarm rate, and zero/nonzero MSE observables.These predictions require the evolving σt together with ϵ and λ.
  • Testing state evolution: For ρ = .3 and δ = .15, theoretical state-evolution curves closely agree with empirical mean observables across all four panels.Figures 8 and 9 extend the same comparison to additional parameter settings.

O. Coefficient Universality

The study finds substantial empirical invariance across coefficient and matrix ensembles, while AMP retains favorable accuracy and scaling behavior in large problems.

  • Coefficient Universality: The four coefficient ensembles produce visually similar success-probability curves, indicating considerable invariance in AMP performance.The comparison uses unit-amplitude, uniform, Gaussian, and Cauchy nonzero coefficients at δ = 0.10 and δ = 0.3.
  • Matrix Ensembles: AMP phase transitions are compared across Uniform Spherical, Rademacher, and partial Fourier matrix ensembles against the ℓ1 transition.The study considers the χ = ± setting and reports results in Figure 11.
  • Timing Results: AMP is substantially faster than LARS in all studied situations, with speedups reaching factors of 10 or more as sparsity complexity increases.For very sparse cases, LARS loses by less than a factor of 3; when k reaches several hundred or one thousand, AMP’s advantage is at least tenfold.
  • Timing Results: For dense matrices, AMP requires at most c2·n·N operations while LARS requires at least c1·k·n·N operations.The denominator’s logarithmic dependence is weak, so the operation-count ratio scales roughly directly with k.
  • Scaling with Problem Size: AMP iteration counts remain nearly constant for large N at fixed (δ, ρ), while each additional five iterations approximately halves the actual MSE.This behavior agrees with the state-evolution prediction that iteration counts are independent of N and each iteration yields proportional MSE reduction.
  • Scaling with Problem Size: With partial Fourier operators, AMP timing appears practically linear in N, and each doubling of accuracy costs about 30% more computation time.Applying A and A∗ costs O(N log(N)), substantially less than applying a general dense matrix.
Loading 0907.3574v1…