Compressive Imaging using Approximate Message Passing and a Markov-Tree Prior

While exact computation of $p(\boldsymbol{\theta}\,|\,\boldsymbol{y})$ is computationally prohibitive, the marginal posteriors $\{p(\theta_{n}\,|\,\boldsymbol{y})\}$ can be efficiently approximated using loopy belief propagation (LBP) [10] on the factor graph of Fig. 2, which uses round nodes to denote variables and square nodes to denote the factors in (10). In doing so, we also obtain the marginal posteriors $\{p(s_{n}\,|\,\boldsymbol{y})\}$. For now, we treat statistical parameters $\rho,\pi_{-1}^{1},\pi_{0}^{1},\{\rho_{j},\pi_{j}^{11},\pi_{j}^{00}\}$, as if they were fixed and known, and we detail the procedure by which they are learned in Section III-D.

In LBP, messages are exchanged between the nodes of the factor graph until they converge. Messages take the form of pdfs (or pmfs), and the message flowing to/from a variable node can be interpreted as a local belief about that variable. According to the sum-product algorithm [12, 13] the message emitted by a variable node along a given edge is (an appropriate scaling of) the product of the incoming messages on all other edges. Meanwhile, the message emitted by a function node along a given edge is (an appropriate scaling of) the integral (or sum) of the product of the node’s constraint function and the incoming messages on all other edges, where the integration (or summation) is performed over all variables other than the one directly connected to the edge along which the message travels. When the factor graph has no loops, exact marginal posteriors result from two (i.e., forward and backward) passes of the sum-product algorithm [12, 13]. When the factor graph has loops, however, exact inference is known to be NP hard [17] and so LBP is not guaranteed to produce correct posteriors. Still, LBP has shown state-of-the-art performance in many applications, such as inference on Markov random fields [18], turbo decoding [19], LDPC decoding [20], multiuser detection [21], and compressive sensing [22, 14, 23, 15].

With loopy belief propagation, there exists some freedom in how messages are scheduled. In this work, we adopt the “turbo” approach recently proposed in [11]. For this, we split the factor graph in Fig. 2 along the dashed line and obtain the two decoupled subgraphs in Fig. 3. We then alternate between belief propagation on each of these two subgraphs, treating the likelihoods on $\{s_{n}\}$ generated from belief propagation on one subgraph as priors for subsequent belief propagation on the other subgraph. We now give a more precise description of this turbo scheme, referring to one full round of alternation as a “turbo iteration.” In the sequel, we use $\nu_{A\to B}^{\scriptscriptstyle(t)}(.)$ to denote the message passed from node $A$ to node $B$ during the $t^{\mathrm{th}}$ turbo iteration.

The procedure starts at $t=1$ by setting the “prior” pmfs $\{h_{n}^{\scriptscriptstyle(1)}(.)\}$ in accordance with the apriori activity rates $\Pr\{s_{n}=1\}$ described in Section II. LBP is then iterated (to convergence) on the left subgraph in Fig. 3, finally yielding the messages $\{\nu_{f_{n}\to s_{n}}^{\scriptscriptstyle(1)}(.)\}$. We note that the message $\nu_{f_{n}\to s_{n}}^{\scriptscriptstyle(1)}(s_{n})$ can be interpreted as the current estimate of the likelihood²²2 In turbo decoding parlance, the likelihood $\nu_{f_{n}\to s_{n}}^{\scriptscriptstyle(t)}(s_{n})$ would be referred to as the “extrinsic” information about $s_{n}$ produced by the left “decoder”, since it does not directly involve the corresponding prior $h_{n}^{\scriptscriptstyle(t)}(s_{n})$. Similarly, the message $\nu_{h\to s_{n}}^{(t)}(s_{n})$ would be referred to as the extrinsic information about $s_{n}$ produced by the right decoder. on $s_{n}$, i.e., $p(\boldsymbol{y}\,|\,s_{n})$ as a function of $s_{n}$. These likelihoods are then treated as priors for belief propagation on the right subgraph, as facilitated by the assignment $d_{n}^{\scriptscriptstyle(1)}(.)=\nu_{f_{n}\to s_{n}}^{\scriptscriptstyle(1)}(.)$ for each $n$. Due to the tree structure of HMT, there are no loops in right subgraph (i.e., inside the “$h$” super-node in Fig. 3), and thus it suffices to perform only one forward-backward pass of the sum-product algorithm [12, 13]. The resulting leftward messages $\nu_{h\to s_{n}}^{(1)}(.)$ are subsequently treated as priors for belief propagation on the left subgraph at the next turbo iteration, as facilitated by the assignment $h_{n}^{(2)}(.)=\nu_{h\to s_{n}}^{(1)}(.)$. The process then continues for turbo iterations $t=2,3,4,\dots$, until the likelihoods converge or a maximum number of turbo iterations has elapsed. Formally, the turbo schedule is summarized by

	$\displaystyle d_{n}^{(t)}(s_{n})$	$\displaystyle=$	$\displaystyle\nu_{f_{n}\to s_{n}}^{(t)}(s_{n})$		(11)
	$\displaystyle h_{n}^{(t+1)}(s_{n})$	$\displaystyle=$	$\displaystyle\nu_{h\to s_{n}}^{(t)}(s_{n}).$		(12)

