PaperlayerView exact revision on arXiv

Source-linked AI summary

Expectation-Maximization Gaussian-Mixture Approximate Message Passing

Jeremy P. Vila, Philip Schniter

arXiv:1207.3107v3cs.IT

TL;DR

Sparse recovery can suffer when the nonzero-coefficient distribution and noise variance are unknown, despite MMSE methods benefiting from known prior statistics. The paper learns a Gaussian-mixture prior and noise variance with EM, using AMP to approximate the E step and recover the signal. Experiments report state-of-the-art reconstruction error and runtime across many signal and matrix types, with difficulties for certain sensing matrices.

  • Problem

    Unknown nonzero-coefficient distributions and noise variance limit direct use of efficient approximate-MMSE recovery, while LASSO cannot exploit easier signal distributions.

  • Method

    EM-GM-AMP models the active-coefficient distribution as a Gaussian mixture, learns its parameters and noise variance with EM, and uses AMP for approximate posterior inference and recovery.

  • Results

    The approach achieves a state-of-the-art combination of reconstruction error and runtime across a wide range of signal and matrix types in the high-dimensional regime.

  • Takeaways & Limitations

    Learning the prior can substantially improve recovery over competing methods for signal classes that the Gaussian-mixture model represents accurately.

  • Takeaways & Limitations

    Certain non-zero-mean and super-Gaussian sensing matrices cause difficulties for the AMP-based method.

Abstract

from arXiv · show

When recovering a sparse signal from noisy compressive linear measurements, the distribution of the signal's non-zero coefficients can have a profound effect on recovery mean-squared error (MSE). If this distribution was apriori known, then one could use computationally efficient approximate message passing (AMP) techniques for nearly minimum MSE (MMSE) recovery. In practice, though, the distribution is unknown, motivating the use of robust algorithms like LASSO---which is nearly minimax optimal---at the cost of significantly larger MSE for non-least-favorable distributions. As an alternative, we propose an empirical-Bayesian technique that simultaneously learns the signal distribution while MMSE-recovering the signal---according to the learned distribution---using AMP. In particular, we model the non-zero distribution as a Gaussian mixture, and learn its parameters through expectation maximization, using AMP to implement the expectation step. Numerical experiments on a wide range of signal classes confirm the state-of-the-art performance of our approach, in both reconstruction error and runtime, in the high-dimensional regime, for most (but not all) sensing operators.

I. INTRODUCTION

The paper targets sparse recovery when the nonzero-coefficient distribution is unknown, combining Gaussian-mixture modeling, EM learning, and AMP-based recovery. It aims to approach MMSE performance without known prior statistics while retaining computational efficiency.

  • Motivation: LASSO’s phase transition is invariant to the active-coefficient distribution, making it robust but unable to exploit easier signal classes.For nonnegative signals, polynomial-complexity algorithms can achieve better phase transitions than LASSO.
  • Motivation: Known signal and noise statistics enable MMSE recovery with a phase-transition boundary K/M = 1, motivating prior learning when those statistics are unavailable.Bayesian AMP provides an efficient approximate-MMSE option when the marginal pdf is known.
  • Proposed approach: EM-GM-AMP models active coefficients with an L-term Gaussian mixture and learns mixture parameters and AWGN variance while recovering the signal.GAMP supplies approximate-MMSE signal estimates and posterior activity probabilities for EM and support recovery.
  • Signal model: The signal prior combines a point mass at zero with Gaussian-mixture components, with λ denoting sparsity and ωk, θk, and φk denoting component weight, mean, and variance.The model assumes independent Gaussian noise with variance ψ; parameters and L are treated as fixed and known within GM-GAMP.
  • Computational framework: GAMP approximates coefficient and output posteriors, enabling closed-form computations for Gaussian-mixture priors and efficient EM updates.For arbitrary priors such as Student’s-t, numerical integration could substantially increase GAMP complexity.
  • Scope: The framework can extend to non-Gaussian noise and non-additive channels because GAMP supports arbitrary output-channel models.The paper notes examples including quantized outputs and logistic regression, while also describing mixture-based learning of channel parameters.

