\h^{(j)}$ // Pseudo data\
$\widehat{\a}^{(j)}=\widehat{\x}^{(j)}$ // Save current estimate \[line:endForLoop\]\
 $\Delta_v=\sum_{j=1}^J \Delta_v^{(j)}$\[line:mean\_delta\]\
$\widehat{\x}_n=f_{a_n}(\Delta_v,\q_n)$ // Estimate\[line:mean\]\
$\s_n=[s_n^{(1)},\ldots,s_n^{(J)}]=f_{v_n}(\Delta_v,\q_n)$ // Variance\[line:var\]\
  $t=t+1$ // Increment iteration index\
  $\delta=\frac{1}{NJ}\sum_{n=1}^N\sum_{j=1}^J\l(\widehat{x}^{(j)}_n-\widehat{a}^{(j)}_n\r)^2$ // Change\
[**Outputs:**]{} Estimate $\widehat{\x}^{(j)}$, pseudo data $\q^{(j)},\ \forall j$, and scalar channel noise variance $\Delta_v$

[**MP algorithm:**]{} For the first part, we modify the generalized approximate message passing (GAMP) algorithm [@RanganGAMP2011ISIT], which is an implementation of MP, and list the pseudo code in Algorithm \[algo:AMP\_MMV\]. The notation diag$(\boldsymbol{\Theta}^{(j)})$ denotes a diagonal matrix whose entries along the diagonal are $\boldsymbol{\Theta^{(j)}}$, and the power-of-two in Lines \[line:Theta\] and \[line:scalarNoiseVar\] is applied element-wise. The function $g_{\text{out}}$ in Lines \[line:g\_out\] and \[line:deriv\_g\_out\] is given by $$\label{eq:g_out}
g_{\text{out}}\l(k,y,\Theta\r)=\frac{1}{\Theta}(\mathbb{E}[w|k,y,\Theta]-k),$$ where we omit the subscripts and super-scripts for brevity, and the expectation is taken over the pdf, $$\label{eq:prob_w}
f(w|k,y,\Theta)\propto f(y|w)\text{exp}\l[-\frac{(w-k)^2}{2\Theta}\r].$$

For the special case of AWGN channels, $$\label{eq:AWGN}
y=w+z,$$ where $z\sim \mathcal{N}(0,\Delta_z)$, we obtain $g_{\text{out}}(k,y,\Theta)=\frac{y-k}{\Delta_z+\Theta}$ [@RanganGAMP2011ISIT]. In Appendix \[app:logit\], we also briefly present the derivation for an i.i.d. parallel logistic channel, $$\label{eq:logit}
f(y|w)=\delta(y-1)\frac{1}{1+\text{exp}(-aw)}+\delta(y)\frac{\text{exp}(-aw)}{1+\text{exp}(-aw)},$$ where $a$ is a scaling factor.[^5]

For the special case of i.i.d. [*joint Bernoulli-Gaussian signals*]{} where $\phi(\x_n)\sim \mathcal{N}(0,\mathbb{I})$ in  and $\mathbb{I}$ is an identity matrix, $f_{a_n}$ and $f_{v_n}$ in Lines \[line:mean\] and \[line:var\] are given below, $$\label{eq:denoiser}
f_{a_n}(\Delta_v,\q_n)=\frac{\rho}{C(\Delta_v+1)}\q_n,
$$ 