In the sequel, we refer to inference of $\{s_{n}\}$ using compressive-measurement structure (i.e., inference on the left subgraph of Fig. 3) as soft support-recovery (SSR) and inference of $\{s_{n}\}$ using HMT structure (i.e., inference on the right subgraph of Fig. 3) as soft support-decoding (SSD). SSR details are described in the next subsection.

We now discuss our implementation of SSR during a single turbo iteration $t$. Because the operations are invariant to $t$, we suppress the $t$-notation. As described above, SSR performs several iterations of loopy belief propagation per turbo iteration using the fixed priors $\lambda_{n}\triangleq h_{n}(s_{n}=1)$. This implies that, over SSR’s LBP iterations, the message $\nu_{f_{n}\to\theta_{n}}(.)$ is fixed at

$\displaystyle\nu_{f_{n}\to\theta_{n}}(\theta_{n})$

$\displaystyle=$

$\displaystyle\lambda_{n}\mathcal{N}(\theta_{n};0,\sigma_{n}^{2})+(1-\lambda_{n})\delta(\theta_{n}).$

(13)

The dashed box in Fig. 3 shows the region of the factor graph on which messages are updated during SSR’s LBP iterations. This subgraph can be recognized as the one that Donoho, Maleki, and Montanari used to derive their so-called approximate message passing (AMP) algorithm [14]. While [14] assumed an i.i.d Laplacian prior for $\boldsymbol{\theta}$, the approach for generic i.i.d priors was outlined in [23]. Below, we extend the approach of [23] to independent non-identical priors (as analyzed in [24]) and we detail the Bernoulli-Gaussian case. In the sequel, we use a superscript-$i$ to index SSR’s LBP iterations.

According to the sum-product algorithm, the fact that $\nu_{f_{n}\to\theta_{n}}(.)$ is non-Gaussian implies that $\nu^{i}_{\theta_{n}\to g_{m}}(.)$ is also non-Gaussian, which complicates the exact calculation of the subsequent messages $\nu^{i}_{g_{m}\to\theta_{n}}(.)$ as defined by the sum-product algorithm. However, for large $N$, the combined effect of $\{\nu^{i}_{\theta_{n}\to g_{m}}(.)\}_{n=1}^{N}$ at the $g_{m}$ nodes can be approximated as Gaussian using central-limit theorem (CLT) arguments, after which it becomes sufficient to parameterize each message $\nu^{i}_{\theta_{n}\to g_{m}}(.)$ by only its mean and variance:

	$\displaystyle\mu_{mn}^{i}$	$\displaystyle\triangleq$	$\displaystyle\textstyle\int_{\theta_{n}}\theta_{n}\,\nu^{i}_{\theta_{n}\to g_{m}}(\theta_{n})$		(14)
	$\displaystyle v_{mn}^{i}$	$\displaystyle\triangleq$	$\displaystyle\textstyle\int_{\theta_{n}}(\theta_{n}-\mu^{i}_{mn})^{2}\,\nu^{i}_{\theta_{n}\to g_{m}}(\theta_{n}).$		(15)

Combining

$\displaystyle\prod_{q}\mathcal{N}(\theta;\mu_{q},v_{q})$

$\displaystyle\propto$

$\displaystyle\mathcal{N}\left(\theta;\frac{\sum_{q}\mu_{q}/v_{q}}{\sum_{q}1/v_{q}},\frac{1}{\sum_{q}1/v_{q}}\right)$

(16)