III. EM LEARNING OF THE PRIOR PARAMETERS q

The paper uses EM to learn the unknown prior and noise parameters, replacing the intractable exact posterior in the E step with GAMP-approximated posteriors. Incremental parameter updates make the optimization tractable while preserving the EM lower-bound framework.

  • EM framework: EM learns q = [λ, ω, θ, φ, ψ] by iteratively improving a lower bound on the measurement likelihood.The likelihood is guaranteed to converge to a local maximum or at least a saddle point.
  • EM framework: The exact E step would use the posterior pX|Y(x|y; qi), while the M step maximizes the expected complete-data log likelihood.These steps alternate between updating the posterior approximation and the parameter vector.
  • Approximate E step: Because the exact posterior is difficult to calculate, the algorithm constructs the lower bound using products of GAMP-approximated marginal posteriors.The resulting expectations are marked as approximate GAMP expectations.
  • Incremental updates: The algorithm updates q one component at a time rather than solving the joint optimization, using an incremental variant of EM.The notation qi with one element removed denotes the parameters held fixed during each component update.

A. EM Update of the Gaussian Noise Variance ψ

The supplied derivation develops EM updates for Gaussian-mixture signal parameters, including the Bernoulli-Gaussian case. It handles the point mass at zero separately and obtains updates from GAMP-computed posterior quantities.

  • Noise-variance update: The EM update for the Gaussian noise variance ψ is obtained by differentiating the expected likelihood and solving the resulting stationarity condition.The unique solution uses GAMP-derived estimates of the noiseless measurements and their variances.
  • Bernoulli-Gaussian case: The Bernoulli-Gaussian prior is the one-component Gaussian-mixture case, with mixture weight fixed to unity rather than learned.Its prior combines a zero-mass component with one Gaussian component.
  • Sparsity-rate update: The EM update for λ is derived by differentiating the expected log likelihood and treating the neighborhood around zero separately.A Dirac approximation justifies exchanging differentiation and integration in the derivation.
  • Gaussian-parameter updates: The updates for θ and φ are obtained analogously by setting derivatives of the expected log likelihood to zero.The resulting quantities are computed from GAMP outputs and the corresponding posterior expressions.

C. EM Updates of the Signal Parameters: GM Case

For a Gaussian-mixture prior with L>1, EM-GM-AMP derives sequential updates for mixture means, variances, weights, and sparsity, using approximations where exact closed forms are unavailable. Learning means works well for sparse signals but can overfit heavy-tailed signals, motivating a zero-mean mode.

  • For L>1 Gaussian mixtures, exact EM updates are unavailable in closed form, so the derivation uses approximations.The approximations are applied to integrals involving the mixture parameters.
  • The EM procedure incrementally updates each mixture mean, then variance, and finally the full weight vector while holding other parameters fixed.
  • βn,k represents the probability that coefficient xn originated from mixture component k.
  • Learning Gaussian-mixture means yields excellent recovery MSE for sparse signals but can overfit heavy-tailed distributions by gravitating toward outliers.For heavy-tailed signals, fixing the means at zero gives better reconstruction performance.
  • The L=1 approximate Gaussian-mixture updates coincide with the exact EM-Bernoulli-Gaussian updates.

D. EM Initialization

EM-GM-AMP requires careful initialization because EM can converge to a local maximum or saddle point. The procedure uses mode-specific parameter initialization and fixed model-order defaults that generally perform well.

  • EM initialization is essential because the likelihood EM algorithm may converge to a local maximum or saddle point.
  • The heavy-tailed mode initializes means differently from the sparse mode while initializing λ and ψ using the preceding strategy.
  • The default fixed model orders are L = 3 for sparse mode and L = 4 for heavy-tailed mode, and these choices usually perform well.
  • Without user knowledge of the true SNR, the procedure recommends SNR0 = 100 because it works well across a wide range of true SNR values.

E. EM-GM-AMP Summary and Demonstration

