{'fraction_non_alphanumeric': 0.06338387445666592, 'fraction_numerical': 0.018468257045487012, 'mean_word_length': 3.9374307304785896, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Expectation-Maximization Gaussian-Mixture Approximate Message Passing

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
We consider estimating a K-sparse (or compressible) signal x ∈ R N from M < N linear measurements y = Ax + w ∈ R M , where A is known and w is additive white Gaussian noise (AWGN). For this problem, accurate (relative to the noise variance) signal recovery is known to be possible with polynomial-complexity algorithms when x is sufficiently sparse and when A satisfies certain restricted isometry properties [4], or when A is large with i.i.d zeromean sub-Gaussian entries [5] as discussed below.
LASSO [6] (or, equivalently, Basis Pursuit Denoising [7]), is a well-known approach to the sparse-signal recovery problem that solves the convex problem with λ lasso a tuning parameter that trades between the sparsity and measurement-fidelity of the solution. When A is constructed from i.i.d zero-mean sub-Gaussian entries, the performance of LASSO can be sharply characterized in the large system limit (i.e., as K, M, N → ∞ with fixed undersampling ratio M/N and sparsity ratio K/M ) using the so-called phase transition curve (PTC) [5], [8]. When the observations are noiseless, the PTC bisects the M/Nversus-K/M plane into the region where LASSO reconstructs the signal perfectly (with high probability) and the region where it does not. (See Figs. 3-5.) When the observations are noisy, the same PTC bisects the plane into the regions where LASSO's noise sensitivity (i.e., the ratio of estimationerror power to measurement-noise power under the worst-case signal distribution) is either finite or infinite [9]. An important fact about LASSO's noiseless PTC is that it is invariant to the distribution of the nonzero signal coefficients. In other words, if the vector x is drawn i.i.d from the pdf where δ(·) is the Dirac delta, f X (·) is the active-coefficient pdf (with zero probability mass at x = 0), and λ K/N , then the LASSO PTC is invariant to f X (·). While this implies that LASSO is robust to "difficult" instances of f X (·), it also implies that LASSO cannot benefit from the case that f X (·) is an "easy" distribution. For example, when the signal is known apriori to be nonnegative, polynomial-complexity algorithms exist with PTCs that are better than LASSO's [10].
At the other end of the spectrum is minimum mean-squared error (MMSE)-optimal signal recovery under known marginal pdfs of the form (2) and known noise variance. The PTC of MMSE recovery has been recently characterized [11] and shown to be well above that of LASSO. In particular, for any f X (·), the PTC on the M/N -versus-K/M plane reduces to the line K/M = 1 in both the noiseless and noisy cases. Moreover, efficient algorithms for approximate MMSE-recovery have been proposed, such as the Bayesian version of Donoho, Maleki, and Montanari's approximate message passing (AMP) algorithm from [12], which performs loopy belief-propagation on the underlying factor graph using central-limit-theorem approximations that become exact in the large-system limit under i.i.d zero-mean sub-Gaussian A. In fact, in this regime, AMP obeys [13] a state-evolution whose fixed points, when unique, are optimal. To handle arbitrary noise distributions and a wider class of matrices A, Rangan proposed a generalized AMP (GAMP) [14] that forms the starting point of this work. (See Table I.) For more details and background on GAMP, we refer the reader to [14].
In practice, one ideally wants a recovery algorithm that does not need to know p X (·) and the noise variance a priori, yet offers performance on par with MMSE recovery, which (by definition) requires knowing these prior statistics. Towards this goal, we propose a recovery scheme that aims to learn the prior signal distribution p X (·), as well as the variance of the AWGN, while simultaneously recovering the signal vector x from the noisy compressed measurements y. To do so, we model the active component f X (·) in (2) using a generic L-term Gaussian mixture (GM) and then learn the GM parameters and noise variance using the expectation-maximization (EM) algorithm [15]. As we will see, all of the quantities needed for the EM updates are already computed by the GAMP algorithm, making the overall process very computationally efficient. Moreover, GAMP provides approximately MMSE estimates of x that suffice for signal recovery, as well as posterior activity probabilities that suffice for support recovery.
Since, in our approach, the prior pdf parameters are treated as deterministic unknowns, our proposed EM-GM-AMP algorithm can be classified as an "empirical-Bayesian" approach [16]. Compared with previously proposed empirical-Bayesian approaches to compressive sensing (e.g., [17]- [19]), ours has a more flexible signal model, and thus is able to better match a wide range of signal pdfs p X (·), as we demonstrate through a detailed numerical study. In addition, the complexity scaling of our algorithm is superior to that in [17]- [19], implying lower complexity in the high dimensional regime, as we confirm numerically. Supplemental experiments demonstrate that our excellent results hold for a wide range of sensing operators A, with some exceptions. Although this paper does not contain any convergence guarantees or a rigorous analysis/justification of the proposed EM-GM-AMP, Kamilov et al. showed (after the submission of this work) in [20] that a generalization of EM-GM-AMP yields asymptotically (i.e., in the large system limit) consistent parameter estimates when A is i.i.d zeromean Gaussian, when the parameterized signal and noise distributions match the true signal and noise distributions, and when those distributions satisfy certain identifiability conditions. We refer interested readers to [20] for more details.
Notation: For matrices, we use boldface capital letters like A, and we use tr(A) and A F to denote the trace and Frobenius norm, respectively. Moreover, we use (·) T , (·) * , and (·) H to denote transpose, conjugate, and conjugate transpose, respectively. For vectors, we use boldface small letters like x, and we use x p = ( n |x n | p ) 1/p to denote the ℓ p norm, with x n representing the n th element of x. For a Gaussian random vector x with mean m and covariance matrix Q, we denote the pdf by N (x; m, Q), and for its circular complex Gaussian counterpart, we use CN (x; m, Q). Finally, we use E{·}, δ(·), R, and C to denote the expectation operation, the Dirac delta, the real field, and the complex field, respectively.

II. GAUSSIAN-MIXTURE GAMP
We first introduce Gaussian-mixture (GM) GAMP, a key component of our overall approach, where the coefficients in x = [x 1 , . . . , x N ] T are assumed to be i.i.d with marginal pdf where δ(·) is the Dirac delta, λ is the sparsity rate, and, for the k th GM component, ω k , θ k , and φ k are the weight, mean, and variance, respectively. In the sequel, we use ω [ω 1 , . . . , ω L ] T and similar definitions for θ and φ. By definition, L ℓ=1 ω ℓ = 1. The noise w = [w 1 , . . . , w M ] T is assumed to be i.i.d Gaussian, with mean zero and variance ψ, i.e., and independent of x. Although above and in the sequel we assume real-valued quantities, all expressions in the sequel can be converted to the circular-complex case by replacing N with CN and removing the 1 2 's from (25), (44), and (58). We note that, from the perspective of GM-GAMP, the prior parameters q [λ, ω, θ, φ, ψ] and the number of mixture components, L, are treated as fixed and known.
GAMP models the relationship between the m th observed output y m and the corresponding noiseless output z m a T m x, where a T m denotes the m th row of A, using the conditional pdf p Y |Z (y m |z m ; q). It then approximates the true marginal posterior p(z m |y; q) by using quantitiesp m and µ p m that change with iteration t (see Table I), although here we suppress the t notation for brevity. Under the AWGN assumption 1 (4) we have p Y |Z (y|z; q) = N (y; z, ψ), and thus the pdf (5) has moments [14] GAMP then approximates the true marginal posterior p(x n |y; q) by p X|Y (x n |y;r n , µ r n , q) p X (x n ; q) N (x n ;r n , µ r n ) x p X (x; q) N (x;r n , µ r n ) where againr n and µ r n vary with the GAMP iteration t. Plugging the sparse GM prior (3) into (8) and simplifying, one can obtain 2 the GM-GAMP approximated posterior p X|Y (x n |y;r n , µ r n , q) = 1 − π n δ(x n ) + π n L ℓ=1 β n,ℓ N x n ; γ n,ℓ , ν n,ℓ 1 Because GAMP can handle an arbitrary p Y |Z (·|·), the extension of EM-GM-AMP to additive non-Gaussian noise, and even non-additive measurement channels (such as with quantized outputs [21] or logistic regression [14]), is straightforward. Moreover, the parameters of the pdf p Y |Z (·|·) could be learned using a method similar to that which we propose for learning the AWGN variance ψ, as will be evident from the derivation in Section III-A. Finally, one could even model p Y |Z (·|·) as a Gaussian mixture and learn the corresponding parameters. 2 Both (10) and (12) can be derived from (9) via the Gaussian-pdf multiplication rule: with normalization factor ζ n x p X (x; q) N (x;r n , µ r n ) = (1−λ)N (0;r n , µ r n )+λ L ℓ=1 ω ℓ N (0;r n −θ ℓ , µ r n +φ ℓ ) (12) and (r n , µ r n , q)-dependent quantities β n,ℓ λω ℓ N (r n ; θ ℓ , φ ℓ + µ r n ) The posterior mean and variance of p X|Y are given in steps (R9)-(R10) of Table I, and (10) makes it clear that π n is GM-GAMP's approximation of the posterior support probability Pr{x n = 0 | y; q}.
In principle, one could specify GAMP for an arbitrary signal prior p X (·). However, if the integrals in (R9)-(R10) are not computable in closed form (e.g., when p X (·) is Student's-t), then they would need to be computed numerically, thereby drastically increasing the computational complexity of GAMP. In contrast, for GM signal models, we see above that all steps can be computed in closed form. Thus, a practical approach to the use of GAMP with an intractable signal prior p X (·) is to approximate p X (·) using an L-term GM, after which all GAMP steps can be easily implemented. The same approach could also be used to ease the implementation of intractable output priors p Y |Z (·|·).

III. EM LEARNING OF THE PRIOR PARAMETERS q
We now propose an expectation-maximization (EM) algorithm [15] to learn the prior parameters q [λ, ω, θ, φ, ψ]. The EM algorithm is an iterative technique that increases a lower bound on the likelihood p(y; q) at each iteration, thus guaranteeing that the likelihood converges to a local maximum or at least a saddle point [22]. In our case, the EM algorithm manifests as follows. Writing, for arbitrary pdfp(x), Lp(y; q) where  iterates over two steps: E) choosingp to maximize the lower bound for fixed q = q i , and M) choosing q to maximize the lower bound for fixedp =p i . For the E step, since Lp(y; q i ) = ln p(y; q i ) − D(p p X|Y (·|y; q i )), the maximizing pdf would clearly bep i (x) = p X|Y (x|y; q i ), i.e., the true posterior under prior parameters q i . Then, for the M step, since Lpi (y; q) = Epi (x) {ln p(x, y; q)} + H(p i ), the maximizing q would clearly be q i+1 = arg max q E{ln p(x, y; q) | y; q i }.
In our case, because the true posterior is very difficult to calculate, we instead construct our lower-bound Lp(y; q) using the GAMP approximated posteriors, i.e., we setp i (x) = n p X|Y (x n |y; q i ) for p X|Y defined in (8), resulting in where "Ê" indicates the use of the GAMP's posterior approximation. Moreover, since the joint optimization in (21) is difficult to perform, we update q one component at a time (while holding the others fixed), which is the well known "incremental" variant on EM from [23]. In the sequel, we use "q i \λ " to denote the vector q i with the element λ removed (and similar for the other parameters).