with $g_{m}(\boldsymbol{\theta})=\mathcal{N}(y_{m};\boldsymbol{a}^{T}_{m}\boldsymbol{\theta},\sigma^{2})$, the CLT then implies that

	$\displaystyle\nu^{i}_{g_{m}\to\theta_{n}}(\theta_{n})$	$\displaystyle\approx$	$\displaystyle\mathcal{N}\left(\theta_{n};\frac{z_{mn}^{i}}{A_{mn}},\frac{c^{i}_{mn}}{A_{mn}^{2}}\right)$		(17)
	$\displaystyle z^{i}_{mn}$	$\displaystyle\triangleq$	$\displaystyle\textstyle y_{m}-\sum_{q\neq n}A_{mq}\mu^{i}_{qm}$		(18)
	$\displaystyle c^{i}_{mn}$	$\displaystyle\triangleq$	$\displaystyle\textstyle\sigma^{2}+\sum_{q\neq n}A_{mq}^{2}v^{i}_{qm}.$		(19)

The updates $\mu^{i+1}_{mn}$ and $v^{i+1}_{mn}$ can then be calculated from

$\displaystyle\nu^{i+1}_{\theta_{n}\to g_{m}}(\theta_{n})$

$\displaystyle\propto$

$\displaystyle\textstyle\nu_{f_{n}\to\theta_{n}}(\theta_{n})\prod_{l\neq m}\nu^{i}_{g_{l}\to\theta_{n}}(\theta_{n}),$

(20)

where, using (16), the product term in (20) is

$\displaystyle\propto\mathcal{N}\left(\theta_{n};\frac{\sum_{l\neq m}A_{ln}z^{i}_{ln}/c^{i}_{ln}}{\sum_{l\neq m}A^{2}_{ln}z^{i}_{ln}/c^{i}_{ln}},\frac{1}{\sum_{l\neq m}A^{2}_{ln}/c^{i}_{ln}}\right).$

(21)

Assuming that the values $A_{ln}^{2}$ satisfy

$\displaystyle\sum_{l\neq m}A^{2}_{ln}\approx\sum_{l=1}^{M}A^{2}_{ln}\approx 1,$

(22)

which occurs, e.g., when $M$ is large and $\{A_{ln}\}$ are generated i.i.d with variance $1/M$, we have $c^{i}_{ln}\approx\textstyle c^{i}_{n}\triangleq\frac{1}{M}\sum_{m=1}^{M}c^{i}_{mn}$, and thus (20) is well approximated by

	$\displaystyle\nu^{i+1}_{\theta_{n}\to g_{m}}(\theta_{n})$	$\displaystyle\propto$	$\displaystyle\left(\lambda_{n}\mathcal{N}(\theta_{n};0,\sigma_{n}^{2})+(1-\lambda_{n})\delta(\theta_{n})\right)$		(23)
			$\displaystyle\mbox{}\times\mathcal{N}(\theta_{n};\xi_{nm}^{i},c_{n}^{i})$
	$\displaystyle\xi^{i}_{nm}$	$\displaystyle\triangleq$	$\displaystyle\textstyle\sum_{l\neq m}A_{ln}z^{i}_{ln}.$		(24)

In this case, the mean and variance of $\nu^{i+1}_{\theta_{n}\to g_{m}}(.)$ become

	$\displaystyle\mu_{nm}^{i+1}$	$\displaystyle=$	$\displaystyle\alpha_{n}(c_{n}^{i})\xi_{nm}^{i}/(1+\gamma_{nm}^{i})$		(25)
	$\displaystyle v_{nm}^{i+1}$	$\displaystyle=$	$\displaystyle\gamma_{nm}^{i}(\mu_{nm}^{i+1})^{2}+\mu_{nm}^{i+1}c_{n}^{i}/\xi_{nm}^{i}$		(26)
	$\displaystyle\gamma_{nm}^{i}$	$\displaystyle\triangleq$	$\displaystyle\beta_{n}(c_{n}^{i})\exp(-\zeta_{n}(c_{n}^{i})(\xi_{nm}^{i})^{2}),$		(27)

where

	$\displaystyle\alpha_{n}(c)$	$\displaystyle\triangleq$	$\displaystyle(c/\sigma_{n}^{2}+1)^{-1}$
	$\displaystyle\beta_{n}(c)$	$\displaystyle\triangleq$	$\displaystyle(-1+1/\lambda_{n})\sqrt{1+\sigma_{n}^{2}/c}$
	$\displaystyle\zeta_{n}(c)$	$\displaystyle\triangleq$	$\displaystyle(2c(1+c/\sigma_{n}^{2}))^{-1}.$

