that SE allows to compute the highest equilibrium of Gibbs free energy [@DMM2011; @Bayati2011; @JavanmardMontanari2012; @Donoho2013; @Bayati2015], which corresponds to the local optimum $D_2$ in Section \[sec:phaseTrans\]. Hence, AMP often achieves the same MSE as BP and we use AMP simulation results to demonstrate that the MMSE can often be achieved.[^11] Considering the structure of $\F$, we simplify the AMP algorithm in Barbier and Krzakala [@Barbier2015] to obtain Algorithm \[algo:AMP\_MMV\],[^12] where $\{\Sigma_j\}_{j=1}^J$, $\{R^j_l\}_{j=1}^J$, $\{a_l^j\}_{j=1}^J$ and $\{v_l^j\}_{j=1}^J$ refer to sets of all intermediate variables $\Sigma_j$, pseudodata $R^j_l$, estimates $a_l^j$, and variances $v^j_l,\ j\in\{1,...,J\},\ l\in\{1,...,N\}$, respectively. The current iteration $t$, change in the estimate $\delta$, and intermediate variables $\Theta_j,\ j\in\{1,...,J\}$, are scalars. The intermediate variables $\q^j$ and $\w^j$ are vectors of length $M$. The functions $f_{a_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)$ and $f_{v_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)$ are given by $$\begin{split}
&f_{a_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)=\\
&\frac{\rho\frac{1}{\Sigma_j+1}\{R^j_l\}_{j=1}^J}{\rho+(1-\rho){\prod_{j=1}^{J}}\l\{\sqrt{1+\frac{1}{\Sigma_j}}\exp\l[-\frac{(R^j_l)^2}{2\Sigma_j(\Sigma_j+1)}\r]\r\}},
\end{split}$$ $$\begin{split}
&f_{v_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)=-\l[f_{a_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)\r]^2\\
&+\frac{\rho\frac{1}{\Sigma_j+1}\l[(\{R^j_l\}_{j=1}^J)^2\frac{1}{\Sigma_j+1}+\Sigma_j\r]}{\rho+(1-\rho){\prod_{j=1}^{J}}\l\{\sqrt{1+\frac{1}{\Sigma_j}}\exp\l[-\frac{(R^j_l)^2}{2\Sigma_j(\Sigma_j+1)}\r]\r\}},
\end{split}$$ for $J$-dimensional Bernoulli-Gaussian signals .

\
[**Inputs:**]{} Maximum number of iterations $t_{max}$, threshold $\epsilon$, sparsity rate $\rho$, noise variance $\Delta$, measurements $\y^j$, and measurement matrices $\F^j, \forall j$\
[**Initialize:**]{} $t=1,\delta=\infty,\w^j=\y^j,\Theta_j=0,v^j_l=\rho\Delta,a^j_l=0,\forall l,j$\
$\q^j=\frac{\y^j-\w^j}{\Delta+\Theta_j}$\
$\Theta_j=\frac{1}{N}\sum_{l=1}^N v^j_l$\
$\w^j=\F^j \a^j-\Theta_j \q^j$\
$\Sigma_j=\frac{N(\Delta+\Theta_j)}{M}$ // Scalar channel noise variance\
$\R^j=\a^j+\Sigma_j (\F^j)^T \frac{\y^j-\w^j}{\Delta+\Theta_j}$ // Pseudodata\
$\widehat{\a}^j=\a^j$ // Save current estimate\
$\{v^j_l\}_{j=1}^J=f_{v_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)$ // Variance\
$\{a^j_l\}_{j=1}^J=f_{a_l}(\{\Sigma_j\}_{j=1}^J,\{R^j_l\}_{j=1}^J)$ // Estimate\
  $t=t+1$ // Increment iteration 