EM-GM-AMP estimates the signal distribution by combining a Gaussian-mixture model with AMP-based recovery. In an example using noisy compressed measurements, it learned reasonable approximations in both sparse and heavy-tailed modes.

  • Fixed-L EM-GM-AMP is summarized as the algorithmic procedure, while EM-BG-AMP is obtained by running EM-GM-AMP with L = 1.
  • The distribution-learning demonstration generated signals, compressed noisy measurements, and learned Gaussian-mixture parameters before plotting the resulting approximation.The experiment used N = 2000 coefficients, M = 1000 measurements, and SNR = 25 dB.
  • The learned model provided a reasonable approximation of the unknown true distribution in both sparse and heavy-tailed modes.

F. Selection of GM Model Order L

The method selects the Gaussian-mixture order L by iteratively maximizing a BIC-penalized likelihood lower bound. A demonstration starting from L0 = 1 selected L = 3, with the best MSE at the final order.

  • Model-order selection estimates L by maximizing a BIC-penalized likelihood lower bound rather than the difficult marginal likelihood directly.
  • The penalty is necessary because the unpenalized likelihood lower bound is non-decreasing as L increases.
  • The implementation approximates the posterior with GM-GAMP outputs, uses a point-mass approximation for the remaining log integral, and sets the BIC sample size to U = Σn πn.
  • The procedure starts at L = 1, increments L until the penalized metric decreases, and iterates until the order stabilizes or Jmax is reached.
  • Starting from L0 = 1, the triangular-mixture demonstration selected L = 3 at iteration j = 1, and the final model order had the best MSE.The experiment used N = 1000 coefficients, M = 500 noisy measurements, and SNR = 20 dB.

IV. NUMERICAL RESULTS

The numerical study evaluates EM-GM-AMP under noiseless and noisy recovery settings using fixed algorithmic tolerances, iteration limits, and mixture orders.

  • Experiments assess EM-GM-AMP under both noiseless and noisy settings.
  • The GAMP and EM procedures use tolerances of 10^-5 and maximum iteration counts of 20.
  • Fixed-L EM-GM-AMP uses L = 3 in sparse mode and L = 4 in heavy-tailed mode.

A. Noiseless Phase Transitions

Noiseless experiments compare empirical phase transitions across sparse signal classes, matrix constructions, and algorithms. EM-GM-AMP generally improves phase-transition performance over LASSO-style AMP and LASSO, while matrix structure can alter its behavior.

  • Noiseless Phase Transitions: Empirical PTCs are estimated on a 30 × 30 (M, K) grid with N = 1000 and R = 100 realizations per grid point.Success uses NMSE < 10^-6, and the empirical PTC is the S = 0.5 contour over the sparsity-undersampling grid.
  • Noiseless Phase Transitions: Across Bernoulli-Gaussian, Bernoulli, and Bernoulli-Rademacher signals, EM-GM-AMP significantly improves empirical PTCs over DMM-AMP and theoretical LASSO PTCs.
  • Noiseless Phase Transitions: For Bernoulli-Gaussian signals, EM-GM-AMP, its MOS variant, and EM-BG-AMP have PTCs nearly identical to genie-GM-AMP.This suggests that the EM-learning procedures work well for this signal class.
  • Noiseless Phase Transitions: At M ≪ N with N = 8192, EM-GM-AMP and its MOS variant have nearly identical PTCs and slightly outperform EM-BG-AMP and genie-GM-AMP.The advantage is especially evident at very small measurement counts.
  • Measurement-Matrix Effects: For measurement matrices, the i.i.d-N PTC is preserved by other i.i.d zero-mean sub-Gaussian matrices and randomly row-sampled DCT matrices.
  • Measurement-Matrix Effects: The EM-GM-AMP PTC can degrade with non-zero-mean i.i.d matrices or super-Gaussian matrices, although i.i.d-BR matrices with λA = 0.15 preserve the i.i.d-N PTC.The BR sub-Gaussian condition requires λA ≥ 1/3, according to the cited background passage.

B. Noisy Sparse Signal Recovery