According to the sum-product algorithm, $\hat{p}^{i+1}(\theta_{n}\,|\,\boldsymbol{y})$, the posterior on $\theta_{n}$ after SSR’s $i^{\mathrm{th}}$-LBP iteration, obeys

$\displaystyle\hat{p}^{i+1}(\theta_{n}\,|\,\boldsymbol{y})$

$\displaystyle\propto$

$\displaystyle\textstyle\nu_{f_{n}\to\theta_{n}}(\theta_{n})\prod_{l=1}^{M}\nu^{i}_{g_{l}\to\theta_{n}}(\theta_{n}),$

(28)

whose mean and variance determine the $i^{\mathrm{th}}$-iteration MMSE estimate of $\theta_{n}$ and its variance, respectively. Noting that the difference between (28) and (20) is only the inclusion of the $m^{\mathrm{th}}$ product term, these MMSE quantities become

	$\displaystyle\mu_{n}^{i+1}$	$\displaystyle=$	$\displaystyle\alpha_{n}(c_{n}^{i})\xi_{n}^{i}/(1+\gamma_{n}^{i})$		(29)
	$\displaystyle v_{n}^{i+1}$	$\displaystyle=$	$\displaystyle\gamma_{n}^{i}(\mu_{n}^{i+1})^{2}+\mu_{n}^{i+1}c_{n}^{i}/\xi_{n}^{i}$		(30)
	$\displaystyle\xi^{i}_{n}$	$\displaystyle\triangleq$	$\displaystyle\textstyle\sum_{l=1}^{M}A_{ln}z^{i}_{ln}$		(31)
	$\displaystyle\gamma_{n}^{i}$	$\displaystyle\triangleq$	$\displaystyle\beta_{n}(c_{n}^{i})\exp(-\zeta_{n}(c_{n}^{i})(\xi_{n}^{i})^{2}).$		(32)

Similarly, the posterior on $s_{n}$ after the $i^{\mathrm{th}}$ iteration obeys

$\displaystyle\hat{p}^{i+1}(s_{n}\,|\,\boldsymbol{y})$

$\displaystyle\propto$

$\displaystyle\nu^{i}_{f_{n}\to s_{n}}(s_{n})\nu_{h_{n}\to s_{n}}(s_{n}),$

(33)

where

$\displaystyle\nu^{i}_{f_{n}\to s_{n}}(s_{n})$

$\displaystyle\propto$

$\displaystyle\int_{\theta_{n}}f_{n}(\theta_{n},s_{n})\prod_{l=1}^{M}\nu^{i}_{g_{l}\to\theta_{n}}(\theta_{n}).$

(34)

Since $f_{n}(\theta_{n},s_{n})=s_{n}\mathcal{N}(\theta_{n};0,\sigma_{n}^{2})+(1-s_{n})\delta(\theta_{n})$, it can be seen that the corresponding log-likelihood ratio (LLR) is

$\displaystyle L^{i+1}_{n}\triangleq\ln\frac{\nu_{f_{n}\to s_{n}}^{i}(s_{n}\!=\!1)}{\nu_{f_{n}\to s_{n}}^{i}(s_{n}\!=\!0)}=\ln\frac{1-\lambda_{n}^{i}}{\gamma_{n}^{i}\lambda_{n}^{i}}.\quad$

(35)

Clearly, the LLR $L_{n}^{i+1}$ and the likelihood function $\nu_{f_{n}\to s_{n}}^{i}(.)$ express the same information, but in different ways.

The procedure described thus far updates $\mathcal{O}(MN)$ variables per LBP iteration, which is impractical since $MN$ can be very large. In [23], Donoho, Maleki, and Montanari proposed, for the i.i.d case, further approximations that yield a “first-order” approximate message passing (AMP) algorithm that allows the update of only $\mathcal{O}(N)$ variables per LBP iteration, essentially by approximating the differences among the outgoing means/variances of the $g_{m}(\cdot)$ nodes (i.e., $z^{i}_{mn}$ and $c^{i}_{mn}$) as well as the differences among the outgoing means/variances of the $\theta_{n}$ nodes (i.e., $\mu_{n}^{i}$ and $v_{n}^{i}$). These resulting algorithm was then rigorously analyzed by Bayati and Montanari in [15]. We now summarize a straightforward extension of the i.i.d AMP algorithm from [23] to the case of an independent but non-identical Bernoulli-Gaussian prior (13):

	$\displaystyle\xi^{i}_{n}$	$\displaystyle=$	$\displaystyle\textstyle\sum_{m=1}^{M}A_{mn}z^{i}_{m}+\mu^{i}_{n}$		(36)
	$\displaystyle\mu^{i+1}_{n}$	$\displaystyle=$	$\displaystyle F_{n}(\xi^{i}_{n};c^{i})$		(37)
	$\displaystyle v^{i+1}_{n}$	$\displaystyle=$	$\displaystyle G_{n}(\xi^{i}_{n};c^{i})$		(38)
	$\displaystyle z^{i+1}_{m}$	$\displaystyle=$	$\displaystyle\textstyle y_{m}-\sum_{n=1}^{N}A_{mn}\mu^{i}_{n}+\frac{z^{i}_{m}}{M}\sum_{n=1}^{N}F^{\prime}_{n}(\xi^{i}_{n};c^{i})$		(39)
	$\displaystyle c^{i+1}$	$\displaystyle=$	$\displaystyle\textstyle\sigma^{2}+\frac{1}{M}\sum_{n=1}^{N}v^{i+1}_{n},$		(40)

