Message Passing Algorithms for Compressed
Sensing: II. Analysis and Validation
David L. Donoho
Department of Statistics
Stanford University
Arian Maleki
Department of Electrical Engineering
Stanford University
Andrea Montanari
Department of Electrical Engineering
and Department of Statistics
Stanford University
Abstract—In a recent paper, the authors proposed a new class
of low-complexity iterative thresholding algorithms for recon-
structing sparse signals from a small set of linear measurements
[1]. The new algorithms are broadly referred to as AMP, for
approximate message passing. This is the second of two conference
papers describing the derivation of these algorithms, connection
with related literature, extensions of original framework, and
new empirical evidence.
This paper describes the state evolution formalism for analyz-
ing these algorithms, and some of the conclusions that can be
drawn from this formalism. We carried out extensive numerical
simulations to confirm these predictions. We present here a few
representative results.
I. GENERAL AMP AND STATE EVOLUTION
We consider the model
y = A s
o
+ w
o
, s
o
∈ R
N
, y, w
o
∈ R
n
, (1)
with s
o
a vector that is ‘compressible’ and w
o
a noise vector.
We will assume that the entries of w
o
are centered independent
gaussian random variables with variance v.
The general AMP (approximate message passing) algorithm
reads
x
t+1
= η
t
(x
t
+ A
∗
z
t
) , (2)
z
t
= y − Ax
t
+
1
δ
z
t−1
hη
0
t−1
(x
t−1
+ A
∗
z
t−1
)i , (3)
with initial condition x
0
= 0. Here, for a vector u =
(u
1
, . . . , u
N
) we write hui ≡
P
N
i=1
u
i
/N , and η
0
( · ; · )
indicates the derivative of η with respect to its first argument.
Further δ ≡ n/N and {η
t
( · )}
t≥0
is a sequence of scalar
non-linearities (see Section III), a typical example being soft
thresholding, which contracts its argument towards zero.
A. Structure of the Algorithm
This algorithm is interesting for its low complexity: its
implementation is dominated at each step by the cost of
applying A and A
∗
to appropriate vectors. In some important
settings, matrices A of interest can be applied to a vector
implicitly by a pipeline of operators requiring N log(N ) flops;
an example would be A whose rows are randomly chosen from
among the rows of a Fourier matrix; then Ax can be computed
by FFT and subsampling.
Even more, the algorithm is interesting for the message
passing term
1
δ
z
t−1
hη
0
t−1
(x
t−1
+A
∗
z
t−1
)i. Similar algorithms
without this term are common in the literature of so-called iter-
ative thresholding algorithms. As discussed in the companion
paper, the message passing term approximates the combined
effect on the reconstruction of the passing of nN messages in
the the full message passing algorithm.
The message passing term completely changes the statis-
tical properties of the reconstruction, and it also makes the
algorithm amenable to analysis by a technique we call State
Evolution. Such analysis shows that the algorithm converges
rapidly, much more rapidly than any known result for the
IST algorithm. Furthermore, it allows us to make a variety
of theoretical predictions about performance characteristics of
the algorithm which are much stronger than any predictions
available for competing methods.
B. State Evolution
In the following we will assume that the columns of A are
normalized to unit Euclidean length. We define the effective
variance
σ(x
t
)
2
≡ v +
1
Nδ
||x
t
− s
0
||
2
2
. (4)
The effective variance combines the observational variance
v with an additional term
1
Nδ
||x
t
− s
0
||
2
2
that we call the
interference term. Notice that v is merely the squared recon-
struction error of the naive ‘matched filter’ for the case where
s
0
contains all zeros and a single nonzero in a given position
i and the matched filter is just the i-th column of A.
The interference term measures the additional error in
estimating a single component of s
o,i
that is caused by the
many small errors in other components j 6= i. The formula
states that the effective variance at iteration t is caused by the
observational noise (invariant across iteration) and the current
errors at iteration t (changing from iteration to iteration). The
interference concept is well known in digital communications,
where phrases like mutual access interference are used for
what is algebraically the same phenomenon.
We will let bσ
t
denote any estimate of σ
t
, and we will assume
that bσ
t
≈ σ
t
; see [1] for more careful discussion. Suppose that
the nonlinearity takes the form η
t
( · ) = η( · ; θ
t
) where θ is
a tuning parameter, possibly depending on bσ
t
; see below for
more. Let F denote the collection of CDFs on R and F be the
CDF of s
0
(i). Define the MSE map Ψ : R
+
× R
3
× F 7→ R
+
arXiv:0911.4222v1 [cs.IT] 22 Nov 2009
by
Ψ(σ; v, δ, θ
t
, F ) = v +
1
δ
E
n
η
t
X + σ Z
− X
2
o
where X has distribution F and Z ∼ N (0, 1) is independent
of X. We suppose that a rule Θ(σ; v, δ, θ, F ) for the update
of θ
t
is also known.
Definition I.1. The state is a 5-tuple S = (σ; v, δ, θ, F ); state
evolution is the evolution of the state by the rule
(σ
2
t
; v, δ, θ
t
, F ) 7→ (Ψ(σ
2
t
); v, δ, θ
t+1
, F )
t 7→ t + 1
As the parameters (v, δ, F ) remain fixed during evolution, we
usually omit mention of them and think of state evolution
simply as the iterated application of Ψ and Θ:
σ
2
t
7→ σ
2
t+1
≡ Ψ(σ
2
t
)
θ
t
7→ θ
t+1
≡ Θ(S
t
)
t 7→ t + 1
The initial state is taken to have σ
2
0
= v + ||s
0
||
2
2
/N δ.
As described, State Evolution is a purely analytical con-
struct, involving sequential application of rules Ψ and Θ. The
crucial point is to know whether this converges to a fixed point,
and to exploit the properties of the fixed point. We expect that
such properties are reflected in the properties of the algorithm.
To make this precise, we need further notation.
Definition I.2. State-Conditional Expectation. Given a func-
tion ζ : R
4
7→ R, its expectation in state S
t
is
E(ζ|S
t
) = E
ζ(U, V, W, η(U + V + W ))
,
where U ∼ F , V ∼ N(0, v) and W ∼ N(0, σ
2
t
− v).
Different choices of ζ allow to monitor the evolution of
different metrics under the AMP algorithm. For instance, ζ =
(u − x)
2
corresponds to the mean square error (MSE). The
False Alarm Rate is tracked by ζ = 1
{η(v+w)6=0}
and the
Detection Rate by ζ = 1
{η(u+v+w)6=0}
.
Definition I.3. Large-System Limits. Let ζ : R
4
7→ R be
a function of real 4-tuples (s, u, w, x). Suppose we run the
iterative algorithm A for a sequence of problem sizes (n, N)
at a the value (v, δ, F ) of underlying implicit parameters,
getting outputs x
t
, t = 1, 2, 3, . . . The large-system limit
ls.lim(ζ, t, A) of ζ at iteration t is
ls.lim(ζ, t, A) = p.lim
N→∞
hζ(s
o,i
, u
t,i
, w
o,i
, x
t,i
)i
N
,
where h · i
N
denotes the uniform average over i ∈
{1, . . . , N} ≡ [N ], and p.lim denotes limit in probability.
Hypothesis I.4. Correctness of State Evolution for AMP.
Run an AMP algorithm for t iterations with implicit state
variables v, δ, F . Run state evolution, obtaining the state S
t
at
time t. Then for any bounded continuous function ζ : R
4
7→ R
of the real 4-tuples (s, u, w, x), and any number of iterations
t,
1) The large-system limit ls.lim(ζ, t, A) exists for the ob-
servable ζ at iteration t.
2) This limit coincides with the expectation E(ζ|S
t
) com-
puted at state S
t
.
State evolution, where correct, allows us to predict the
performance of AMP algorithms and tune them for optimal
performance. In particular, SE can help us to choose the non-
linearities {η
t
} and their tuning. The objective of the rest of
this paper is twofold: (1) Provide evidence for state evolution;
(2) Describe some guidelines towards the choice of the non-
linearities {η
t
}.
II. AMP-BASED ALGORITHMS
Already in [1] we showed that a variety of algorithms
can be generated by varying the choice of η. We begin with
algorithms based on soft thresholding. Here η
t
(x) = η(x; θ
t
)
is given by the soft threshold function
η(x; θ) =
x − θ if θ < x,
0 if −θ ≤ x ≤ θ,
x + θ if x < −θ.
(5)
This function shrinks its argument towards the origin. Several
interesting AMP-Based algorithms are obtained by varying the
choice of the sequence {θ
t
}
t∈N
.
A. AMP.M(δ)
The paper [1] considered the noiseless case v = 0 where
the components of s
o
are iid with common distribution F that
places all but perhaps a fraction = ρ(δ) · δ, ρ ∈ (0, 1) of its
mass at zero. That paper proposed the choice
θ
t
= τ(δ)bσ
t
. (6)
where an explicit formula for τ(δ) is derived in the online
supplement [2]. As explained in that supplement, this rule has
a minimax interpretation, namely, to give the smallest MSE
guaranteed across all distributions F with mass at zero larger
than or equal to 1 − .
B. AMP.T(τ )
Instead of taking a worst case viewpoint, we can think of
specifically tuning for the case at hand. Consider general rules
of the form:
θ
t
= τ bσ
t
. (7)
Such rules have a very convenient property for state evolution;
namely, if we suppose that bσ
t
≡ σ
t
, we can redefine the
state as (σ
2
t
; v, δ, τ, F ), with (v, δ, τ, F ) invariant during the
iteration, and then the evolution is effectively one-dimensional:
σ
2
t
7→ σ
2
t+1
≡ Ψ(σ
2
t
). The dynamics are then very easy to
study, just by looking for fixed points of a scalar function Ψ.
(This advantage is also shared by AMP.M(δ), of course).
While the assumption bσ
t
≡ σ
t
does not hold, strictly
speaking, at any finite size, it will hold asymptotically in the
large system limit for many good estimators of the effective
variance.
It turns out that, depending on F and δ, different values
of τ lead to very different performance characteristics. It is
natural to ask for the fixed value τ = τ
∗
(v, δ, F ) which, under
state evolution gives the smallest equilibrium MSE. We have
developed software to compute such optimal tuning; results
are discussed in [5].
C. AMP.A(λ)
In much current work on compressed sensing, it is desired
to solve the `
1
-penalized least squares problem
minimize
1
2
ky − Axk
2
2
+ λkxk
1
. (8)
In different fields this has been called Basis Pursuit denoising
[6] or Lasso [7]. Large scale use of general convex solvers is
impractical when A is of the type interesting from compressed
sensing, but AMP-style iterations are practical. And, surpris-
ingly an AMP-based algorithm can effectively compute the
solution by letting the threshold ‘float’ to find the right level
for solution of the above problem. The threshold recursion is:
θ
t+1
= λ +
θ
t
δ
hη
0
(x
t
+ A
∗
z
t
; θ
t
)i . (9)
D. AMP.0
It can also be of interest to solve the `
1
-minimization
problem
min
x
kxk
1
subject to y = Ax. (10)
This has been called Basis Pursuit [6] in the signal processing
literature. While formally it can be solved by linear program-
ming, standard linear program codes are far too slow for many
of the applications interesting to us.
This is formally the λ = 0 case of AMP.A(λ). In fact it can
be advantageous to allow λ to decay with the iteration number
θ
t+1
= λ
t
+
θ
t
δ
hη
0
(x
t
+ A
∗
z
t
; θ
t
)i . (11)
Here, we let λ
t
↓ 0 as t → ∞.
E. Other Nonlinearities
The discussion above has focused entirely on soft thresh-
olding, but both the AMP algorithm and SE formalism make
perfect sense with many other nonlinearities. Some case of
specific interest include
• The Bayesian conditional mean: η(x) = E{s
0
|s
0
+ U +
V = x}, where U and V are just as in Definition I.2. This
is indeed discussed in the companion paper [3], Section
V.
• Scalar nonlinearities associated to various nonconvex op-
timization problems, such as minimizing `
p
pseudonorms
for p < 1.
III. CONSEQUENCES OF STATE EVOLUTION
A. Exponential Convergence of the Algorithm
When State Evolution is correct for an AMP-type algorithm,
we can be sure that the algorithm converges rapidly to its
limiting value – exponentially fast. The basic point was
shown in [1]. Suppose we are considering either AMP.M(δ)
or AMP.T(τ ). In either case, as explained above, the state
evolution is effectively one-dimensional. Then the following
is relevant.
Definition III.1. Stable Fixed Point. The Highest Fixed Point
of the continuous function Ψ is
HFP(Ψ) = sup{m : Ψ(m) ≥ m}.
The stability coefficient of the continuously differentiable func-
tion Ψ is
SC(Ψ) =
d
dm
Ψ(m)
m=HFP(Ψ)
.
We say that HFP(Ψ) is a stable fixed point if 0 ≤ SC(Ψ) < 1.
Let µ
2
(F ) =
R
x
2
dF denote the second-moment functional
of the CDF F .
Lemma III.2. Let Ψ( · ) = Ψ( · ; v, δ, F ). Suppose that
µ
2
(F ) > HFP(Ψ). The sequence of iterates σ
2
t
defined by
starting from σ
2
0
= µ
2
(F ) and σ
2
t+1
= Ψ(σ
2
t
) converges:
σ
2
t
→ HFP(Ψ), t → ∞.
Suppose that the stability coefficient 0 < SC(Ψ) < 1. Then
(σ
2
t
− HFP(Ψ)) ≤ SC(Ψ)
t
· (µ
2
(F ) − HFP(Ψ)).
In short, when F and v are such that the highest fixed point
is stable, state evolution converges exponentially fast to that
fixed point.
Other iterative thresholding algorithms have theoretical
guarantees which are far weaker. For example, FISTA [8] has
a theoretical guarantee of O(1/t
2
), while SE evolution implies
O(exp(−ct)).
B. Phase Transitions For `
1
minimization
Consider the special setting where the noise is absent w
o
=
0 and the object s
o
obeys a strict sparsity condition; namely
the distribution F places a fraction ≥ 1 − of its mass at the
origin; and thus, if s
o
is iid F , approximately N · (1 − ) of
its entries are exactly zero.
A phase transition occurs in this setting when using `
1
minimization for reconstruction. Namely, if we write = ρ · δ
then there is a critical value ρ(δ) such that, for < ρ(δ)· δ, `
1
minimization correctly recovers s
o
, while for > ρ(δ) · δ, `
1
minimization fails to correctly recover s
o
, with probability ap-
proaching one in the large size limit. State Evolution predicts
this phenomenon, because, for < ρ
SE
(δ)·δ, the highest fixed
point is at σ
2
t
= 0, while above this value, the highest fixed
point is at σ
2
t
> 0. Previously, the exact critical value ρ(δ) at
which this transition occurs was computed by combinatorial
geometry, with a rigorous proof; however, it was shown in [1]
that the algorithm AMP.M(δ) has ρ(δ) = ρ
SE
(δ), validating
the correctness of SE.
0 10 20 30 40
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
MSE on the non−zero elements
iteration
empirical
theoretical
0 10 20 30 40
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
MSE on zero elements
iteration
empirical
theoretical
0 10 20 30 40
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Missed Detection Rate
iteration
empirical
theoretical
0 10 20 30 40
0.6985
0.699
0.6995
0.7
0.7005
0.701
0.7015
0.702
False Alarm Rate
iteration
empirical
theoretical
Fig. 1. Observables versus iteration, and predictions by state evolution. Panels
(a)-(d): MSENZ, MSE, MDR, FAR. Curve in red: theoretical prediction. Curve
in blue: mean observable. For this experiment, N = 5000, δ = n/N = .3.
F = 0.955δ
0
+ 0.045δ
1
C. Operating Characteristics of `
1
penalized Least-squares.
State evolution predicts the following relationships between
AMP.T(τ) and BPDN(λ). AMP.T(τ) has, according to SE,
for its large-t limit an equilibrium state characterized by its
equilibrium noise plus interference level σ
∞
(τ). In that state
AMP.T(τ) uses an equilibrium threshold θ
∞
(τ). Associated
to this equilibrium NPI and Threshold, there is an equilibrium
detection rate
EqDR(τ) = P{η(U + V + W ; θ
∞
) 6= 0}
where U ∼ F, V is N (0, v) and W is N(0, σ
2
∞
− v), with
U,V ,W independent. Namely, for all sufficiently large τ (i.e
τ > τ
0
(δ, F, v)) we have
λ = (1 − EqDR(τ)/δ) · θ
∞
(τ);
this creates a one-one relationship λ ↔ τ(λ; v, δ, F ) cal-
ibrating the two families of procedures. SE predicts that
observables of the `
1
-penalized least squares estimator with
penalty λ will agree with the calculations of expectations for
AMP.T(τ(λ; v, δ, F )) made by state evolution.
IV. EMPIRICAL VALIDATION
The above-mentioned consequences of State Evolution can
be tested as follows. In each case, we can use SE to make
a fixed prediction in advance of an experiment and then we
can run a simulation experiment to test the accuracy of the
prediction.
A. SE Predictions of Dynamics of Observables
Exponential convergence of AMP-based algorithms is
equivalent to saying that a certain observable – Mean-squared
error of reconstruction – decays exponentially in t. This is but
one observable of the algorithm’s output; and we have tested
not only the SE predictions of MSE but also the SE predictions
of many other quantities.
In Figure 1 we present results from an experiment with
signal length N = 5000, noise level v = 0, indeterminacy
δ = n/N = 0.30 and sparsity level = 0.045. The distribution
F places 95.5% of its mass at zero and 4.5% of its mass at
1. the fit between predictions and observations is extremely
good – so much so that it is hard to tell the two curves apart.
For more details, see [2].
B. Phase Transition Calculations
Empirical observations of Phase transitions of `
1
minimiza-
tion and other algorithms have been made in [4], [9], and we
follow a similar procedure. Specifically, to observe a phase
transition in the performance of a sparsity-seeking algorithm,
we perform 200 reconstructions on randomly-generated prob-
lem instances with the same underlying situation (v = 0, δ,
F ) and we record the fraction of successful reconstructions in
that situation. We do this for each member of a large set of
situations by varying the undersampling ratio δ and varying
sparsity of F . More specifically, we define a (δ, ρ) phase
diagram [0, 1] and consider a grid of sites in this domain with
δ = .05, .10, . . . and ρ = .03, .06, . . . , .99. For each δ, ρ pair
in this grid, we generate random problem instances having a
k-sparse solution s
0
, i.e. a vector having k ones and n − k
zeros; here k = ρ · δ · N.
Defining success as exact recovery of s
0
to within a small
fixed error tolerance, we define the empirical phase transition
as occurring at the ρ value where the success fraction drops
below 50%. For more details, see [2].
Figure 2 depicts the theoretical phase transition predicted
by State Evolution as well as the empirical phase transition
of AMP.M(δ) and a traditional iterative soft thresholding
algorithm. In this figure N = 1000, and AMP.M(δ) was run
for T = 1, 000 iterations. One can see that empirical phase
transition of AMP.M(δ) matches closely the state evolution
prediction. One can also see that the empirical phase transition
of iterative soft thresholding, without the message passing
term, is substantially worse than that for the AMP-based
method with the message passing term.
C. Operating Characteristics of `
1
penalized Least-squares
The phase transition study gives an example of SE’s accu-
racy in predicting AMP-based algorithms in a strictly sparse
setting, i.e. where only a small fraction of entries in s
0
are
nonzero. For a somewhat different example, we consider the
generalized Gaussian family, i.e. distribution functions F
α
with
densities
f
α
(x) = exp(−|x|
α
)/Z
α
.
In the case α = 1 there is a very natural connection with `
1
-
minimization algorithms, which then become MAP estimation
schemes. In the case α = 1, an iid realization from f
α
,
properly rescaled to unit `
1
norm, will be uniformly distributed
on the surface of the `
1
ball, and in that sense this distribution
samples all of the `
1
ball, unlike the highly sparse distributions
used in the phase transition study, which sample only the low-
dimensional faces. When α < 1, the sequence is in a sense
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
δ
ρ
Comparison of Different Algorithms
IST
L
1
AMP
Fig. 2. Phase transitions of reconstruction algorithms. Blue Curve:
Phase Transition predicted by SE; Red Curve empirical phase transition
for AM P M(δ) as observed in simulation; Green Curve, empirical phase
transition for Iterative Soft Thresholding as observed in simulation.
more sparse than when α = 1. The case α = .7 has been found
useful in modelling wavelet coefficients of natural images.
We considered exponents α ∈ {0.35, 0.50, 0.65, 0.75, 1.0}.
At each such case we considered incompleteness ratios δ ∈
{0.1, 0.2, 0.3, 0.4, 0.5}. The set of resulting (α, δ) pairs gives
a collection of 25 experimental conditions. At each such
experimental condition, we considered 5 or so different values
of λ for which SE-predicted MSE’s were available. In total,
simulations were run for 147 different combinations of α, δ
and λ. At each such combination, we randomly generated 200
problem instances using the problem specification, and then
computed more than 50 observables of the solution. In this
subsection, we used N = 500.
To solve an instance of problem (8) we had numerous
options. Rather than a general convex optimizer, we opted to
use the LARS/LASSO algorithm.
Figure 3 shows a scatterplot comparing MSE values for
the LARS/LASSO solution of (8) with predictions by State
Evolution, as decribed in section III.C. Each data point cor-
responds to one experimental combination of α, δ, λ, and the
datapoint presents the median MSE across 200 simulations
under that combination of circumstances. Even though the
observed MSE’s vary by more than an order of magnitude, it
will be seen that the SE predictions track them accurately. It
should be recalled that the problem size here is only N = 500,
and that only 200 replications were made at each experimental
situation. In contrast, the SE prediction is designed to match
large-system limit. In a longer paper, we will consider a much
wider range of observables and demonstrate that, at larger
problem sizes N , we get successively better fits between
observables and their SE predictions.
ACKNOWLEDGMENT
The authors would like to thank NSF for support in grants
DMS-05-05303 and DMS-09-06812, and CCF-0743978 (CA-
0 5 10 15
0 5 10 15
MSE(Ell1PenEll2) and SE pred.; GenGauss Vary delta & p
1D SE prediction
observed MSE,
Fig. 3. Mean-squared Error for `
1
-penalized Least-squares estimate versus
predicted error according to State Evolution.
REER) and DMS-0806211 (AM).
REFERENCES
[1] D. L. Donoho, A. Maleki and A. Montanari, “Message-passing algorithms
for compressed sensing,” PNAS 2009 106:18914-18919; November 10,
2009
[2] D. L. Donoho, A. Maleki and A. Montanari, on-line supplement to
accompany [1]. See also appendix to arXiv:0907.3574v1.
[3] D. L. Donoho, A. Maleki, A. Montanari, “Message passing algorithms for
compressed sensing: I. motivation and construction,” Proc. ITW, Cairo,
Egypt, 2010.
[4] A. Maleki and D. L. Donoho. Optimally tuned iterative thresholding
algorithms for compressed sensing. arXiv:0909.0777, 2009
[5] D. L. Donoho, A. Maleki and A. Montanari, “Theoretical Calculations of
LASSO operating Characteristics in the case p > n” Manuscript.
[6] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition
by basis pursuit,” SIAM Journal on Scientific Computing, vol. 20, pp.
33-61, 1998.
[7] R. Tibshirani, “ Regression Shrinkage and Selection via the Lasso,”
Journal of the Royal Statistical Society. Series B (Methodological), Vol.
58, No. 1. (1996), pp. 267-288, 1996.
[8] A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding
Algorithm for Linear Inverse Problems,” SIAM journal on imaging
sciences, vol. 2, no. 1, pp. 183-202, 2009.
[9] D. L. Donoho and J. Tanner. Observed universality of phase transitions in
high-dimensional geometry, with implications for modern data analysis
and signal processing. Phil. Trans. Roy. Soc. A, 2009, 367, 4273-4293.
[10] D. L. Donoho, A. Maleki and A. Montanari, “Minimax MSE over `
p
balls in Compressed Sensing” Manuscript.