Noisy experiments evaluate NMSE across sparse signal classes and SNRs. EM-GM-AMP variants achieve the strongest reported recovery performance, with gains varying by signal distribution and a dependence on noise-variance knowledge for T-MSBL.

  • Noisy Sparse Signal Recovery: Noisy NMSE experiments fix N = 1000, K = 100, and SNR = 25 dB while varying M and averaging R = 500 realizations.
  • Noisy Sparse Signal Recovery: For Bernoulli-Gaussian signals, EM-GM-AMP variants reduce the M/N breakpoint from 0.3 to 0.26 and improve NMSE by ≈1 dB over BCS.
  • Noisy Sparse Signal Recovery: For Bernoulli signals, EM-GM-AMP variants reduce the M/N breakpoint from 0.4 to 0.32, or 0.3 with MOS, and reduce NMSE by ≈8 dB over T-MSBL.
  • Noisy Sparse Signal Recovery: For Bernoulli-Rademacher signals, EM-GM-AMP-MOS reduces NMSE by 10 dB for M/N ≥ 0.36 and lowers the breakpoint from 0.38 to 0.35.
  • SNR Robustness: Across the full SNR range, EM-GM-AMP and EM-GM-AMP-MOS yield significantly lower NMSE than the other evaluated algorithms for Bernoulli-Rademacher recovery.
  • SNR Robustness: T-MSBL must receive some knowledge of the true noise variance to perform well, unlike the proposed algorithms.

C. Heavy-Tailed Signal Recovery

EM-GM-AMP performs strongly on heavy-tailed Student’s-t and log-normal signals, including cases where the signal is non-compressible. Its performance is attributed to adapting the learned model to the signal and realization.

  • Student’s-t signals: EM-GM-AMP outperformed all tested algorithms for Student’s-t signals with rate q = 1.67.The experiment used N = 1000, SNR = 25 dB, and R = 500 realizations.
  • Student’s-t signals: EM-GM-AMP variants also showed similarly good performance for other Student’s-t rates and centered Cauchy signals.
  • Log-normal signals: EM-GM-AMP-MOS, EM-GM-AMP, and EM-BG-AMP showed excellent performance for positive log-normal signals with μ = 0 and σ2 = 1.
  • Experimental setup: The heavy-tailed experiments evaluate NMSE as a function of undersampling ratio M/N for Student’s-t and log-normal signals.
  • Overall result: Across the signal types tested, EM-GM-AMP and its MOS variant achieved the best recovery performance.The authors attribute this to tuning itself to the signal and realization at hand.

D. Runtime and Complexity Scaling with N

EM-GM-AMP has favorable high-dimensional complexity, scaling as O(MN) with dense sensing matrices and O(N log N) when fast FFT-based multiplications are available. Runtime advantages become more pronounced for long signals, while audio experiments show strong reconstruction performance and improving scalability.

  • Complexity scaling: O(MN) is the expected dense-matrix complexity of EM-GM-AMP, becoming O(N^2) when M/N is fixed.The runtime curves exhibit the corresponding slope-2 scaling at sufficiently large N.
  • Complexity scaling: O(N^3) scaling makes T-MSBL, BCS, SL0, OMP, and SP orders of magnitude slower than EM-GM-AMP for long signals such as N ≥ 10^4.For short signals, including N < 1300, several competing algorithms remain faster.
  • NMSE scaling: With fixed undersampling and sparsity ratios, most algorithms’ NMSE is relatively insensitive to signal length N, while EM-GM-AMP improves with N.
  • Fast sensing operators: O(N log N) complexity is available for EM-GM-AMP and SPGL1 when multiplication by A and A^T uses a fast FFT-based implementation.This also avoids storing a dense M×N matrix in memory.
  • Audio recovery: The audio experiment reconstructed DCT transform coefficients blockwise from compressed measurements, using either i.i.d. Gaussian or randomly selected sensing operators.Additional multitime or multidimensional coefficient structure was outside the paper’s scope.
  • Audio recovery: In compressive audio recovery, EM-GM-AMP methods and SL0 achieved the best TNMSE, while EM-GM-AMP runtimes improved as block length increased.At N = 8192 with i.i.d. Gaussian sensing, EM-BG-AMP and fixed-L EM-GM-AMP became the two fastest tested algorithms.
Loading 1207.3107v3…