where $F_{n}(.;.)$, $G_{n}(.;.)$ and $F^{\prime}_{n}(.;.)$ are defined as

	$\displaystyle F_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\frac{\alpha_{n}(c)}{1+\tau_{n}(\xi;c)}\xi$		(41)
	$\displaystyle G_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\tau_{n}(\xi;c)F_{n}(\xi;c)^{2}+c\,\xi^{-1}\!F_{n}(\xi;c)$		(42)
	$\displaystyle F^{\prime}_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\frac{\alpha_{n}(c)\big{[}1+\tau_{n}(\xi;c)\left(1+2\zeta(c_{n})\xi^{2}\right)\!\!\big{]}}{[1+\tau_{n}(\xi;c)]^{2}}\!\quad$		(43)
	$\displaystyle\tau_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\beta_{n}(c)\exp(-\zeta_{n}(c)\xi^{2}).$		(44)

For the first turbo iteration (i.e., $t\!=\!1$), we initialize AMP using $z^{i=1}_{m}=y_{m}$, $\mu^{i=1}_{n}=0$, and $c^{i=1}_{n}\gg\sigma_{n}^{2}$ for all $m,n$. For subsequent turbo iterations (i.e., $t\!>\!1$), we initialize AMP by setting $z^{i=1}_{m},\mu^{i=1}_{n},c^{i=1}_{n}$ equal to the final values of $z^{i}_{m},\mu^{i}_{n},c^{i}_{n}$ generated by AMP at the previous turbo iteration. We terminate the AMP iterations as soon as either $\|\boldsymbol{\mu}^{i}-\boldsymbol{\mu}^{i-1}\|_{2}<10^{-5}$ or a maximum of $10$ AMP iterations have elapsed. Similarly, we terminate the turbo iterations as soon as either $\|\boldsymbol{\mu}^{\scriptscriptstyle(t)}-\boldsymbol{\mu}^{\scriptscriptstyle(t-1)}\|_{2}<10^{-5}$ a maximum of $10$ turbo iterations have elapsed. The final value of $\boldsymbol{\mu}=[\mu_{1},\dots,\mu_{N}]^{T}$ is output at the signal estimate $\hat{\boldsymbol{\theta}}$.

We now describe how the precisions $\{\rho_{j}\}$ are learned. First, we recall that $\rho_{j}$ describes the apriori precision on the active coefficients at the $j^{th}$ level, i.e., on $\{\theta_{n}\}_{n\in\mathcal{S}_{j}}$, where the corresponding index set $\mathcal{S}_{j}\triangleq\{n\in\mathcal{W}_{j}:s_{n}=1\}$ is of size $K_{j}\triangleq|\mathcal{S}_{j}|$. Furthermore, we recall that the prior on $\rho_{j}$ was chosen as in (4). Thus, if we had access to the true values $\{\theta_{n}\}_{n\in\mathcal{S}_{j}}$, then (2) implies that

$\displaystyle p(\theta_{n}\,|\,n\in\mathcal{S}_{j})$

$\displaystyle=$

$\displaystyle\mathcal{N}(\theta_{n};0,\rho_{j}),$

(45)