A. EM Update of the Gaussian Noise Variance ψ
We first derive the EM update for the noise variance ψ given a previous parameter estimate q i . For this, we write p(x, y; q) = Cp(y|x; ψ) = C M m=1 p Y |Z (y m |a T m x; ψ) for a ψ-invariant constant C, so that The maximizing value of ψ in (23) is necessarily a value of ψ that zeroes the derivative of the sum, i.e., that satisfies 3 M m=1 zm Because p Y |Z (y m |z m ; ψ) = N (y m ; z m , ψ), we can obtain which, when plugged into (24), yields the unique solution where the use ofẑ m and µ z m follows from (R3)-(R4) in Table I.

B. EM Updates of the Signal Parameters: BG Case
Suppose that the signal distribution p X (·) is modeled using an L = 1-term GM, i.e., a Bernoulli-Gaussian (BG) pdf. In this case, the marginal signal prior in (3) reduces to Note that, in the BG case, the mixture weight ω is, by definition, unity and does not need to be learned. We now derive the EM update for λ given previous param- Because we can write p(x, y; q) = C N n=1 p X (x n ; λ, θ, φ) for a λ-invariant constant C, The maximizing value of λ in (29) is necessarily a value of λ that zeroes the derivative of the sum, i.e., that satisfies 4 N n=1 xn For the BG p X (x n ; λ, θ, φ) in (28), it is readily seen that Plugging (32) and (9) into (30), it becomes evident that the neighborhood around the point x n = 0 should be treated differently than the remainder of R. Thus, we define the closed ball B ǫ [−ǫ, ǫ] and its complement B ǫ R \ B ǫ , and note that, in the limit ǫ → 0, the following is equivalent to (30): where the values taken by the integrals are evident from (10). Finally, the EM update for λ is the unique value satisfying (33) as ǫ → 0, which is readily shown to be Conveniently, the posterior support probabilities {π n } N n=1 are easily calculated from the GM-GAMP outputs via (15).
Similar to (29), the EM update for θ can be written as The maximizing value of θ in (35) is again a necessarily a value of θ that zeroes the derivative, i.e., that satisfies N n=1 xn For the BG p X (x n ; λ, θ, φ) given in (28), Splitting the domain of integration in (36) into B ǫ and B ǫ as before, and then plugging in (38), we find that the following is equivalent to (36) in the limit of ǫ → 0: The unique value of θ satisfying (39) as ǫ → 0 is then where {γ n,1 } N n=1 defined in (16) are easily computed from the GM-GAMP outputs. The equality in (41) can be verified by plugging the GAMP posterior expression (10) into (40).
Similar to (29), the EM update for φ can be written aŝ The maximizing value of φ in (42) is again necessarily a value of φ that zeroes the derivative, i.e., that satisfies N n=1 xn For the p X (x n ; λ, θ, φ) given in (28), it is readily seen that Splitting the domain of integration in (43) into B ǫ and B ǫ as before, and then plugging in (44), we find that the following is equivalent to (43) in the limit of ǫ → 0: The unique value of φ satisfying (45) as ǫ → 0 is then where {ν n,1 } N n=1 from (17) are easily computed from the GAMP outputs. The equality in (47) can be readily verified by plugging (10) into (46).