which implies³³3 This posterior results because the chosen prior is conjugate [16] for the likelihood in (45). that the posterior on $\rho_{j}$ would take the form of $\mathrm{Gamma}(\hat{a}_{j},\hat{b}_{j})$ where $\hat{a}_{j}=a_{j}+\frac{1}{2}K_{j}$ and $\hat{b}_{j}=b_{j}+\frac{1}{2}\sum_{n\in\mathcal{S}_{j}}\theta_{n}^{2}$. In practice, we don’t have access to the true values $\{\theta_{n}\}_{n\in\mathcal{S}_{j}}$ nor to the set $\mathcal{S}_{j}$, and thus we propose to build surrogates from the SSR outputs. In particular, to update $\rho_{j}$ after the $t^{\mathrm{th}}$ turbo iteration, we employ

	$\displaystyle\mathcal{S}^{\scriptscriptstyle(t)}_{j}$	$\displaystyle\triangleq$	$\displaystyle\{n\in\mathcal{W}_{j}:L^{\scriptscriptstyle(t)}_{n}>0\}$		(46)
	$\displaystyle K^{\scriptscriptstyle(t)}_{j}$	$\displaystyle\triangleq$	$\displaystyle|\mathcal{S}^{\scriptscriptstyle(t)}_{j}|,$		(47)

and $\{\mu_{n}^{\scriptscriptstyle(t)}\}_{n\in\mathcal{S}^{\scriptscriptstyle(t)}_{j}}$, where $L^{\scriptscriptstyle(t)}_{n}$ and $\mu^{\scriptscriptstyle(t)}_{n}$ denote the final LLR on $s_{n}$ and the final MMSE estimate of $\theta_{n}$, respectively, at the $t^{\mathrm{th}}$ turbo iteration. These choices imply the hyperparameters

	$\displaystyle\hat{a}_{j}^{\scriptscriptstyle(t+1)}$	$\displaystyle=$	$\displaystyle\textstyle a_{j}+\frac{1}{2}K_{j}^{\scriptscriptstyle(t)}$		(48)
	$\displaystyle\hat{b}_{j}^{\scriptscriptstyle(t+1)}$	$\displaystyle=$	$\displaystyle\textstyle b_{j}+\frac{1}{2}\sum_{n\in\mathcal{S}^{\scriptscriptstyle(t)}_{j}}(\mu_{n}^{\scriptscriptstyle(t)})^{2}.$		(49)

Finally, to perform SSR at turbo iteration $t+1$, we set the variances $\{\sigma_{n}^{2}\}_{n\in\mathcal{W}_{j}}$ equal to the inverse of the expected precisions, i.e., $1/\operatorname{E}\{\rho_{j}\}=\hat{b}_{j}^{\scriptscriptstyle(t+1)}/\hat{a}_{j}^{\scriptscriptstyle(t+1)}$. The noise variance $\sigma^{2}$ is learned similarly from the SSR-estimated residual.

Next, we describe how the transition probabilities $\{\pi_{j}^{11}\}$ are learned. First, we recall that $\pi_{j}^{11}$ describes the probability that a child at level $j+1$ is active (i.e., $s_{n}=1$) given that his parent (at level $j$) is active. Furthermore, we recall that the prior on $\pi_{j}^{11}$ was chosen as in (7). Thus if we knew that there were $K_{j}$ active coefficients at level $j$, of which $C_{j}$ had active children, then⁴⁴4 This posterior results because the chosen prior is conjugate to the Bernoulli likelihood [16]. the posterior on $\pi_{j}^{11}$ would take the form of $\mathrm{Beta}(\hat{c}_{j},\hat{d}_{j})$, where $\hat{c}_{j}=c_{j}+C_{j}$ and $\hat{d}_{j}=d_{j}+K_{j}-C_{j}$. In practice, we don’t have access to the true values of $K_{j}$ and $C_{j}$, and thus we build surrogates from the SSR outputs. In particular, to update $\pi_{j}^{11}$ after the $t^{\mathrm{th}}$ turbo iteration, we approximate $s_{n}=1$ by the event $L^{\scriptscriptstyle(t)}_{n}>0$, and based on this approximation set $K^{\scriptscriptstyle(t)}_{j}$ (as in (47)) and $C^{\scriptscriptstyle(t)}_{j}$. The corresponding hyperparameters are then updated as

	$\displaystyle\hat{c}_{j}^{\scriptscriptstyle(t+1)}$	$\displaystyle=$	$\displaystyle c_{j}+C_{j}^{\scriptscriptstyle(t)}$		(50)
	$\displaystyle\hat{d}_{j}^{\scriptscriptstyle(t+1)}$	$\displaystyle=$	$\displaystyle d_{j}+K_{j}^{\scriptscriptstyle(t)}-C_{j}^{\scriptscriptstyle(t)}.$		(51)

Finally, to perform SSR at turbo iteration $t+1$, we set the transition probabilities ${\pi}_{j}^{11}$ equal to the expected value $\hat{c}_{j}^{\scriptscriptstyle(t+1)}/(\hat{c}_{j}^{\scriptscriptstyle(t+1)}+\hat{d}_{j}^{\scriptscriptstyle(t+1)})$. The parameters $\pi_{0}^{1}$, $\pi_{-1}^{1}$, and $\{\pi_{j}^{00}\}$ are learned similarly.

Until now, we have focused on the Bernoulli-Gaussian (BG) signal model (2). In this section, we describe the modifications needed to handle the Gaussian mixture (GM) model

$\displaystyle p(\theta_{n}\,|\,s_{n})$

$\displaystyle=$

$\displaystyle s_{n}\mathcal{N}(\theta_{n};0,\sigma_{n,\text{\sf L}}^{2})+(1-s_{n})\mathcal{N}(\theta_{n};0,\sigma_{n,\text{\sf S}}^{2}),$

where $\sigma_{n,\text{\sf L}}^{2}$ denotes the variance of “large” coefficients and $\sigma_{n,\text{\sf S}}^{2}$ denotes the variance of “small” ones. For either the BG or GM prior, AMP is performed using the steps (36)-(40). For the BG case, the functions $F_{n}(.;.)$, $G_{n}(.;.)$, $F^{\prime}_{n}(.;.)$, and $\tau_{n}(.;.)$ are given in (41)–(44), whereas for the GM case, they take the form

	$\displaystyle F_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\frac{\bar{\alpha}_{n,\text{\sf L}}(c)+\bar{\alpha}_{n,\text{\sf S}}(c)\bar{\tau}_{n}(\xi,c)}{1+\bar{\tau}_{n}(\xi,c)}\xi$		(53)
	$\displaystyle G_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\frac{\bar{\tau}_{n}(\xi,c)\xi^{2}[\bar{\alpha}_{n,\text{\sf L}}(c)-\bar{\alpha}_{n,\text{\sf S}}(c)]^{2}}{[1+\bar{\tau}_{n}(\xi,c)]^{2}}+\,c\,\xi^{-1}\!F_{n}(\xi;c)$
	$\displaystyle F^{\prime}_{n}(\xi;c)$	$\displaystyle=$	$\displaystyle\frac{\bar{\tau}_{n}(\xi;c)\bar{\alpha}_{n,\text{\sf L}}(c)(1+\bar{\tau}_{n}(\xi;c)-2\xi^{2}\bar{\zeta}_{n}(c))}{[1+\bar{\tau}_{n}(\xi;c)]^{2}}$		(55)
			$\displaystyle+\frac{\bar{\alpha}_{n,\text{\sf S}}(c)(1+\bar{\tau}_{n}(\xi;c)+2\xi^{2}\bar{\zeta}_{n}(c)\bar{\tau}_{n}(\xi;c))}{[1+\bar{\tau}_{n}(\xi;c)]^{2}}\quad$
	$\displaystyle\bar{\tau}_{n}(\xi,c)$	$\displaystyle=$	$\displaystyle\bar{\beta}_{n}(c)\exp(-\bar{\zeta}_{n}(c)\xi^{2}),$		(56)

where

	$\displaystyle\bar{\alpha}_{n,\text{\sf L}}(c)$	$\displaystyle\triangleq$	$\displaystyle\frac{\sigma_{n,\text{\sf L}}^{2}}{c+\sigma_{n,\text{\sf L}}^{2}}$		(57)
	$\displaystyle\bar{\alpha}_{n,\text{\sf S}}(c)$	$\displaystyle\triangleq$	$\displaystyle\frac{\sigma_{n,\text{\sf S}}^{2}}{c+\sigma_{n,\text{\sf S}}^{2}}$		(58)
	$\displaystyle\bar{\beta}_{n}(c)$	$\displaystyle\triangleq$	$\displaystyle\frac{1-\lambda_{n}}{\lambda_{n}}\sqrt{\frac{c+\sigma_{n,\text{\sf L}}^{2}}{c+\sigma_{n,\text{\sf S}}^{2}}}$		(59)
	$\displaystyle\bar{\zeta}_{n}(c)$	$\displaystyle\triangleq$	$\displaystyle\frac{\sigma_{n,\text{\sf L}}^{2}-\sigma_{n,\text{\sf S}}^{2}}{2(c+\sigma_{n,\text{\sf L}}^{2})(c+\sigma_{n,\text{\sf S}}^{2})}.$		(60)