C. EM Updates of the Signal Parameters: GM Case
We now generalize the EM updates derived in Section III-B to the GM prior given in (3) for L ≥ 1. As we shall see, it is not possible to write the exact EM updates in closed-form when L > 1, and so some approximations will be made.
We begin by deriving the EM update for λ given the previous parameters The first two steps are identical to the steps (29) and (30) presented for the BG case, and for brevity we do not repeat them here. In the third step, use of the GM prior (3) yields which coincides with the BG expression (32). The remaining steps also coincide with those in the BG case, and so the final EM update for λ, in the case of a GM, 5 is given by (34). We next derive the EM updates for the GM parameters ω, θ, and φ. For each k = 1, . . . , L, we incrementally update θ k , then φ k , and then the entire vector ω, while holding all other parameters fixed. The EM updates are thus Following (36), the maximizing value of θ k in (49) is again necessarily a value of θ k that zeros the derivative, i.e., N n=1 xn Plugging in the derivative and the version of p X|Y (x n |y; q i ) from (9), integrating (52) separately over B ǫ and B ǫ as in (33), and taking ǫ → 0, we find that the B ǫ portion vanishes, giving the necessary condition N n=1 xn (54) Since this integral cannot be evaluated in closed form, we apply the approximation N (x n ; θ k , φ i k ) ≈ N (x n ; θ i k , φ i k ) in both the numerator and denominator, and subsequently exploit the fact that p(x n |x n = 0, y; q i ) = N (x n ;r n , µ r n ) ℓ ω i ℓ N (x n ; θ i ℓ , φ i ℓ ) from (9) to cancel terms, and so obtain the (approximated) necessary condition N n=1 xn (55) We then simplify (55) using the Gaussian-pdf multiplication rule, and set θ i+1 k equal to the value of θ k that satisfies (55), which can be found to be Note from (10) that π n β n,k can be interpreted as the probability that x n originated from the k th mixture component.
For sparse signals x, we find that learning the GM means {θ k } using the above EM procedure yields excellent recovery MSE. However, for "heavy-tailed" signals (i.e., whose pdfs have tails that are not exponentially bounded, such as Student's-t), our experience indicates that the EM-learned values of {θ k } tend to gravitate towards the outliers in {x n } N n=1 , resulting in an overfitting of p X (·) and thus poor reconstruction MSE. For such heavy-tailed signals, we find that better reconstruction performance is obtained by fixing the means at zero (i.e., θ i k = 0 ∀k, i). Thus, in the remainder of the paper, we consider two modes of operation: a "sparse" mode where θ is learned via the above EM procedure, and a "heavy-tailed" mode that fixes θ = 0. Following (52), the maximizing value of φ k in (50) is necessarily a value of φ k that zeroes the derivative, i.e., N n=1 xn As for the derivative in the previous expression, we find .
Integrating (57) separately over B ǫ and B ǫ , as in (33), and taking ǫ → 0, we find that the B ǫ portion vanishes, giving N n=1 xn Similar to (54), this integral is difficult to evaluate, and so we again apply the approximation N ( in the numerator and denominator, after which several terms cancel, yielding the necessary condition To find the value of φ k satisfying (60), we expand |x n − θ i k | 2 = |x n | 2 − 2 Re(x * n θ i k ) + |θ i k | 2 and apply the Gaussian-pdf multiplication rule, which gives Finally, the value of the positive ω maximizing (51) under the pmf constraint L k=1 ω k = 1 can be found by solving the unconstrained optimization problem max ω,ξ J(ω, ξ), where ξ is a Lagrange multiplier and We start by setting d dω k J(ω, ξ) = 0, which yields N n=1 xn Like in (54) and (59), the above integral is difficult to evaluate, and so we approximate ω ≈ ω i , which reduces the previous equation to Multiplying both sides by ω i k for k = 1, . . . , L, summing over k, employing the fact 1 = k ω i k , and simplifying, we obtain the equivalent condition Plugging (67) into (65) and multiplying both sides by ω k , the derivative-zeroing value of ω k is seen to be where, if we use ω k ≈ ω i k on the right of (68), then we obtain Although, for the case of GM priors, approximations were used in the derivation of the EM updates (56), (61), and (69), it is interesting to note that, in the case of L = 1 mixture components, these approximate EM-GM updates coincide with the exact EM-BG updates derived in Section III-B. In particular, the approximate-EM update of the GM parameter θ 1 in (56) coincides with the exact-EM update of the BG parameter θ in (41), the approximate-EM update of the GM parameter φ 1 in (61) coincides with the exact-EM update of the BG parameter φ in (47), and the approximate-EM update of the GM parameter ω 1 in (69) reduces to the fixed value 1. Thus, one can safely use the GM updates above in the BG setting without any loss of optimality.

D. EM Initialization
Since the EM algorithm may converge to a local maximum or at least a saddle point of the likelihood function, proper initialization of the unknown parameters q is essential. Here, we propose initialization strategies for both the "sparse" and "heavy-tailed" modes of operation, for a given value of L. Regarding the value of L, we prescribe a method to learn it in Section III-F. However, the fixed choices L = 3 for "sparse" mode and L = 4 for "heavy tailed" mode usually perform well, as shown in Section IV.
For the "sparse" mode, we set the initial sparsity rate λ 0 equal to the theoretical noiseless LASSO PTC, i.e., describes the maximum value of K M supported by LASSO for a given M N , and where Φ(·) and φ(·) denote the cdf and pdf of 7 the N (0, 1) distribution, respectively. Using the energies ||y|| 2 2 and ||A|| 2 F and an assumed value of SNR 0 , we initialize the noise and signal variances, respectively, as where, in the absence of (user provided) knowledge about the true SNR Ax 2 2 / w 2 2 , we suggest SNR 0 = 100, because in our experience this value works well over a wide range of true SNR. Then, we uniformly space the initial GM means θ 0 over [ −L+1 2L , L−1 2L ], and subsequently fit the mixture weights ω 0 and variances φ 0 to the uniform pdf supported on [−0.5, 0.5] (which can be done offline using the standard approach to EM-fitting of GM parameters, e.g., [24, p. 435]). Finally, we multiply θ 0 by 12ϕ 0 and φ 0 by 12ϕ 0 to ensure that the resulting signal variance equals ϕ 0 .

E. EM-GM-AMP Summary and Demonstration
The fixed-L EM-GM-AMP 6 algorithm developed in the previous sections is summarized in Table II. For EM-BG-AMP (as previously described in [2]), one would simply run EM-GM-AMP with L = 1.
To demonstrate EM-GM-AMP's ability to learn the underlying signal distribution, Fig. 1 shows examples of the GM-modeled signal distributions learned by EM-GM-AMP in both "sparse" and "heavy-tailed" modes. To create the figure, we first constructed the true signal vector x ∈ R N using N = 2000 independent draws of the true distribution p X (·) shown in each of the subplots. Then, we constructed measurements y = Ax + w by drawing A ∈ R M×N with i.i.d N (0, M −1 ) elements and w ∈ R M with i.i.d N (0, σ 2 ) elements, with M = 1000 and σ 2 chosen to achieve SNR = 25 dB. Finally, we ran EM-GM-AMP according to Table II, and plotted the GM approximation p X (x; q i ) from (3) using the learned pdf parameters Figure 1 confirms that EM-GM-AMP is successful in learning a reasonable approximation of the unknown true pdf p X (·) from the noisy compressed observations y, in both sparse and heavy-tailed modes.

F. Selection of GM Model Order L
We now propose a method to learn the number of GM components, L, based on standard maximum likelihood (ML)based model-order-selection methodology [25], i.e., arg max whereq L is the ML estimate of q under the hypothesis L and η(L) is a penalty term. For η(L), there are several possibilities, but we focus on the Bayesian information criterion (BIC) [25]: 6 Matlab code at http://www.ece.osu.edu/ ∼ schniter/EMturboGAMP.
The top subplot shows "sparse" mode EM-GM-AMP run using GM-order L = 3 on a sparse signal whose nonzero components were generated according to a triangular mixture, whereas the bottom subplot shows "heavy-tailed" EM-GM-AMP run using L = 4 on a Student's-t signal with rate parameter q = 1.67 (defined in (82)). The density of the continuous component λf X (x) is marked on the left axis, while the mass of the discrete component (1 − λ)δ(x) is marked on the right axis.
Initialize L and q 0 as described in Section III-D. Table I).
Compute λ i from π i−1 as described in (34). for k = 1 to L do if sparse mode enabled then Compute θ i k from π i−1 , γ i−1 k , {β i−1 l } L l=1 as described in (56). else if heavy-tailed mode enabled then Compute ω i from π i−1 and {β i−1 l } L l=1 as described in (69). end for Compute ψ i fromẑ i and (µ z ) i as in (27). end for

TABLE II THE EM-GM-AMP ALGORITHM (FIXED-L CASE)
where |q L | denotes the number 7 of real-valued parameters affected by L, and U is the sample size (see below).
Because ln p(y;q L ) is difficult to evaluate, we work with 7 In our case, the parameters affected by L are the GM means, variances, and weights, so that, for real-valued signals, we use |q L | = 3L − 1 in "sparse" mode and |q L | = 2L − 1 in heavy-tailed mode, and for complex-valued signals, we use |q L | = 4L − 1 in "sparse" mode and |q L | = 2L − 1 in heavy-tailed mode. the lower bound (where for now L j ,q L , andq L j are arbitrary) ln p(y;q L ) = ln x p(x|y;q L j ) p(x, y;q L ) p(x|y;q L j ) (75) ≥ x p(x|y;q L j ) ln p(x, y;q L ) p(x|y;q L j ) (76) = x p(x|y;q L j ) ln p(x, y;q L ) + const (77) = N n=1 xn p(x n |y;q L j ) ln p X (x n ;q L ) + const (78) p(x n |y;q L j ) ln f X (x n ;q L ) where (76) applies Jensen's inequality, "const" denotes a constant term w.r.t L, and (78) holds because ln p(x, y;q L ) = ln p(x;q L ) + ln p(y|x;ψ) = N n=1 ln p X (x n ;q L ) + const. Equation (79) can then be obtained integrating (78) separately over B ǫ and B ǫ and taking ǫ → 0, as done several times in Section III-B. Using this lower bound in place of ln p(y;q L ) in (73), we obtain the BIC-inspired model order estimate (where nowq L is specifically the ML estimate of q L ) We in fact propose to perform (80) iteratively, with j = 0, 1, 2, . . . denoting the iteration index. Notice that (80) can be interpreted as a "penalized" EM update for L; if we neglect the penalty term η(L), then (75)-(79) becomes a standard derivation for the EM-update of L (recall, e.g., the EM derivation in Section III). The penalty term is essential, though, because the unpenalized log-likelihood lower bound L L j (y;q L ) is non-decreasing 8 in L.
We now discuss several practical aspects of our procedure. First, we are forced to approximate the integral in (79). To start, we use GM-GAMP's approximation of the posterior p(x n |y;q L j ) from (9), and the EM approximations of the MLestimatesq L j andq L outlined in Section III-C. In this case, the integral in (79) takes the form xn π n L j l=1 β n,l N (x n ; γ n,l , ν n,l ) ln which is still difficult due to the log term. Hence, we evaluate (81) using the point-mass approximation N (x n ; γ n,l , ν n,l ) ≈ δ(x n − γ n,l ). Second, for the BIC penalty (74), we use the sample size U = N n=1 π n , which is the effective number of terms in the sum in (79). Third, when maximizing L over Z + in (80), we start with L = 1 and increment L in steps of one until the penalized metric decreases. Fourth, for the initial model order L 0 , we recommend using L 0 = 3 in "sparse" mode and L 0 = 4 in "heavy-tailed" mode, i.e., the fixed-L defaults from Section III-D. Finally, (80) is iterated until either L j+1 = L j or a predetermined maximum number of allowed model-order iterations J max has been reached.
As a demonstration of the proposed model-order selection procedure, we estimated a realization of x with N = 1000 coefficients drawn i.i.d from the triangular mixture pdf shown in Fig. 1 (top, red) with λ = 0.1, from the M = 500 noisy measurements y = Ax + w, where A was i.i.d N (0, M −1 ), and w was AWGN such that SNR = 20 dB. For illustrative purposes, we set the initial model order at L 0 = 1. Iteration j = 1 yielded the metric L L j (y;q L ) − η BIC (L) shown at the top of Fig. 2, which was maximized by L = 3 L 1 . The metric resulting from iteration j = 2 is shown in the middle of Fig. 2, which was maximized by L = 2 L 2 . At iteration j = 3, we obtained the metric at the bottom of Fig. 2, which is also maximized by L = 2 L 3 . Since L 3 = L 2 , the algorithm terminates with final model order estimate L = 2. Figure 2 also indicates the per-iteration MSE, which is best at the final model order.

IV. NUMERICAL RESULTS
In this section we report the results of a detailed numerical study that investigate the performance of EM-GM-AMP under both noiseless and noisy settings. For all experiments, we set the GM-GAMP tolerance to τ gamp = 10 −5 and the maximum GAMP-iterations to T max = 20 (recall Table I), and we set the EM tolerance to τ em = 10 −5 and the maximum EM-iterations to I max = 20 (recall Table II). For fixed-L EM-GM-AMP, we set L = 3 in "sparse" and L = 4 in "heavy-tailed" modes. From the noiseless measurements y = Ax, we recovered the signal x using several algorithms. A recoveryx from realization r ∈ {1, . . . , R} was defined a success if the NMSE x −x 2 2 / x 2 2 < 10 −6 , and the average success rate was defined as S 1 R R r=1 S r , where S r = 1 for a success and S r = 0 otherwise. The empirical PTC was then plotted, using Matlab's contour command, as the S = 0.5 contour over the sparsity-undersampling grid. Figures 3-5 show the empirical PTCs for five recovery algorithms: the proposed EM-GM-AMP algorithm (in "sparse" mode) for both L fixed and L learned through model-order selection (MOS), the proposed EM-BG-AMP algorithm, a genie-tuned 9 GM-AMP that uses the true parameters q = [λ, ω, θ, φ, ψ], and the Donoho/Maleki/Montanari (DMM) LASSO-style AMP from [10]. For comparison, Figs. 3-5 also display the theoretical LASSO PTC (70). The signals were generated as Bernoulli-Gaussian (BG) in Fig. 3 (using mean θ = 0 and variance φ = 1 for the Gaussian component), as Bernoulli in Fig. 4 (i.e., all non-zero coefficients set equal to 1), and as Bernoulli-Rademacher (BR) in Fig. 5.

A. Noiseless Phase Transitions
For all three signal types, Figs. 3-5 show that the empirical PTC of EM-GM-AMP significantly improves on the empirical PTC of DMM-AMP as well as the theoretical PTC of LASSO. (The latter two are known to converge in the large system limit [10].) For BG signals, Fig. 3 shows that EM-GM-AMP-MOS, EM-GM-AMP, and EM-BG-AMP all yield PTCs that are nearly identical to that of genie-GM-AMP, suggesting that our EM-learning procedures are working well. For Bernoulli signals, Fig. 4 shows EM-GM-AMP-MOS performing very close to genie-GM-AMP, and both EM-GM-AMP and EM-BG-AMP performing slightly worse but far better than DMM-AMP. Finally, for BR signals, Fig. 5 shows EM-GM-AMP performing significantly better than EM-BG-AMP, since the former is able to accurately model the BR distribution (with L ≥ 2 mixture components) whereas the latter (with a single mixture component) is not, and on par with genie-GM-AMP, whereas EM-GM-AMP-MOS performs noticeably better than genie-GM-AMP. The latter is due to EM-GM-AMP-MOS doing per-realization parameter tuning, while genie-GM-AMP employs the best set of fixed parameters over all realizations.
To better understand the performance of EM-GM-AMP when M N ≪ 1, we fixed N = 8192 and constructed a 12 × 9 grid of (M, K) values spaced uniformly in the log domain. At each grid point, we generated R = 100 independent realizations of a K-sparse BG signal and an i.i.d N (0, M −1 ) matrix A. We then recovered x from the noiseless measurements using EM-GM-AMP-MOS, EM-GM-AMP, EM-BG-AMP, genie-GM-AMP, and the Lasso-solver 10 FISTA 11 [26]. Figure 6 shows that the PTCs of EM-GM-AMP-MOS and EM-GM-AMP are nearly identical, slightly better than those of EM-BG-AMP and genie-GM-AMP (especially at very small 9 For genie-tuned GM-AMP, for numerical reasons, we set the noise variance at ψ = 10 −6 and, with Bernoulli and BR signals, the mixture variances at φ k = 10 −2 . 10 For this experiment, we also tried DMM-AMP but found that it had convergence problems, and we tried SPGL1 but found performance degradations at small M . 11 For FISTA, we used the regularization parameter λ FISTA = 10 −5 , which is consistent with the values used for the noiseless experiments in [26]. M ), and much better than FISTA's. Next, we studied the effect of the measurement matrix construction on the performance of EM-GM-AMP in "sparse" mode with fixed L = 3. For this, we plotted EM-GM-AMP empirical PTCs for noiseless recovery of a length-N = 1000 BG signal under several types of measurement matrix A: i.i. d N (0, 1) . This is not surprising given that AMP itself has rigorous guarantees for i.i.d zeromean sub-Gaussian matrices [5]. Figure 7 shows that the i.i.d-N PTC is also preserved with randomly row-sampled DCT matrices, which is not surprising given AMP's excellent empirical performance with many types of deterministic A [27] even in the absence of theoretical guarantees. Figure 7 shows, however, that EM-GM-AMP's PTC can degrade with non-zero-mean i.i.d matrices (as in the Bernoulli example) or with super-Gaussian i.i.d matrices (as in the BR example with sparsity rate λ A = 0.05 and the Cauchy example). Surprisingly, the i.i.d-N PTC is preserved by i.i.d-BR matrices with sparsity rate λ A = 0.15, even though λ A > 1 3 is required for a BR matrix to be sub-Gaussian [28]. Figures 8-10 show NMSE for noisy recovery of BG, Bernoulli, and BR signals, respectively. To construct these plots, we fixed N = 1000, K = 100, SNR = 25 dB, and varied M . Each data point represents NMSE averaged over R = 500 realizations, where in each realization we drew an A with i.i.d N (0, M −1 ) elements, an AWGN noise vector, and a random signal vector. For comparison, we show the performance of the proposed EM-GM-AMP (in "sparse" mode) for both MOS and L = 3 versions, EM-BG-AMP, genietuned 13 Orthogonal Matching Pursuit (OMP) [29], genietuned 13 Subspace Pursuit (SP) [30], Bayesian Compressive Sensing (BCS) [19], Sparse Bayesian Learning [18] (via the 12 For the Bernoulli and BR matrices, we ensured that no two columns of a given realization A were identical. 13 We ran both OMP (using the implementation from http://sparselab.stanford.edu/OptimalTuning/code.htm) and SP under 10 different sparsity assumptions, spaced uniformly from 1 to 2K, and reported the lowest NMSE among the results. For BG signals, Fig. 8 shows that EM-GM-AMP-MOS, EM-GM-AMP, and EM-BG-AMP together exhibit the best performance among the tested algorithms, reducing the M/N breakpoint (i.e., the location of the knee in the NMSE curve, which represents a sort of phase transition) from 0.3 down to 0.26, but also improving NMSE by ≈ 1 dB relative to the next best algorithm, which was BCS. Relative to the other EM-AMP variants, MOS resulted in a slight degradation of performance for M N between 0.26 and 0.31, but was otherwise identical. For Bernoulli signals, Fig. 9 shows much more significant gains for EM-GM-AMP-MOS, EM-GM-AMP and EM-BG-AMP over the other algorithms: the M/N breakpoint was reduced from 0.4 down to 0.32 (and even 0.3 with MOS), and the NMSE was reduced by ≈ 8 dB relative to the next best algorithm, which was T-MSBL in this case. Finally, for BR signals, Fig. 10 shows a distinct advantage for EM-GM-AMP and EM-GM-AMP-MOS over the other algorithms, including EM-BG-AMP, due to the formers' ability to accurately model the BR signal prior. In particular, for M/N ≥ 0.36, EM-GM-AMP-MOS reduces the NMSE by 10 dB relative to the best of the other algorithms (which was either EM-BG-AMP or T-MSBL depending on the value of M/N ) and reduces the M/N breakpoint from 0.38 down to 0.35.

B. Noisy Sparse Signal Recovery
To investigate each algorithm's robustness to AWGN, we plotted the NMSE attained in the recovery of BR signals with N = 1000, M = 500, and K = 100 as a function of SNR in Fig. 11, where each point represents an average over R = 100 problem realizations, where in each realization we drew an A with i.i.d N (0, M −1 ) elements, an AWGN noise vector, and a random signal vector. All algorithms were under the same conditions as those reported previously, except that T-MSBL used noise=small when SNR > 22dB and noise=mild when SNR ≤ 22 dB, as recommended in [34]. From Fig. 11, we see that the essential behavior observed in the fixed-SNR BR plot Fig. 10 holds over a wide range of SNRs. In particular, Fig. 11 shows that EM-GM-AMP and EM-GM-AMP-MOS yield significantly lower NMSE than all other algorithms over the full SNR range, while EM-BG-AMP and T-MSBL yield the second lowest NMSE (also matched by BCS for SNRs between 30 and 40 dB). Note, however, than T-MSBL must be given some knowledge about the true noise variance in order to perform well [34], unlike the proposed algorithms.

C. Heavy-Tailed Signal Recovery
In many applications of compressive sensing, the signal to be recovered is not perfectly sparse, but instead contains a few large coefficients and many small ones. While the literature often refers to such signals as "compressible," there are many real-world signals that do not satisfy the technical definition of compressibility (see, e.g., [35]), and so we refer to such signals more generally as "heavy tailed." To investigate algorithm performance for these signals, we first consider an i.i.d Student's-t signal, with prior pdf under the (non-compressible) rate q = 1.67, which has been shown to be an excellent model for wavelet coefficients of natural images [35]. For such signals, Fig. 12 plots NMSE versus the number of measurements M for fixed N = 1000, SNR = 25 dB, and an average of R = 500 realizations, where in each realization we drew an A with i.i.d N (0, M −1 ) elements, an AWGN noise vector, and a random signal vector. Figure 12 shows both variants of EM-GM-AMP (here run in "heavy-tailed" mode) outperforming all other algorithms under test. 15 We have also verified (in experiments not shown here) that "heavy-tailed" EM-GM-AMP exhibits similarly good performance with other values of the Student's-t rate parameter q, as well as for i.i.d centered Cauchy signals.
To investigate the performance for positive heavy-tailed signals, we conducted a similar experiment using i.i.d lognormal x, generated using the distribution with location parameter µ = 0 and scale parameter σ 2 = 1. It may be interesting to notice that, with the perfectly sparse signals examined in Figs. 8-10, SL0 and SPGL1 performed relatively poorly, the relevance-vector-machine (RVM)-based approaches (i.e., BCS, T-MSBL) performed relatively well, and the greedy approaches (OMP and SP) performed inbetween. With the heavy-tailed signals in Figs. 12-13, it is more difficult to see a consistent pattern. For example, with the Student's-t signal, the greedy approaches performed the worse, the RVM approaches were in the middle, and SL0 and SPGL1 performed very well. But with the log-normal signal, the situation was very different: the greedy approaches performed very well, SPGL1 performed moderately well, but SL0 and the RVM approaches performed very poorly.
In conclusion, for all of the many signal types tested above, the best recovery performance came from EM-GM-AMP and its MOS variant. We attribute this behavior to EM-GM-AMP's ability to tune itself to the signal (and in fact the realization) at hand.

D. Runtime and Complexity Scaling with N
Next we investigated how complexity scales with signal length N by evaluating the runtime of each algorithm on a typical personal computer. For this, we fixed K/N = 0.1, M/N = 0.5, SNR = 25 dB and varied the signal length N . Figure 14 shows the runtimes for noisy recovery of a Bernoulli-Rademacher signal, while Fig. 15 shows the corresponding NMSEs. In these plots, each datapoint represents an average over R = 50 realizations. The algorithms that we tested are the same ones that we described earlier. However, to fairly evaluate runtime, we configured some a bit differently than before. In particular, for genie-tuned SPGL1, in order to yield a better runtime-vs-NMSE tradeoff, we reduced the tolerance grid (recall footnote 14) to σ 2 ∈ {0.6, 0.8, . . . , 1.4}×M ψ and turned off debiasing. For OMP and SP, we used the fixed support size K lasso = M ρ SE ( M N ) rather than searching for the size that minimizes NMSE over a grid of 10 hypotheses, as before. Otherwise, all algorithms were run under the suggested defaults, with T-MSBL run under noise=small and EM-GM-AMP run in "sparse" mode.
The complexities of the proposed EM-GM-AMP methods are dominated by one matrix multiplication by A and A T per iteration. Thus, when these matrix multiplications are explicitly implemented and A is dense, the total complexity of EM-GM-AMP should scale as O(M N ). This scaling is indeed visible in the runtime curves of Fig. 14. There, O(M N ) becomes O(N 2 ) since the ratio M/N was fixed, and the horizontal axis plots N on a logarithmic scale, so that this complexity scaling manifests, at sufficiently large values of N , as a line with slope 2. Figure 14 confirms that genietuned SPGL1 also has the same complexity scaling, albeit with longer overall runtimes. Meanwhile, Fig. 14 shows T-MSBL, BCS, SL0, OMP, and SP exhibiting a complexity scaling of O(N 3 ) (under fixed K/N and M/N ), which results in orders-of-magnitude larger runtimes for long signals (e.g., N ≥ 10 4 ). With short signals (e.g., N < 1300), though, OMP, SP, SL0, and SPGL1 are faster than EM-GM-AMP. Finally, Fig. 15 verifies that, for most of the algorithms, the NMSEs are relatively insensitive to signal length N when the undersampling ratio M/N and sparsity ratio K/M are both fixed, although the performance of EM-GM-AMP improves with N (which is not surprising in light of AMP's largesystem-limit optimality properties [13]) and the performance of BCS degrades with N .
Both the proposed EM-GM-AMP methods and SPGL1 can exploit the case where multiplication by A and A T is implemented using a fast algorithm like the fast Fourier transform (FFT) 16 , which reduces the complexity to O (N log N ), and avoids the need to store A in memory-a potentially serious problem when M N is large. The dashed lines in Figs. 14-15 (labeled "fft") show the average runtime and NMSE of the proposed algorithms and SPGL1 in case that A was a randomly row-sampled FFT. As expected, the runtimes are dramatically reduced. While EM-BG-AMP retains its place as the fastest algorithm, SPGL1 now runs 1.5× faster than EM-GM-AMP (at the cost of 14 dB higher NMSE). The MOS version of EM-GM-AMP yields slightly better NMSE, but takes ≈ 2.5 times as long to run as the fixed-L version.

E. Example: Compressive Recovery of Audio
As a practical example, we experimented with the recovery of an audio signal from compressed measurements. The full length-81920 audio signal was first partitioned into T blocks {u t } T t=1 of length N . Noiseless compressed measurements y t = Φu t ∈ R M were then collected using M = N/2 samples per block. Rather than reconstructing u t directly from y t , we first reconstructed 17 the transform coefficients x t = Ψ T u t , using the (orthogonal) discrete cosine transform (DCT) Ψ ∈ R N ×N , and later reconstructed u t via u t = Ψx t . Our effective sparse-signal model can thus be written as y t = Ax t with A = ΦΨ. We experimented with two types of measurement matrix Φ: i.i.d zero-mean Gaussian and random selection (i.e., containing rows of the identity matrix selected uniformly at 17 Although one could exploit additional structure among the multipletimestep coefficients {xt} T t=1 for improved recovery (e.g., sparsity clustering in the time and/or frequency dimensions, as well as amplitude correlation in those dimensions) as demonstrated in [36], such techniques are outside the scope of this paper. random), noting that the latter allows a fast implementation of A and A T . Table III shows the resulting time-averaged NMSE, i.e., TNMSE 1 T T t=1 ||u t −û t || 2 /||u t || 2 , and total runtime achieved by the previously described algorithms at block lengths N = 1024, 2048, 4096, 8192, which correspond to T = 80, 40, 20, 10 blocks, respectively. The numbers reported in the table represent an average over 50 realizations of Φ. For these experiments, we configured the algorithms as described in Section IV-C for the heavy-tailed experiment except that, for genie-SPGL1, rather than using ψ = 0, we used ψ = 10 −6 for the tolerance grid (recall footnote 14) because we found that this value minimized TNMSE and, for T-MSBL, we used the setting prune_gamma = 10 −12 as recommended in a personal correspondence with the author. For certain combinations of algorithm and blocklength, excessive runtimes prevented us from carrying out the experiment, and thus no result appears in the table. Table III shows that, for this audio experiment, the EM-GM-AMP methods and SL0 performed best in terms of TNMSE. As in the synthetic examples presented earlier, we attribute EM-GM-AMP's excellent TNMSE to its ability to tune itself to whatever signal is at hand. As for SL0's excellent TNMSE, we reason that it had the good fortune of being particularly well-tuned to this audio signal, given that it performed relatively poorly with the signal types used for Figs. 8-11 and Fig. 13. From the runtimes reported in Table III, we see that, with i.i.d Gaussian Φ and the shortest block length (N = 1024), genie-OMP is by far the fastest, whereas the EM-GM-AMP methods are the slowest. But, as the block length grows, the EM-GM-AMP methods achieve better and better runtimes as a consequence of their excellent complexity scaling, and eventually EM-BG-AMP and fixed-L EM-GM-AMP become the two fastest algorithms under test (as shown with i.i.d Gaussian Φ at N = 8192). For this audio example, the large-block regime may be the more important, because that is where all algorithms give their smallest TNMSE. Next, looking at the runtimes under random-selection Φ, we see dramatic speed improvements for the EM-GM-AMP methods and SPGL1, which were all able to leverage Matlab's fast DCT. In fact, the total runtimes of these four algorithms decrease as N is increased from 1024 to 8192. We conclude by noting that EM-BG-AMP (at N = 8192 with random selection Φ) achieves the fastest runtime in the entire table while yielding a TNMSE that is within 1.3 dB of the best value in the entire table. Meanwhile, fixed-L EM-GM-AMP (at N = 8192 with random selection Φ) gives TNMSE only 0.3 dB away from the best in the entire table with a runtime of only about twice the best in the entire table. Finally, the best TNMSEs in the entire table are achieved by EM-GM-AMP-MOS (at N = 8192), which takes ≈ 2.5 times as long to run as its fixed-L counterpart.

V. CONCLUSIONS
Those interested in practical compressive sensing face the daunting task of choosing among literally hundreds of signal reconstruction algorithms (see, e.g., [37]). In testing these algorithms, they are likely to find that some work very well  with particular signal classes, but not with others. They are also likely to get frustrated by those algorithms that require the tuning of many parameters. Finally, they are likely to find that some of the algorithms that are commonly regarded as "very fast" are actually very slow in high-dimensional problems. Meanwhile, those familiar with the theory of compressive sensing know that the workhorse LASSO is nearly minimax optimal, and that its phase transition curve is robust to the nonzero-coefficient distribution of sparse signals. However, they also know that, for most signal classes, there is a large gap between the MSE performance of LASSO and that of the MMSE estimator derived under full knowledge of the signal and noise statistics [11]. Thus, they may wonder whether there is a way to close this gap by designing a signal reconstruction algorithm that both learns and exploits the signal and noise statistics.
With these considerations in mind, we proposed an empirical Bayesian approach to compressive signal recovery that merges two powerful inference frameworks: expectation maximization (EM) and approximate message passing (AMP). We then demonstrated-through a detailed numerical studythat our approach, when used with a flexible Gaussian-mixture signal prior, achieves a state-of-the-art combination of reconstruction error and runtime on a very wide range of signal and matrix types in the high-dimensional regime. However, certain non-zero-mean and super-Gaussian sensing matrices give our AMP-based method trouble. Making AMP robust to these matrices remains a topic of importance for future research.