Likewise, for the BG case, the extrinsic LLR is given by (35), whereas for the GM case, it becomes

$\displaystyle L^{i+1}_{n}\triangleq\ln\frac{\nu_{f_{n}\to s_{n}}^{i}(s_{n}\!=\!1)}{\nu_{f_{n}\to s_{n}}^{i}(s_{n}\!=\!0)}=\ln\frac{1-\lambda_{n}^{i}}{\bar{\tau}_{n}(\xi_{n}^{i};c_{n}^{i})\lambda_{n}^{i}}.\quad$

(61)

The proposed turbo⁵⁵5An implementation of our algorithm can be downloaded from http://www.ece.osu.edu/~schniter/turboAMPimaging approach to compressive imaging was compared to several other tree-sparse reconstruction algorithms: ModelCS [3], HMT+IRWL1 [4], MCMC [5], variational Bayes (VB) [6]; and to several simple-sparse reconstruction algorithms: CoSaMP [7], SPGL1 [25], and Bernoulli-Gaussian (BG) AMP. All numerical experiments were performed on $128\!\times\!128$ (i.e., $N\!=\!16384$) grayscale images using a $J=4$-level 2D Haar wavelet decomposition, yielding $8^{2}\!=\!64$ approximation coefficients and $3\!\times\!8^{2}\!=\!192$ individual Markov trees. In all cases, the measurement matrix $\boldsymbol{\Phi}$ had i.i.d Gaussian entries. Unless otherwise specified, $M\!=\!5000$ noiseless measurements were used. We used normalized mean squared error (NMSE) $\|\boldsymbol{x}-\boldsymbol{\hat{x}}\|_{2}^{2}/\|\boldsymbol{x}\|_{2}^{2}$ as the performance metric.

We now describe how the hyperparameters were chosen for the proposed Turbo schemes. Below, we use $N_{j}$ to denote the total number of wavelet coefficients at level $j$, and $N_{-1}$ to denote the total number of approximation coefficients. For both Turbo-BG and Turbo-GM, the Beta hyperparameters were chosen so that $c\!+\!d\!=\!N_{0}$, $\underline{c}\!+\!\underline{d}\!=\!N_{-1}$ and $c_{j}\!+\!d_{j}\!=\!N_{j}~\forall j$ with $\operatorname{E}\{p_{0}^{1}\}\!=\!1/N$, $\operatorname{E}\{p_{-1}^{1}\}\!=\!1-10^{-6}$, $\operatorname{E}\{p_{j}^{00}\}\!=\!1/N~\forall j$, and $\operatorname{E}\{p_{j}^{11}\}\!=\!0.5~\forall j$. These informative hyperparameters are similar to the “universal” recommendations in [26] and, in fact, identical to the ones suggested in the MCMC work [5]. For Turbo-BG, the hyperparameters for the signal precisions $\{\sigma^{-2}_{n}\}_{n\in\mathcal{W}_{j}}$ were set to $a_{j}\!=\!1~\forall j$ and $[b_{-1},\dots,b_{3}]\!=\![10,1,1,0.1,0.1]$. This choice is motivated by the fact that wavelet coefficient magnitudes are known to decay exponentially with scale $j$ (e.g., [26]). Meanwhile, the hyperparameters for the noise precision $\sigma^{-2}$ were set to $(a,b)=(1,10^{-6})$. Although the measurements were noiseless, we allow Turbo-BG a nonzero noise variance in order to make up for the fact that the wavelet coefficients are not exactly sparse, as assumed by the BG signal model. (We note that the same was done in the BG-based work [6, 5].) For Turbo-GM, the hyperparameters $(a_{j,\text{\sf L}},b_{j,\text{\sf L}})$ for the signal precisions $\{\sigma^{-2}_{n,\text{\sf L}}\}_{n\in\mathcal{W}_{j}}$ were set at the values of $(a_{j},b_{j})$ for the BG case, while the hyperparameters $(a_{j,\text{\sf S}},b_{j,\text{\sf S}})$ for $\{\sigma^{-2}_{n,\text{\sf S}}\}_{n\in\mathcal{W}_{j}}$ were set as $a_{j,\text{\sf S}}\!=\!1$ and $b_{j,\text{\sf S}}=10^{-6}b_{j,\text{\sf L}}$. Meanwhile, the noise variance $\sigma^{2}$ was assumed to be exactly zero, because the GM signal prior is capable of modeling non-sparse wavelet coefficients.