---
abstract: |
    Compressed sensing (CS) demonstrates that sparse signals can be estimated from under-determined linear systems. Distributed CS (DCS) further reduces the number of measurements by considering joint sparsity within signal ensembles. DCS with jointly sparse signals has applications in multi-sensor acoustic sensing, magnetic resonance imaging with multiple coils, remote sensing, and array signal processing. Multi-measurement vector (MMV) problems consider the estimation of jointly sparse signals under the DCS framework. Two related MMV settings are studied. In the first setting, each signal vector is measured by a different independent and identically distributed (i.i.d.) measurement matrix, while in the second setting, all signal vectors are measured by the same i.i.d. matrix. Replica analysis is performed for these two MMV settings, and the minimum mean squared error (MMSE), which turns out to be identical for both settings, is obtained as a function of the noise variance and number of measurements. To showcase the application of MMV models, the MMSE’s of complex CS problems with both real and complex measurement matrices are also analyzed. Multiple performance regions for MMV are identified where the MMSE behaves differently as a function of the noise variance and the number of measurements.

    Belief propagation (BP) is a CS signal estimation framework that often achieves the MMSE asymptotically. A phase transition for BP is identified. This phase transition, verified by numerical results, separates the regions where BP achieves the MMSE and where it is suboptimal. Numerical results also illustrate that more signal vectors in the jointly sparse signal ensemble lead to a better phase transition.
author:
- 'Junan Zhu, , Dror Baron, , and Florent Krzakala [^1] [^2] [^3]'
title: |
    Performance Limits\
    for Noisy Multi-Measurement Vector Problems
---

[*Keywords*]{}:Approximate message passing, multi-measurement vector problem, replica analysis.

Introduction
============

Compressed sensing (CS) [@CandesRUP; @DonohoCS; @BaraniukCS2007] demonstrates that sparse signals can be estimated from under-determined linear measurements. Owing to the potential for radically reduced measurement rates, CS has become an active research area within signal processing. CS has many application areas including magnetic resonance imaging [@JuYeKi07; @JuSuNaKiYe09], communication [@Cotter2002scemp], and remote sensing [@Ma2009deblur].

Distributed CS (DCS) [@HN05; @DuarteWakinBaronSarvothamBaraniuk2013] is based on the premise that joint sparsity within signal ensembles enables a further reduction in the number of measurements. Motivated by sensor networks [@Pottie2000], preliminary work in DCS [@Duarte2006IPSN; @HN05; @DuarteWakinBaronSarvothamBaraniuk2013] showed that the number of measurements required per sensor must account for the minimum features unique to that sensor while features that are common to multiple sensors are amortized. DCS led to a proliferation of research on the multi-measurement vector (MMV) problem [@chen2006trs; @cotter2005ssl; @Mishali08rembo; @Berg09jrmm; @LeeKimBreslerYe2011; @LeeBreslerJunge2012; @YeKimBresler2015]. The MMV problem considers the estimation of a set of sparse signal vectors that share common supports, and has applications such as radar array signal processing, acoustic sensing with multiple speakers, magnetic resonance imaging with multiple coils [@JuYeKi07; @JuSuNaKiYe09], and diffuse optical tomography using multiple illumination patterns. In MMV, thanks to the common support, the number of sparse coefficients that can be successfully estimated increases with the number of measurements. This property was evaluated rigorously for noiseless measurements using $l_0$ minimization [@DuarteWakinBaronSarvothamBaraniuk2013]. To address measurement noise, estimation approaches for MMV problems have included greedy algorithms such as SOMP [@tropp2006ass; @chen2006trs], $l_1$ convex relaxation [@malioutov2005ssr; @tropp2006ass2], and M-FOCUSS [@cotter2005ssl]. REduce MMV and BOost (ReMBo) has been shown to outperform conventional methods [@Mishali08rembo], and subspace methods have also been used to solve MMV problems [@LeeBreslerJunge2012; @YeKimBresler2015]. Statistical approaches [@ZinielSchniter2011] often achieve the oracle minimum mean squared error (MMSE). However, the performance limits of MMV signal estimation in the presence of measurement noise have not been studied.

Replica analysis is a statistical physics method that can be used to analyze the MMSE and phase transition for inverse problems [@Tanaka2002; @GuoVerdu2005; @Montanari2006; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015; @Lesieur2015]. Barbier and Krzakala [@Barbier2015] studied the MMSE for estimating superposition codes using replica analysis. In this paper, we extend the derivation in Barbier and Krzakala [@Barbier2015] to two related yet different MMV settings: ([*i*]{}) $J$ jointly sparse signals are measured by $J$ different dense matrices that are independent and identically distributed (i.i.d.), and ([*ii*]{}) $J$ jointly sparse signals are measured by $J$ identical i.i.d. matrices. We only consider dense i.i.d. Gaussian matrices in this work, while our analysis can be extended to other i.i.d. matrices easily.

We make several contributions in this paper. First, we obtain the information theoretic MMSE for the two MMV settings above under the Bayesian setting. Second, we show that in the large system limit the MMSE’s for these two settings are identical to the single measurement vector (SMV) problem with a dense measurement matrix and a block sparse signal with fixed length blocks. Third, we derive the MMSE for SMV complex CS problems by noticing that SMV complex CS is essentially an MMV problem. Fourth, we identify several performance regions for MMV, where the MMSE has different characteristics based on channel noise variance and measurement rate. Finally, we find a phase transition for belief propagation algorithms (BP) [@DMM2009; @CSBP2010; @Bayati2011; @Montanari2012; @Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015] applied to MMV problems, which separates regions where BP achieves the MMSE asymptotically and where it is suboptimal. BP simulation results confirm the phase transition results.

The remainder of the paper is organized as follows. We introduce our signal and measurement models in Section \[sec:model\], followed by replica analyses for two MMV settings as well as two SMV complex CS problems in Section \[sec:main\]. Section \[sec:proof\] proves the results of Section \[sec:main\]. Numerical results are discussed in Section \[sec:numeric\] and we conclude in Section \[sec:conclusion\].

[**Notations:**]{} In this paper, bold capital letters represent matrices, bold lower case letters represent vectors, and normal font lower case letters represent scalars. The entry (scalar) in the $\mu$-th row, $l$-th column of a matrix $\F$ is denoted by $F_{\mu,l}$, where the comma is often omitted. The $\mu$-th entry (scalar) in a vector $\z$ is denoted by $z_{\mu}$.

Signal and Measurement Models {#sec:model}
=============================

[**Signal model**]{}: We consider an ensemble of $J$ signal vectors, $\underline{\s}^j\in\mathbb{R}^N,\ j\in\{1,...,J\}$, where $j$ is the index of the signal. Consider a [*super symbol*]{} $\s_l=[\underline{s}_l^1,...,\underline{s}_l^J]^T,\ l\in\{1,...,N\}$, where $[\cdot]^T$ denotes the transpose. The super symbol $\s_l$ follows a $J$-dimensional Bernoulli-Gaussian distribution, $$\label{eq:jsm}
P(\s_l)=\rho \phi(\s_l)+(1-\rho)\delta(\s_l),$$ where $\rho$ is the sparsity rate, $\phi(\s_l)$ is a $J$-dimensional Gaussian distribution with zero mean and identity covariance matrix, and $\delta(\s_l)$ is the delta function for $J$-dimensional vectors.

\[def:jointly\_sparse\] [*Ensembles of signals that obey  are called jointly sparse.*]{}

[**Measurement models**]{}: Each signal $\underline{\s}^j$ is measured by an i.i.d. Gaussian measurement matrix $\underline{\F}^j\in\mathbb{R}^{M\times N}$, $\underline{F}_{\mu l}^j \sim \mathcal{N}(0,1/N)$, where $\mu$ refers to the row index and $l$ is the column index. The measurements $\underline{\y}^j$ are corrupted by i.i.d. Gaussian noise $\underline{\z}^j$ consisting of entries $\underline{z}_{\mu}^j\sim \mathcal{N}(0,\Delta)$, $$\label{eq:MMVmodel}
\underline{\y}^j=\underline{\F}^j\underline{\s}^j+\underline{\z}^j,\quad j\in\{1,\cdots,J\}.$$ When the number of signal vectors $J=1$, this MMV model  becomes an SMV problem. Our analyses in this paper are readily extended to other i.i.d. matrices, jointly sparse signals , and other i.i.d. noise distributions.

\[def:MMV\_set1\] [*The setting MMV-1 refers to the measurement model in  with all matrices $\underline{F}^j$ being different.*]{}

\[def:MMV\_set2\] [*The setting MMV-2 refers to the measurement model in  with all matrices $\underline{F}^j$ being equal.*]{}

In the signal model  and measurement model , the sparsity rate $\rho$, channel noise variance $\Delta$, and number of channels $J$ are constant.

\[def:largeSystemLimit\] The signal length $N$ scales to infinity, and the number of measurements $M=M(N)$ depends on $N$ and also scales to infinity, where the ratio approaches a positive constant $R$ for practical problems, $$\label{eq:measurementRate}
\lim_{N\rightarrow\infty} \frac{M(N)}{N} = R>0.$$

We call $R$ the measurement rate.

Replica Analyses for MMV Settings {#sec:main}
=================================

Section \[sec:model\] discussed two MMV settings. Both settings have applications in real-world problems such as magnetic resonance imaging [@JuYeKi07; @JuSuNaKiYe09] and sensor networks [@Pottie2000]. Although numerous algorithms for MMV signal estimation have been proposed [@tropp2006ass; @chen2006trs; @malioutov2005ssr; @tropp2006ass2; @cotter2005ssl; @Mishali08rembo; @ZinielSchniter2011], what is missing is an information theoretic analysis of the best possible mean squared error (MSE) performance. Throughout this paper, we only consider the MSE as our performance metric.

Statistical physics background and replica method {#sec:set1}
-------------------------------------------------

![Illustration of MMV channel  with $J=3$ signal vectors (left), and one of its possible SMV forms (right). Different background patterns differentiate entries from different channels, and blank space denotes zeros.[]{data-label="fig:channel"}](MMV_channel_v2.png){width="8.5cm"}

In order to express  using a single channel, we transform it to an SMV form. One possible way to do so is illustrated in Fig. \[fig:channel\]. The equivalent SMV problem is $$\label{eq:MMVchannel}
  \y=\F\s+\z,$$ where $\F\in\mathbb{R}^{MJ\times NJ}$ is the matrix, $\y\in\mathbb{R}^{MJ}$ are the measurements, and the noise is $\z\in\mathbb{R}^{MJ}$. Entries of the signal vectors $\underline{\s}^j$, measurement vectors $\underline{\y}^j$, and noise vectors $\underline{\z}^j$  form the SMV signal $\s$, measurements $\y$, and noise $\z$  with $$s_{(l-1)J+j}=\underline{s}^{j}_l,\ y_{(j-1)M+\mu}=\underline{y}^j_{\mu},\ \text{and}\ z_{(j-1)M+\mu}=\underline{z}^j_{\mu},$$ respectively. Entries of the matrix $\underline{\F}^j$  form the SMV matrix $\F$  with $F_{(j-1)M+\mu,(l-1)J+j}=\underline{F}^j_{\mu l}$; other entries of $\F$ are zeros. The posterior for the estimate $\x\in\mathbb{R}^{NJ}$, comprised of super symbols $\x_l=[x_{(l-1)J+1},...,x_{lJ}]^T,\ l\in\{1,...,N\}$, is $$\label{eq:pxy2}
  P(\x|\y)=\frac{1}{Z}{\prod_{l=1}^{N}}P(\x_l){\prod_{\mu=1}^{MJ}}\l[\frac{\operatorname{e}^{-\frac{1}{2\Delta}(\\y_{\mu}-{\sum_{l=1}^{N}}\F_{\mu l}\x_l)^2}}{\sqrt{2\pi\Delta}}\r],$$ where $\F_{\mu l}=[F_{\mu,(l-1)J+1}, \ldots, F_{\mu,lJ}]$ is a super symbol highlighted by the dashed area in Fig. \[fig:channel\], and the denominator $Z$ is the partition function [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015], $$\label{eq:partition}
  Z=\int {\prod_{l=1}^{N}}P(\x_l){\prod_{\mu=1}^{MJ}}\l[\frac{\operatorname{e}^{-\frac{1}{2\Delta}(y_{\mu}-{\sum_{l=1}^{N}}\F_{\mu l}\x_l)^2}}{\sqrt{2\pi\Delta}}\r]{\prod_{l=1}^{N}}d\x_l.$$ Note that multi-dimensional integrations such as  are denoted by a single $\int$ operator for brevity. Confining our attention to the Bayesian setting [@Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015], $P(\x_l)$ follows the true distribution , $P(\x_l)=\rho \phi(\x_l)+(1-\rho)\delta(\x_l)$.

By creating an analogy between the channel  and a many-body thermodynamic system [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015], the posterior  can be interpreted as the Boltzmann measure on a disordered system with the following Hamiltonian, $$\label{eq:Hamiltonian}
H(\x)=\sum_{l=1}^N \log [P(\x_l)]+\sum_{\mu=1}^{MJ} \frac{1}{2\Delta}\l(y_{\mu}-\sum_{l=1}^N\F_{\mu l} \x_l\r)^2.$$

The averaged free energy of the disordered system given by  characterizes the thermodynamic properties of the system. Evaluating the fixed points (local maxima) in the free energy expression provides the MMSE for the channel  [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015]. [*Under the assumption of self-averaging*]{} [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015], the free energy is defined as[^4] $$\label{eq:free_energy}
  \mathcal{F}=\lim_{N\rightarrow\infty}\frac{1}{N}\mathbb{E}_{\F,\s,\z}[\log (Z)],$$ which is difficult to evaluate. The replica method [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015] introduces $n$ replicas of the estimate $\x$ as $\x^a,\ a\in\{1,...,n\}$, and the free energy  can be approximated by the replica trick [@Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015], $$\label{eq:replicaTrick}
  \mathcal{F}=\lim_{N\rightarrow\infty}\lim_{n\rightarrow 0} \frac{\mathbb{E}_{\F,\s,\z}[ Z^n]-1}{Nn}.$$ Note that the self-averaging property that leads to  and the replica trick , as well as the replica symmetry assumptions that appear in latter parts of this paper, are assumed to be valid in this work, and their rigorous justification is still an open problem in mathematical physics [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015].[^5]

**Evaluating the free energy**: To evaluate the free energy , we calculate $\mathbb{E}_{\F,\s,\z}\l[Z^n\r]$, where $\cdot_{\F,\s,\z}$ denotes expectation with respect to (w.r.t.) $\F,\s$, and $\z$, and $Z$ is given in : $$\label{eq:EZn1}
  \mathbb{E}_{\F,\s,\z}\l[Z^n\r]\!=\!\frac{\mathbb{E}_{\s}\!\l[\displaystyle{\int {\prod_{l=1}^{N}}{\prod_{a=1}^{n}}\! P(\x_l^a)\!{\prod_{\mu=1}^{M}}\!\mathbb{X}_{\mu}\!{\prod_{l=1}^{N}}{\prod_{a=1}^{n}}d\x_l^a}\r]}{(2\pi\Delta)^{\frac{nMJ}{2}}},$$ where $$\label{eq:Xmu}
  \mathbb{X}_{\mu}=\mathbb{E}_{\F,\z}\l[\operatorname{e}^{-\frac{1}{2\Delta}{\sum_{j=1}^{J}}{\sum_{a=1}^{n}}(v_{\mu j}^a)^2}\r],$$ $a$ is the replica index, $\x^a_l$ is the $l$-th super symbol of $\x^a$, and $$\label{eq:v_mu_a}
v_{\mu j}^a={\sum_{l=1}^{N}}\F_{\mu+M(j-1),l}(\s_l-\x_l^a)+z_{\mu+M(j-1)}.$$

\[lemma:covIsSame\] In the large system limit, the quantity $\mathbb{X}_{\mu}$  is the same for both MMV-1 and MMV-2.

Lemma \[lemma:covIsSame\] is proved in Section \[sec:proof\]. Because of Lemma \[lemma:covIsSame\], the free energy expressions for MMV-1 and MMV-2 should be identical in the large system limit. We state the result as a theorem and the detailed derivations appear in the Appendix.

\[th:free\_energy\] For settings MMV-1 and MMV-2, the free energy expressions as functions of $E$ are identical in the large system limit and are given in .[^6]

$$\begin{aligned}
  \mathcal{F}(E)&=&-\frac{J}{2}R\l\{\log[2\pi(\Delta+E)]+\frac{\rho+\Delta}{E+\Delta}\r\}\!+\!\int\!P(\s_1)\! \int\! \log\! \l[ \int P(\x_1)\operatorname{e}^{-\frac{\widehat{Q}+\widehat{q}}{2}\x_1^T\x_1+\widehat{m}\x_1^T\s_1+
  \sqrt{\widehat{q}}\h^T\x_1}\!d\x_1\!\r]\!\mathcal{D}\h \ d\s_1\label{eq:free_energy3}\\
  &=&-\frac{J}{2}R\l\{\log[2\pi(\Delta+E)]+\frac{\Delta}{E+\Delta}\r\}+\frac{JR(1-\rho)}{2(R+E+\Delta)}+\rho\int \log \Bigg[ \rho \l(\frac{E+\Delta}{R+E+\Delta}\r)^{J/2}+\nonumber\\
  &\ &(1-\rho)\operatorname{e}^{-\frac{R}{2(E+\Delta)}\g^T\g}\Bigg]\mathcal{D}\g+(1-\rho)\int \log \l[ \rho \l(\frac{E+\Delta}{R+E+\Delta}\r)^{J/2}+(1-\rho)\operatorname{e}^{-\frac{R}{2(R+E+\Delta)}\h^T\h}\r]\mathcal{D}\h.\label{eq:free_energy4}\end{aligned}$$

**MMSE**: The $E$ that maximizes the free energy  the MMSE [@Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015]. After finding the $E_0$ that maximizes the free energy , we obtain the MMSE, $D_0=E_0$, in the large system limit.

The MMSE for MMV-1 and MMV-2 is the same for the same measurement rate $R$, noise variance $\Delta$, and number of signal vectors $J$.

[**Remark 1:**]{} As the reader can see from the proof of Lemma \[th:free\_energy\] in Section \[sec:proof\], the key reason that both MMV-1 and MMV-2 have an identical MMSE is that the entries in the super symbols $\s_l$ and $\x_l^{\cdot}$ are i.i.d. That said, we suspect that the MMSE for MMV-1 and MMV-2 could differ by some higher order terms. If the entries of these super symbols are not i.i.d., which is true in some practical MMV applications [@ZinielSchniter2013MMV], then it becomes more difficult to analyze the covariance matrix $\G_{\mu}$ as in Section \[sec:proof\]. Therefore, we do not have an analysis for non-i.i.d. entries within $\s_l$ and $\x_l^{\{\cdot\}}$. However, we speculate that MMV-1 might have lower MMSE than MMV-2 in that case.

[**Link to SMV with block sparse signal:**]{} The signal $\s$ in  is a block sparse signal comprised of $N$ blocks of length $J$. We study a single measurement vector (SMV) problem by replacing the measurement matrix $\F$ in  with an i.i.d. Gaussian matrix ${\mathbf A}\in\mathbb{R}^{MJ\times NJ}$, i.e., $\y={\mathbf A}\s+\z$. The entries of ${\mathbf A}$ follow the distribution, $A_{\mu l}\sim \mathcal{N}(0,\frac{1}{NJ})$. This SMV is similar to the setting in Barbier and Krzakala [@Barbier2015], except for the different priors and different $\ell_2$ norms in each row of ${\mathbf A}$. We consider these differences while following their derivation [@Barbier2015], and obtain the same free energy expression as . We have also shown that MMV-1 and MMV-2 have the same MMSE in the large system limit. Hence, the three settings have the same free energy expression and their MMSE’s are the same under the same noise variance $\Delta$ and measurement rate $R$ in the large system limit.

Extension to complex SMV {#sec:complex}
------------------------

MMV with jointly sparse signals is a versatile model that can be adapted to other problems. As an example, we show how the MMV model can be used to analyze the MMSE of a complex SMV.[^7] Consider the complex CS channel, $\y^{\mathcal{C}}=\F^{\mathcal{C}}\s^{\mathcal{C}}+\z^{\mathcal{C}}$, where $\s^{\mathcal{C}}=\s^{\mathcal{R}}+i\s^{\mathcal{I}}\in\mathbb{C}^N$, $\F^{\mathcal{C}}=\F^{\mathcal{R}}+i\F^{\mathcal{I}}\in\mathbb{C}^{M\times N}$, $\z^{\mathcal{C}}=\z^{\mathcal{R}}+i\z^{\mathcal{I}}\in\mathbb{C}^M$, $\y^{\mathcal{C}}=\y^{\mathcal{R}}+i\y^{\mathcal{I}}\in\mathbb{C}^M$, $i=\sqrt{-1}$, and $\mathcal{R}$ and $\mathcal{I}$ refer to the real and imaginary parts, respectively. The real and imaginary parts of the entries of $\z^{\mathcal{C}}$ both follow a Gaussian distribution, $z_l^{\mathcal{R}}, z_l^{\mathcal{I}}\sim\mathcal{N}(0,\Delta), l\in\{1,...,M\}$. Assume that the complex signal $\s^{\mathcal{C}}$ is comprised of two jointly sparse signals, $\s^{\mathcal{R}}$ and $\s^{\mathcal{I}}$, that satisfy the $J=2$ dimensional Bernoulli-Gaussian distribution . We can extend the analysis of Section \[sec:set1\] to two settings of complex CS: ([*i*]{}) the measurement matrix $\F^{\mathcal{C}}$ is real, and ([*ii*]{}) $\F^{\mathcal{C}}$ is complex.[^8]

[**Real measurement matrix:**]{} Suppose that $\F^{\mathcal{C}}$ is real, $\F^{\mathcal{C}}=\F^{\mathcal{R}}\in \mathbb{R}^{M\times N}$, and the entries of $\F^{\mathcal{R}}$ follow a Gaussian distribution, $F^{\mathcal{R}}_{\mu l}\sim \mathcal{N}(0,\frac{1}{N})$. Complex CS with a real measurement matrix can be written as real-valued MMV, $$\label{eq:complexRealMat}
    \y^{\mathcal{R}}=\F^{\mathcal{R}}\s^{\mathcal{R}}+\z^{\mathcal{R}}\ \text{and}\     \y^{\mathcal{I}}=\F^{\mathcal{R}}\s^{\mathcal{I}}+\z^{\mathcal{I}},$$ where $\s^R$ and $\s^I$ are jointly sparse and follow . This formulation  fits into MMV-2 for $J=2$. Hence, we can obtain the MMSE according to .[^9]

[**Complex measurement matrix:**]{} Consider a complex $\F^{\mathcal{C}}=\F^{\mathcal{R}}+i\F^{\mathcal{I}}\in\mathbb{C}^{M\times N}$ with entries $F_{\mu l}^{\mathcal{R}}, F_{\mu l}^{\mathcal{I}}\sim \mathcal{N}(0,\frac{1}{2N})$. Expanding out the complex channel, $\y^{\mathcal{C}}=\F^{\mathcal{C}}\s^{\mathcal{C}}+\z^{\mathcal{C}}$, we obtain the equivalent real-valued SMV channel, $$\label{eq.realComplexChannel}
  \begin{bmatrix}
    \y^{\mathcal{R}} \\
    \y^{\mathcal{I}}
  \end{bmatrix}
=
  \begin{bmatrix}
    \F^{\mathcal{R}} & -\F^{\mathcal{I}}  \\
    \F^{\mathcal{I}}  & \F^{\mathcal{R}}
  \end{bmatrix}
    \begin{bmatrix}
    \s^{\mathcal{R}} \\
    \s^{\mathcal{I}}
  \end{bmatrix}
  +
    \begin{bmatrix}
    \z^{\mathcal{R}} \\
    \z^{\mathcal{I}}
  \end{bmatrix}.$$

We re-arrange  as follows, $$\label{eq:complexMatRearrange}
\underbrace{\begin{bmatrix}
    \y^{\mathcal{R}} \\
    \y^{\mathcal{I}}
  \end{bmatrix}}_{\overline{\y}}
\!=\!
\underbrace{\begin{bmatrix}
    \F_{:,1}^{\mathcal{R}},-\F_{:,1}^{\mathcal{I}},...,\F_{:,N}^{\mathcal{R}}, -\F_{:,N}^{\mathcal{I}} \\
    \F_{:,1}^{\mathcal{I}},\ \ \F_{:,1}^{\mathcal{R}},...,\F_{:,N}^{\mathcal{I}},\ \  \F_{:,N}^{\mathcal{R}}
  \end{bmatrix}}_{\overline{\F}}
    \underbrace{\begin{bmatrix}
    s_1^{\mathcal{R}}\\
    s_1^{\mathcal{I}}\\
    \vdots\\
    s_N^{\mathcal{R}}\\
    s_N^{\mathcal{I}}
  \end{bmatrix}}_{\overline{\s}}
  \!+\!
    \underbrace{\begin{bmatrix}
    \z^{\mathcal{R}} \\
    \z^{\mathcal{I}}
  \end{bmatrix}}_{\overline{\z}},$$ where $\{:\}$ refers to all the rows. In the re-arranged channel , the measurement matrix $\overline{\F}$ consists of super symbols, $$\label{eq:SMV_F}
    \overline{\F}_{\mu l}=\left\{
                \begin{array}{ll}
                 &[F_{\mu l}^{\mathcal{R}},-F_{\mu l}^{\mathcal{I}}],\ \mu\in\{1,...,M\}\\
                 &[F_{\mu l}^{\mathcal{I}}, F_{\mu l}^{\mathcal{R}}],\ \mu\in\{M+1,...,2M\}
                \end{array}
            \right.,\\$$ and the signal $\overline{\s}$ consists of $\overline{\s}_{l}=\begin{bmatrix}
    s_l^{\mathcal{R}} \\
    s_l^{\mathcal{I}}
  \end{bmatrix},\ l\in\{1,...,N\}$. The measurements and noise are $\overline{\y}=\begin{bmatrix}
    \y^{\mathcal{R}} \\
    \y^{\mathcal{I}}
  \end{bmatrix}$ and $\overline{\z}=\begin{bmatrix}
    \z^{\mathcal{R}} \\
    \z^{\mathcal{I}}
  \end{bmatrix}$, respectively. Hence, $\overline{y}_{\mu}=\sum_{l=1}^N \overline{\F}_{\mu l}\overline{\s}_l+\overline{z}_{\mu},\ \mu\in\{1,...,2M\}$.

Section \[sec:proof\] shows that the free energy and MMSE for SMV complex CS with complex measurement matrices are the same as MMV-1 with $J=2$. Note that in the free energy expression  of MMV-1, the MSE, $D=E$ , is the average MSE of the $J$ entries of $\s_l$. Therefore, in this complex CS setting, $D$ is the average MSE of the real and imaginary parts of the signal entries.

Proof of Lemma \[lemma:covIsSame\] {#sec:proof}
==================================

In this section, we show that the quantity $\mathbb{X_{\mu}}$ is the same for MMV-1 and MMV-2. Moreover, we show that complex SMV with a complex measurement matrix also yields the same $\mathbb{X_{\mu}}$ with $J=2$.

First, we re-write  in the vector form $$\label{eq:XMuVector}
\mathbb{X}_{\mu}\!=\!\mathbb{E}_{\v_{\mu}}\!\left[\operatorname{e}^{-\frac{1}{2\Delta}{\sum_{j=1}^{J}}{\sum_{a=1}^{n}}(v_{\mu j}^a)^2}\right]\!
=\!\mathbb{E}_{\v_{\mu}}\!\left[\operatorname{e}^{-\frac{1}{2\Delta}\v_{\mu}^T\v_{\mu}}\right],$$ where $\v_{\mu}=[v_{\mu 1}^1,...,v_{\mu 1}^a,...,v_{\mu J}^1$, $...,v_{\mu J}^n]^T$ and $v_{\mu j}^a$ is given in . In order to calculate the expectation w.r.t. $\v_{\mu}$ in , we calculate the distribution of $\v_{\mu}$, which is approximated by a Gaussian distribution, due to the central limit theorem. The mean is $\mathbb{E}_{\F,\z}[v_{\mu j}^a]=0$.

We now calculate the covariance matrix, $\G_{\mu}=\mathbb{E}[\v_{\mu}\v_{\mu}^T]$. The matrix is separated into $J\times J$ blocks of size $n\times n$, as shown in Fig. \[fig.cov\]. The main diagonal of $\G_{\mu}$ consists of entries $w_1=\mathbb{E}_{\F,\z}[(v_{\mu j}^a)^2]$. The entries in the blocks along the main diagonal (other than entries along the main diagonal itself) are $w_3=\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu j}^b]$. The main diagonals of other blocks have entries $w_2=\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu\eta}^a]$, and other entries in these blocks are $w_4=\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu \eta}^b]$. We now calculate each of these values as follows for MMV-1, MMV-2, and complex SMV with a complex measurement matrix.

![Covariance matrix $\G_{\mu}\in\mathbb{R}^{nJ\times nJ}$. Each block in $\G_{\mu}$ has a size of $n\times n$. The entries in the heavily marked blocks take the value $w_3$, except that entries along the dashed diagonal are $w_1$. The entries in the lightly marked blocks take the value $w_4$, except that entries along the dotted diagonal are $w_2$.[]{data-label="fig.cov"}](MMV_cov_mat.png){width="3.5cm"}

[**MMV-1:**]{} We begin by calculating the diagonal entries of the covariance matrix $\G_{\mu}=\mathbb{E}[\v_{\mu}\v_{\mu}^T]$, $$\label{eq:vVar1}
\begin{split}
    &w_1=\mathbb{E}_{\F,\z}\!\l[(v_{\mu j}^a)^2\r]\!=\!{\sum_{l,k=1}^{N,N}}\!\Bigg[(\s_l-\x_l^a)^T\times\\
    & \mathbb{E}_{\F}\l\{\F_{\mu+M(j-1), l}^T\F_{\mu+M(j-1), k}\r\}(\s_k-\x_k^a)\Bigg]+\Delta.
\end{split}$$ In , $\mathbb{E}_{\F}\l\{\F_{\mu+M(j-1), l}^T\F_{\mu+M(j-1), k}\r\}=\frac{\delta_{k,l}}{N}\widetilde{\mathbf{I}}_J$ (cf. Fig. \[fig:channel\]), where $\widetilde{\mathbf{I}}_J$ is a $J\times J$ matrix with only one 1 located at row $j$ and column $j$, and $\delta_{k,l}=1$ when $k=l$, else zero. Hence,  becomes $$\begin{aligned}
 w_1&=& \mathbb{E}_{\F,\z}\l[(v_{\mu j}^a)^2\r]=\frac{1}{N}{\sum_{l=1}^{N}} (s_{l,j}-x_{l,j}^a)^2+\Delta\label{eq:vVar2}\\
  &=&\!\frac{1}{NJ}\!{\sum_{l=1}^{N}} \!(\s_l-\x_l^a)^T (\s_l-\x_l^a)\!+\!\Delta,\label{eq:vVar2_1}\end{aligned}$$ where $s_{l,j}$ and $x_{l,j}^a$  denote the $j$-th entries in super symbols $\s_l$ and $\x_l^a$, respectively, and  holds because all $J$ entries within the same super symbol ($\s_l$ or $\x_l^a$) are i.i.d.

Similarly, we obtain $$\begin{split}
w_2=&\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu\eta}^{a}] = \frac{1}{N}{\sum_{l=1}^{N}}(s_{l,j}-x_{l,j}^a)(s_{l,\eta}-x_{l,\eta}^{a})\\
     &= \frac{1}{NJ}{\sum_{l=1}^{N}}(\s_l-\x_l^a)^T (\s_l^a-\x_l^b),\label{eq:sx_iid1}
\end{split}$$ where entries of $\s_l^a$ follow the same distribution as entries of $\s_l$ given $l$, and  is due to ([*i*]{}) entries of $\s_l$ being i.i.d., ([*ii*]{}) entries of $\x_l^{\{\cdot\}}$ being i.i.d. for fixed $l$, and ([*iii*]{}) the replica symmetry assumption [@Krzakala2012probabilistic; @krzakala2012statistical]. We also obtain $$\label{eq:v_j_eta_a}
\begin{split}
w_3=\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu j}^b]\!&=\!\frac{1}{NJ}\!{\sum_{l=1}^{N}}(\s_l\!-\!\x_l^a)^T \!(\s_l-\x_l^b)+\Delta.\\
w_4=\mathbb{E}_{\F,\z}[v_{\mu j}^a v_{\mu\eta}^{b}]&=\frac{1}{NJ}{\sum_{l=1}^{N}}(\s_l-\x_l^a)^T (\s_l^a-\x_l^b),
\end{split}$$

We now define the following auxiliary parameters $$\label{eq:auxParamsSet1}
\begin{split}
m_a=\frac{\displaystyle{\sum_{l=1}^{N}} (\x_l^a)^T\s_l}{NJ},&\quad Q_a=\frac{\displaystyle{\sum_{l=1}^{N}} (\x_l^a)^T\x_l^a}{NJ}, \\
q_{ab}=\frac{\displaystyle{\sum_{l=1}^{N}} (\x_l^a)^T\x_l^b}{NJ},&\quad
q_0=\frac{1}{NJ}{\sum_{l=1}^{N}}(\s_l^a)^T \s_l,
\end{split}$$ which allow us to express – as $$w_1=\rho-2m_a+Q_a+\Delta,$$ $$\label{eq:ws2}
w_2=q_0-(m_a+m_b)+q_{ab},$$ $$w_3 = \rho-(m_a+m_b)+q_{ab}+\Delta,$$ $$\label{eq:ws4}
w_4=q_0-(m_a+m_b)+q_{ab}.$$

Plugging the distribution of $\v_{\mu}$, approximated by $P(\v_{\mu})=[(2\pi)^n\det (\G_{\mu})]^{-\frac{1}{2}}\exp(-\frac{1}{2}\v_{\mu}^T\G_{\mu}^{-1}\v_{\mu})$, into , we obtain $$\label{eq:Xmu2}
\mathbb{X}_{\mu}=\l[\det(\mathbb{I}_n+\frac{1}{\Delta}\G_{\mu})\r]^{-1/2}.$$

[**MMV-2:**]{} For the matrix $\F$  in this setting, rows $jM+1,...,(j+1)M,\ 2\leq j \leq J$, will be the right-shift of rows $(j-1)M+1,...,jM$. We express $v_{\mu j}^a$  as $$\label{eq:v_mu_j_a_mmv2}
  v_{\mu j}^a={\sum_{l=1}^{N}}\F_{\mu l}{\bf T}_{j}(\s_l-\x_l^a)+z_{\mu+M(j-1)},\ \mu\in \{1,\cdots,M\},$$ where $\T_{j}$ is a $J\times J$ transform matrix with the $j$-th entry of the first row being one and all other entries in $\T_j$ being zeros. Using the same derivations as in MMV-1, it can be proved that the covariance matrix $\G_{\mu}=\mathbb{E}[\v_{\mu}\v_{\mu}^T]$ in MMV-2 is identical to that of MMV-1. Therefore, $\mathbb{X}_{\mu}$ in MMV-1 and MMV-2 are identical in the large system limit.

[**Complex SMV with complex measurement matrix:**]{} The derivations are the same as in MMV-2 above, except that we need to change $\F_{\mu l}$ in  to $\overline{\F}_{\mu l}$  and replace $\T_j$ by $$\label{eq:TransMat}
{\bf T}=\begin{bmatrix}
    0 & 1 \\
    -1 & 0
  \end{bmatrix},$$ because $\overline{\F}_{(\mu+M)l}=\overline{\F}_{\mu l}{\bf T},\ \mu\in\{1,...,M\}$. Using similar steps as above, we obtain that the covariance matrix $\G_{\mu}$ in this case, is also the same as that of MMV-1 with $J=2$.

**Solving $\mathbb{X}_{\mu}$**: For such a structured matrix $\G_{\mu}$ (Fig. \[fig.cov\]), elementary transforms show that the eigen-values (EV’s) are comprised of one EV equal to $\alpha_1=[w_1+(J-1)w_2]+(n-1)[w_3+(J-1)w_4],\ (J-1)$ EV’s equal to $\alpha_2=(w_1-w_2)+(n-1)(w_3-w_4),\ (n-1)$ EV’s equal to $\alpha_3=[w_1+(J-1)w_2]-[w_3+(J-1)w_4]$, and $(J-1)(n-1)$ EV’s equal to $\alpha_4=(w_1-w_2)-(w_3-w_4)$.

Owing to replica symmetry [@Krzakala2012probabilistic; @krzakala2012statistical], we have $m_a=m_b=m$, $Q_a=Q$, and $q_{ab}=q$. Also, in the Bayesian setting, we have $m=q_0=q$ and $Q=\rho$. Thus, $w_2=w_4=0$ ( and ), and $$\label{eq:detSet2}
\begin{split}
\det (\mathbb{I}_{nJ}+&\frac{1}{\Delta}\G_{\mu})= \l(1+\frac{\alpha_1}{\Delta}\r)\l(1+\frac{\alpha_2}{\Delta}\r)^{J-1}\times\\
&\quad \l(1+\frac{\alpha_1}{\Delta}\r)^{n-1}\l(1+\frac{\alpha_1}{\Delta}\r)^{(n-1)(J-1)}\\
&=\l(1+n\frac{w_3}{\Delta+\alpha_4}\r)^J\!\l(1+\frac{1}{\Delta}\alpha_4\r)^{Jn}\!.
\end{split}$$ Considering , we simplify , $$\label{eq:XmuNew}
\lim_{n\rightarrow 0}\mathbb{X}_{\mu}=\operatorname{e}^{-\frac{nJ}{2}\l[\frac{\rho-2m+\Delta+q}{Q-q+\Delta}+\log(Q-q+\Delta)-\log(\Delta)\r]},$$ where we rely on the following Taylor series, $$\label{eq:firstOrder1}
\operatorname{e}^{nk}\approx 1+nk\Rightarrow \operatorname{e}^{-\frac{n}{2}k}\approx (1+nk)^{-1/2},\ n\rightarrow 0.$$

Numerical Results {#sec:numeric}
=================

Given a free energy expression for a CS problem, the MMSE can be obtained by evaluating the largest free energy [@Tanaka2002; @GuoVerdu2005; @Krzakala2012probabilistic; @krzakala2012statistical; @MezardMontanariBook; @Barbier2015]. Having derived the free energy for the two MMV settings in Section \[sec:main\], this section calculates the MMSE under various cases. Different performance regions of MMV are identified, where the MMSE behaves differently as a function of the noise variance $\Delta$ and measurement rate $R$. We identify a phase transition of belief propagation (BP) that separates regions where BP is optimal asymptotically or not. Simulation results match the predicted performance of BP.

Performance regions: Definitions and numerical results {#sec:PerfRegion}
------------------------------------------------------

When calculating the MMSE  for different settings from the free energy expression , four different [*performance regions*]{} will appear, as illustrated in Fig. \[fig:freeEnergyProf\] and discussed below.

![Free energy as a function of MSE for different measurement rates $R$ (number of jointly sparse signal vectors $J=3$ and noise variance $\Delta=-35$ dB). The black circles mark the largest free energy, and so they correspond to the MMSE.[]{data-label="fig:freeEnergyProf"}](freeEnergyProfiles_J3.pdf){width="8cm"}

[**Regions 1 and 4:**]{} The free energy  has one local maximum point w.r.t. the MSE $D$ . This $D$ leads to the globally maximum free energy and is the MMSE.

[**Regions 2 and 3:**]{} There are 2 local maxima in the free energy, $D_1$ and $D_2$, where $D_1<D_2$. In Region 2, the smaller MSE, $D_1$, leads to the larger local maximum free energy  (hence, $\mathcal{F}(D_1)$ is the global maximum), and is the MMSE. In Region 3, the larger MSE, $D_2$, is the MMSE.

[**Boundaries between regions:**]{} We denote the boundary separating regions 1 and 2 by the [*BP threshold*]{} $R_{BP}(\Delta)$, the boundary separating regions 2 and 3 by the [*low noise threshold*]{} $R_l(\Delta)$, and the boundary separating regions 3 and 4 by the [*critical threshold*]{} $R_c(\Delta)$.

[**Numerical results:**]{} Consider $J$-dimensional Bernoulli-Gaussian signals  with sparsity rate $\rho=0.1$. Evaluating the free energy  with the noise variance $\Delta$ from -20 dB to -50 dB and measurement rate $R$ from 0.11 to 0.24, we obtain the MMSE as a function of $\Delta$ and $R$ for $J=1,3$, and $5$, as shown in Fig. \[fig:PerformanceRegions\].[^10] The darkness of the shades represents the natural logarithm of the MMSE, $\ln$(MMSE). In all panels, the critical threshold $R_c(\Delta)$, low noise threshold $R_l(\Delta)$, and BP threshold $R_{BP}(\Delta)$, as well as Regions 1-4, are marked.

In Regions 3 and 4, the best-possible algorithm yields a large MMSE for all noise variances. In contrast, in Regions 1 and 2, the optimal algorithm yields an MMSE that decreases with the noise variance $\Delta$. To summarize, the optimal algorithm yields poor estimation performance below the low noise threshold $R_{l}(\Delta)$, and good performance above $R_{l}(\Delta)$.

We further examine the MMSE as a function of the number of jointly sparse signal vectors $J$ and the measurement rate $R$. We plotted the MMSE in dB scale in Fig. \[fig:MMSE\_R\_J\]. The noise variance is -35 dB. We can see that the MMSE decreases with more signal vectors $J$ and greater measurement rate $R$. However, the MMSE depends less on $J$ as $J$ is increased. Note that the discontinuity in the MMSE surface in Fig. \[fig:MMSE\_R\_J\] is a result of the different performance regions that the various settings (different $J$ and $R$) lie in.

BP phase transition {#sec:phaseTrans}
-------------------

Belief propagation (BP) [@DMM2009; @CSBP2010; @Montanari2012; @Bayati2011; @Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015] is an algorithmic framework motivated by statistical physics, which can often achieve the optimal estimation performance (MMSE). When there are multiple local maxima $D_1<D_2$ in the free energy , BP converges to the local maximum with the larger MSE, $D_2$ [@DMM2009; @Montanari2012; @Bayati2011; @Krzakala2012probabilistic; @krzakala2012statistical]. Hence, $D_2$ characterizes the [*predicted MSE*]{} for BP. Moving from Region 1 to Region 2 by decreasing the measurement rate $R$ with fixed noise variance $\Delta$, the number of local maxima increases from 1 to 2. Therefore, BP estimation performance experiences a sudden deterioration (increase in MSE) when the measurement rate $R$ drops such that the combination of the noise variance $\Delta$ and measurement rate $R$ moves from Region 1 to Region 2. The BP threshold, $R_{BP}(\Delta)$, is the boundary between Regions 1 and 2, and is where the BP phase transition happens. That is, BP achieves poor estimation performance below $R_{BP}(\Delta)$, and good performance above $R_{BP}(\Delta)$.

[**Remark 2**]{}: In Fig. \[fig:PerformanceRegions\], we see that increasing $J$ reduces the BP threshold $R_{BP}(\Delta)$. Since BP achieves the MMSE when $R>R_{BP}(\Delta)$, increasing $J$ is beneficial to applications that use BP as the estimation algorithm.

[**Remark 3**]{}: We further analyzed the low noise ($\Delta\rightarrow 0$) and zero noise ($\Delta=0$) cases. The critical threshold $R_c(\Delta)$ converges to $\rho$ as the noise variance $\Delta$ is decreased for $J=1,3$, and $5$. We believe that this numerical result holds for every $J$. Moreover, this result matches the theoretical robust threshold of Wu and Verd[ú]{} [@WuVerdu2012] for $J=1$ in the low noise limit. Our numerical results also show that the BP threshold $R_{BP}(\Delta)$ converges to some value for different $J$ as $\Delta\rightarrow 0$. Analyzing these observations rigorously is left for future work.

![MMSE in dB as a function of measurement rate $R$ and number of jointly sparse signal vectors $J$ (noise variance $\Delta=-35$ dB).[]{data-label="fig:MMSE_R_J"}](MMSE_vs_R_J.pdf){width="8cm"}

BP simulation {#sec:AMPsim}
-------------

After obtaining the theoretic MMSE for MMV, as well as the predicted MSE for BP, we run some simulations to estimate the $\underline{\s}^j$ of channel  in a Bayesian setting. The algorithm we use is approximate message passing (AMP) [@DMM2009; @Montanari2012; @Bayati2011; @Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015], which is an approximation to the BP algorithm; related algorithms have been proposed by Ziniel and Schniter [@ZinielSchniter2013MMV] and Kim et al. [@KimChangJungBaronYe2011]. In the SMV case, when the measurement matrix and the signal have i.i.d. entries, AMP has the state-evolution (SE) property [@DMM2011; @Bayati2011; @JavanmardMontanari2012; @Donoho2013; @Bayati2015] that tracks the evolution of the MSE at each iteration. Recently, Javanmard and Montanari proved that SE tracks AMP rigorously in an SMV setting with a spatially coupled measurement matrix [@JavanmardMontanari2012]. According to our transform in Fig. \[fig:channel\], we can see that the proof [@JavanmardMontanari2012] could be extended to the MMV setting. Note 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 index.\
  $\delta=\frac{1}{NJ}\sum_{l=1}^N\sum_{j=1}^J(\widehat{a}^j_l-a^j_l)^2$ // Change in estimate\
[**Outputs:**]{} Estimate $\a^j,\forall j$

We simulated the signals in  with $J=3$ signal vectors and sparsity rate $\rho=0.1$ measured by a channel  with measurement rate $R\in[0.11,0.24]$ and noise variance $\Delta\in[-20,-50]$ dB. For each setting, we generated 50 signals of length $N=5000$, and the resulting MSE compared to the predicted BP MSE is shown in Fig. \[fig:AMPoverMSE\].[^13]

![AMP simulation results ($\text{MSE}_{\text{AMP}}$) compared to the predicted BP MSE ($\text{MSE}_{\text{BP}}$) with $J=3$ jointly sparse signal vectors. The dashed curve, solid curve, and the curve comprised of little circles correspond to thresholds $R_c(\Delta),\ R_l(\Delta)$, and $R_{BP}(\Delta)$, respectively. Regions 1-4 are also marked. The shade denotes $\ln \l(\frac{\text{MSE}_{\text{AMP}}}{\text{MSE}_{\text{BP}}}\r)$, which we expect to be 0 (completely dark shades) in the entire $R$ versus $\Delta$ plane. The narrow bright band above the BP threshold indicates the mismatch of AMP simulated MSE to the BP predicted MSE.[]{data-label="fig:AMPoverMSE"}](simMSE_over_predMSE_J3_TSP_v2.pdf){width="8cm"}

The labels of the thresholds are omitted for brevity. We can see that AMP simulation results match with the predicted MSE of BP and BP phase transition from the replica analysis of Section \[sec:phaseTrans\]. Note that there is a narrow band of light shades above the BP threshold, $R_{BP}(\Delta)$ (the top threshold), meaning that the simulated MSE is greater than the predicted MSE; this is due to randomness in our generated signals and channels. Note that we also compared the AMP simulation results to that of the M-SBL algorithm [@YeKimBresler2015], a widely used algorithm to solve the MMV problem. The M-SBL results were not as good. Indeed, because AMP is often an achievable to the MMSE, other algorithms are expected to provide greater MSE.

Conclusion {#sec:conclusion}
==========

We analyzed the minimum mean squared error (MMSE) for two settings of multi-measurement vector (MMV) problems, where the entries in the signal vectors are independent and identically distributed (i.i.d.), and share the same support. One MMV setting has i.i.d. Gaussian measurement matrices, while the other MMV setting has identical i.i.d. Gaussian measurement matrices. Replica analyses yield identical free energy expressions for these two settings in the large system limit when the signal length goes to infinity and the number of measurements scales with the signal length. Because of the identical free energy expressions, the MMSE’s for both MMV settings are identical. By numerically evaluating the free energy expression, we identified different performance regions for MMV where the MMSE as a function of the channel noise variance and the measurement rate behaves differently. We also identified a phase transition for belief propagation algorithms (BP) that separates regions where BP achieves the MMSE asymptotically and where it is suboptimal. Simulation results of an approximated version of BP matched with the MSE predicted by replica analysis. As a special case of MMV, we extended our replica analysis to single measurement vector (SMV) complex CS, so that we can calculate the MMSE for SMV complex CS with real or complex measurement matrices.

Acknowledgments {#acknowledgments .unnumbered}
===============

The work in this paper is based in part on preliminary work with Jong Min Kin, Woohyuk Chang, Bangchul Jung, and Jong Chul Ye [@KimChangJungBaronYe2011]. The authors thank Lenka Zdeborov[á]{} for useful discussions about replica analysis, and Yanting Ma and Ryan Pilgrim for helpful comments. Junan Zhu also thanks Shikai Luo for helpful discussions.

This appendix follows the derivation of Barbier and Krzakala [@Barbier2015], except for some nuances. Our compressed derivation makes the presentation self-contained.

Plugging  and the following identity [@Barbier2015; @Krzakala2012probabilistic], $$\begin{split}
  &1=\int \exp\Bigg\{-{\sum_{a=1}^{n}} \l[\widehat{m}_a\l(m_a NJ-{\sum_{l=1}^{N}}(\x_l^a)^T\s_l\r)\r]+{\sum_{a=1}^{n}}\Bigg[\\
  &\widehat{Q}_a\!\l(Q_a\frac{NJ}{2}\!-\!\frac{1}{2}{\sum_{l=1}^{N}}(\x_l^a)^T\x_l^a\r)\Bigg]\!-\!{\sum_{1\leq a< b\leq n}^{}}\Bigg[\widehat{q}_{ab}\Bigg(q_{ab} NJ-\\
  &{\sum_{l=1}^{N}}(\x_l^a)^T\x_l^b\Bigg)\Bigg]\!\Bigg\}\!{\prod_{a=1}^{n}}dQ_a\  d\widehat{Q}_{a}\  dm_a \ d\widehat{m}_{a}\! {\prod_{1\leq a<b\leq n}^{}}\! dq_{ab}\ d\widehat{q}_{ab},
\end{split}$$ into , we obtain $$\label{eq:EZn3}
\begin{split}
  \mathbb{E}_{\F,\s,\z}&[Z^n]\!=\!(2\pi\Delta)^{-\frac{nMJ}{2}}\bigintsss\! \exp\Bigg[NJ\Bigg(\frac{1}{2}{\sum_{a=1}^{n}}\widehat{Q}_aQ_a\\
 & -\frac{1}{2}{\sum_{\substack{1\leq a,b\leq n\\a\neq b}}^{}}\widehat{q}_{ab}q_{ab}-{\sum_{a=1}^{n}}\widehat{m}_am_a\Bigg)\Bigg]\l[{\prod_{\mu=1}^{M}} \mathbb{X}_{\mu}\!\r]\times\\
  &\Gamma^N{\prod_{a=1}^{n}}dQ_a\  d\widehat{Q}_{a}\  dm_a\  d\widehat{m}_{a}\!{\prod_{\substack{1\leq a,b\leq n\\a\neq b}}^{}}\! dq_{ab}\ d\widehat{q}_{ab},
\end{split}$$ where $$\label{eq:Gamma_original}
\begin{split}
&\Gamma=\!\!\!\int\! P(\s_1)\! \l[{\prod_{a=1}^{n}}P(\x^a_1)\r]\!\exp\!\Bigg[\!-\frac{1}{2}{\sum_{a=1}^{n}}\widehat{Q}_a(\x^a_1)^T\x^a_1+\\
&\frac{1}{2}{\sum_{\substack{1\leq a,b\leq n\\a\neq b}}^{}}\widehat{q}_{ab}(\x^a_1)^T\x^b_1+{\sum_{a=1}^{n}}\widehat{m}_a(\x^a_1)^T\s_1\Bigg]d\s_1 {\prod_{a=1}^{n}}d\x^a_1.
\end{split}$$

$$\label{eq:PhiJ}
\begin{split}
  \widetilde{\Phi}_J&(m,\widehat{m},q,\widehat{q},Q,\widehat{Q})=\frac{J}{2}(Q\widehat{Q}+q\widehat{q}-2m\widehat{m})-\frac{MJ}{2N}\l[\frac{\rho-2m+\Delta+q}{Q-q+\Delta}+\log(Q-q+\Delta)-\log(\Delta)\r]+\\
  &\int P(\s_1) \l\{ \int \log \l\{ \int P(\x_1)\exp\l[-\frac{1}{2}(\widehat{Q}+\widehat{q})\x_1^T\x_1+\widehat{m}\x^T_1\s_1+
  \sqrt{\widehat{q}}\h^T\x_1\r]d\x_1\r\} \mathcal{D}\h \r\} d\s_1-\frac{MJ}{2N}\log(2\pi\Delta).
\end{split}$$

**Further simplification of **: The Stratanovitch transform [@Stratanovitch-Wiki] in $J$ dimensions is given by $$\label{eq:strat}
\begin{split}
  &\exp\!\l[\frac{\widehat{q}}{2}\!{\sum_{\substack{1\!\leq\! a,b\!\leq n\\a\neq b}}^{}}\!(\x^a_1)^T\!\x^b_1\!\r]\!=\!{\prod_{j=1}^{J}}\!\exp\!\l[\frac{\widehat{q}}{2}{\sum_{\substack{1\leq a,b\leq n\\a\neq b}}^{}}x^a_{1,j} x^b_{1,j}\r]\\
  &={\prod_{j=1}^{J}}\int\exp\l[\sqrt{\widehat{q}}h_j{\sum_{a=1}^{n}}x^a_{1,j}-\frac{\widehat{q}}{2}{\sum_{a=1}^{n}}(x^a_{1,j})^2\r]\mathcal{D}h_j\\
  &=\int\exp\l[\sqrt{\widehat{q}}\h^T{\sum_{a=1}^{n}}\x^a_1-\frac{\widehat{q}}{2}{\sum_{a=1}^{n}}(\x^a_1)^T\x^a_1\r]\mathcal{D}\h,
\end{split}$$ where $\h=[h_1,...,h_J]^T$, and the differential $\mathcal{D}h_j=\frac{1}{\sqrt{2\pi}}\operatorname{e}^{-h_j^2/2}d h_j$. With the Stratanovitch transform , we simplify $\Gamma$  as follows, $$\label{eq:Gamma}
\Gamma=\int P(\s_1)\int\l[f(\h)\r]^n \mathcal{D}\h\ d\s_1,$$ where $f(\h)=\int P(\x_1)\operatorname{e}^{-\frac{\widehat{Q}+\widehat{q}}{2}\x_1^T\x_1+\widehat{m}\x_1^T\s_1+\sqrt{\widehat{q}}\h^T\x_1}d\x_1$, and we drop the super-script $a$ of $\x_1^a$ owing to the replica symmetry assumption [@Krzakala2012probabilistic; @krzakala2012statistical]. In the limit of $n\rightarrow 0$, using another Taylor series $[f(\h)]^n\approx 1+n\log [f(\h)]$, we have $\int [f(\h)]^n \mathcal{D}\h\approx 1+n\int \log [f(\h)] \mathcal{D}\h\approx\operatorname{e}^{n\int \log [f(\h)]\mathcal{D}\h}$, so that $\mathbb{E}\{\int[f(\h)]^n\mathcal{D}\h\}\approx\mathbb{E}\{1+n\int \log [f(\h)] \mathcal{D}\h\}\approx\operatorname{e}^{\mathbb{E}\{n\int \log [f(\h)]\mathcal{D}\h\}}$. Hence, we can approximate  as $$\label{eq:Gamma1}
\Gamma=\exp\l\{n\int P(\s_1) \int\log [f(\h)] \mathcal{D}\h \ d\s_1\r\}.$$ Considering , we rewrite as $$\label{eq:EZn4}
    \mathbb{E}_{\F,\s,\z}[Z^n]=\int \operatorname{e}^{nN\widetilde{\Phi}_J(m,\widehat{m},q,\widehat{q},Q,\widehat{Q})}
      dm\ d\widehat{m}\ dq\ d\widehat{q}\ dQ\ d\widehat{Q},$$ where $\widetilde{\Phi}_J(m,\widehat{m},q,\widehat{q},Q,\widehat{Q})$ is given in .

**Free energy expression**: We now substitute  into . Assuming that the limits in  commute and that we only evaluate  at optimum points of $\widetilde{\Phi}_J$  [@Barbier2015; @Krzakala2012probabilistic; @krzakala2012statistical], we have $\mathcal{F}=\widetilde{\Phi}_J(m^*,\widehat{m}^*,q^*,\widehat{q}^*,Q^*,\widehat{Q}^*)$, where the asterisks denote stationary points. Next, we calculate the stationary points: $$\frac{\partial \widetilde{\Phi}_J }{\partial m}=0 \Rightarrow \widehat{m}^*=\frac{R}{Q^*-q^*+\Delta},$$ $$\frac{\partial \widetilde{\Phi}_J }{\partial q}=0 \Rightarrow \widehat{q}^*=R\frac{\Delta+\rho-2m^*+q^*}{(Q^*-q^*+\Delta)^2},$$ $$\frac{\partial \widetilde{\Phi}_J }{\partial Q}=0 \Rightarrow \widehat{Q}^*=R\frac{2m^*-\rho-2q^*+Q^*}{(Q^*-q^*+\Delta)^2},$$ where $R$  is the measurement rate. Because we are analyzing the MMSE, we must assume that the estimated prior matches the true underlying prior, which is a Bayesian setting. Thus, $q^*=m^*$ and $Q^*=\rho$ . Let $E=q^*-2m^*+Q^*=Q^*-q^*$, then we obtain $\widehat{q}^*=\widehat{m}^*=\frac{R}{E+\Delta}$ and $\widehat{Q}^*=0$. Therefore, we solve for the free energy as a function of $E$ in . Using a change of variables, we obtain , which is a function of $E$. Using , the MSE is $$\label{eq:DandE}
D=E+Q-q=E+\rho/N\overset{N\rightarrow \infty}{\longrightarrow} E.$$ Hence, in the large system limit, we can regard the free energy  as a function of the MSE, $D$.

[10]{}

E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: [E]{}xact signal reconstruction from highly incomplete frequency information,” , vol. 52, no. 2, pp. 489–509, Feb. 2006.

D. Donoho, “Compressed sensing,” , vol. 52, no. 4, pp. 1289–1306, Apr. 2006.

R. G. Baraniuk, “A lecture on compressive sensing,” , vol. 24, no. 4, pp. 118–121, July 2007.

H. Jung, J. C. Ye, and E. Y. Kim, “[Improved k-t [BLAST]{} and k-t [SENSE]{} using [FOCUSS]{}]{},” , vol. 52, no. 11, pp. 3201–3226, May 2007.

H. Jung, K. Sung, K. S. Nayak, E. Y. Kim, and J. C. Ye, “[k-t [FOCUSS]{}: a general compressed sensing framework for high resolution dynamic [MRI]{}]{},” , vol. 61, no. 1, pp. 103–116, Jan. 2009.

S. F. Cotter and B. D. Rao, “Sparse channel estimation via matching pursuit with application to equalization,” , vol. 50, no. 3, pp. 374–377, Mar. 2002.

J. Ma and F. X. Le Dimet, “[Deblurring from highly incomplete measurements for remote sensing]{},” , vol. 47, no. 3, pp. 792–802, Mar. 2009.

J. Haupt and R. Nowak, “[S]{}ignal reconstruction from noisy random projections,” , vol. 52, no. 9, pp. 4036–4048, Sept. 2006.

M. F. Duarte, M. B. Wakin, D. Baron, S. Sarvotham, and R. G. Baraniuk, “Measurement bounds for sparse signal ensembles via graphical models,” , vol. 59, no. 7, pp. 4280–4289, July 2013.

G. J. Pottie and W. J. Kaiser, “[Wireless integrated network sensors]{},” , vol. 43, no. 5, pp. 51–58, May 2000.

M. F. Duarte, M. B. Wakin, D. Baron, and R. G. Baraniuk, “[Universal distributed sensing via random projections]{},” in [*Proc. Inf. Process. Sensor Networks (IPSN)*]{}, Nashville, TN, Apr. 2006.

J. Chen and X. Huo, “Theoretical results on sparse representations of multiple measurement vectors,” , vol. 54, no. 12, pp. 4634–4643, Dec. 2006.

S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiple measurement vectors,” , vol. 53, no. 7, pp. 2477–2488, July 2005.

M. Mishali and Y. C. Eldar, “Reduce and boost: recovering arbitrary sets of jointly sparse vectors,” , vol. 56, no. 10, pp. 4692–4702, Oct. 2009.

E. Berg and M. P. Friedlander, “Joint-sparse recovery from multiple measurements,” , Apr. 2009.

O. Lee, J. M. Kim, Y. Bresler, and J. C. Ye, “Compressive diffuse optical tomography: [N]{}oniterative exact reconstruction using joint sparsity,” , vol. 30, no. 5, pp. 1129–1142, May 2011.

K. Lee, Y. Bresler, and M. Junge, “Subspace methods for joint sparse recovery,” , vol. 58, no. 6, pp. 3613–3641, June 2012.

J. C. Ye, J. M. Kim, and Y. Bresler, “Improving [M-SBL]{} for joint sparse recovery using a subspace penalty,” , vol. 63, no. 24, pp. 6595–6605, Dec. 2015.

J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. part [I]{}: Greedy pursuit.,” , vol. 86, no. 3, pp. 572–588, Mar. 2006.

D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” , vol. 53, no. 8, pp. 3010–3022, Aug. 2005.

J. A. Tropp, “Algorithms for simultaneous sparse approximation. part [II]{}: Convex relaxation,” , vol. 86, no. 3, pp. 589–602, Mar. 2006.

J. Ziniel and P. Schniter, “Efficient message passing-based inference in the multiple measurement vector problem,” in [*Proc. IEEE 45th Asilomar Conf. Signals, Syst., and Comput.*]{}, Nov. 2011, pp. 1447–1451.

T. Tanaka, “[A statistical-mechanics approach to large-system analysis of CDMA multiuser detectors]{},” , vol. 48, no. 11, pp. 2888–2910, Nov. 2002.

D. Guo and S. Verd[ú]{}, “Randomly spread [CDMA]{}: [A]{}symptotics via statistical physics,” , vol. 51, no. 6, pp. 1983–2010, June 2005.

A. Montanari and D. Tse, “Analysis of belief propagation for non-linear problems: The example of [CDMA]{} (or: How to prove [T]{}anaka’s formula),” in [*IEEE Inf. Theory Workshop*]{}, Mar. 2006, pp. 160–164.

F. Krzakala, M. M[é]{}zard, F. Sausset, Y. Sun, and L. Zdeborov[á]{}, “Probabilistic reconstruction in compressed sensing: [A]{}lgorithms, phase diagrams, and threshold achieving matrices,” , vol. 2012, no. 08, pp. P08009, Aug. 2012.

F. Krzakala, M. M[é]{}zard, F. Sausset, Y. Sun, and L. Zdeborov[á]{}, “Statistical-physics-based reconstruction in compressed sensing,” , vol. 2, no. 2, pp. 021005, May 2012.

M. M[é]{}zard and A. Montanari, , Oxford University press, 2009.

J. Barbier and F. Krzakala, “Approximate message-passing decoder and capacity-achieving sparse superposition codes,” , Mar. 2015.

T. Lesieur, F. Krzakala, and L. Zdeborov[á]{}, “Phase transitions in sparse [PCA]{},” in [*Proc. Int. Symp. Inf. Theory (ISIT)*]{}, July 2015, pp. 1635–1639.

D. L. Donoho, A. Maleki, and A. Montanari, “[Message passing algorithms for compressed sensing]{},” , vol. 106, no. 45, pp. 18914–18919, Nov. 2009.

D. Baron, S. Sarvotham, and R. G. Baraniuk, “Bayesian compressive sensing via belief propagation,” , vol. 58, pp. 269–280, Jan. 2010.

M. Bayati and A. Montanari, “The dynamics of message passing on dense graphs, with applications to compressed sensing,” , vol. 57, no. 2, pp. 764–785, Feb. 2011.

A. Montanari, “Graphical models concepts in compressed sensing,” , pp. 394–438, 2012.

D. Guo and C. C. Wang, “Multiuser detection of sparsely spread [CDMA]{},” , vol. 26, no. 3, pp. 421–431, Apr. 2008.

J. Barbier, M. Dia, N. Macris, and F. Krzakala, “The mutual information in random linear estimation,” , July 2016.

G. Reeves and H. D. Pfister, “The replica-symmetric prediction for compressed sensing with gaussian matrices is exact,” , July 2016.

J. Ziniel and P. Schniter, “Efficient high-dimensional inference in the multiple measurement vector problem,” , vol. 61, no. 2, pp. 340–354, Jan. 2013.

J. Zhu and D. Baron, “Performance regions in compressed sensing from noisy measurements,” in [*Proc. 2013 Conf. Inference Sci. Syst. (CISS)*]{}, Baltimore, MD, Mar. 2013.

Y. Wu and S. Verd[ú]{}, “Optimal phase transitions in compressed sensing,” , vol. 58, no. 10, pp. 6241 – 6263, Oct. 2012.

J. Kim, W. Chang, B. Jung, D. Baron, and J. C. Ye, “Belief propagation for jointly sparse recovery,” , Feb. 2011.

D. L. Donoho, A. Maleki, and A. Montanari, “The noise-sensitivity phase transition in compressed sensing,” , vol. 57, no. 10, pp. 6920–6941, Oct. 2011.

A. Javanmard and A. Montanari, “State evolution for general approximate message passing algorithms, with applications to spatial coupling,” , Dec. 2012.

D. Donoho, I. Johnstone, and A. Montanari, “Accurate prediction of phase transitions in compressed sensing via a connection to minimax denoising,” , vol. 59, no. 6, pp. 3396–3433, June 2013.

M. Bayati, M. Lelarge, and A. Montanari, “Universality in polytope phase transitions and message passing algorithms,” , vol. 25, no. 2, pp. 753–822, Feb. 2015.

“Hubbard–[S]{}tratonovich transformation,” https://en.wikipedia.org/wiki/\
Hubbard%E2%80%93Stratonovich\_transformation.

[^1]: The work was supported in part by the National Science Foundation under the Grant CCF-1217749, the U.S. Army Research Office under the Contract W911NF-14-1-0314, and the European Research Council under the European Union’s 7th Framework Programme (FP/2007- 2013)/ERC Grant Agreement 307087-SPARCS.

[^2]: Junan Zhu and Dror Baron are with the Department of Electrical and Computer Engineering, NC State University, Raleigh, NC 27695. E-mail: {jzhu9, barondror}@ncsu.edu.

[^3]: Florent Krzakala is with Sorbonne Universit[é]{}s, Universit[é]{} Pierre et Marie Curie Paris 6 and Ecole Normale Superieure, 75005 Paris, France. E-mail: florent.krzakala@ens.fr.

[^4]: Part of the literature [@Tanaka2002; @GuoVerdu2005] defines the free energy as the negative of , so that fixed points of the free energy correspond to local minima.

[^5]: Recently, the replica Gibbs free energy has been proven rigorously for the SMV case by Barbier et al. [@BDMK2016] and Reeves and Pfister [@ReevesPfister2016]. We conjecture that by generalizing these two works [@BDMK2016; @ReevesPfister2016], our MMV analysis can be made rigorous; we leave it for future work.

[^6]: The $J$-dimensional integrals in  can be simplified to one-dimensional integrals using a change of coordinate to $J$-sphere coordinate. Note also that $E$ approaches the MSE in the large system limit; details appear in the appendix.

[^7]: In Section \[sec:complex\], we only deal with SMV CS, and omit the word “SMV."

[^8]: A replica analysis for complex CS with a real measurement matrix appears in Guo and Verd[ú]{} [@GuoVerdu2005]. Their derivation does not cover complex matrices.

[^9]: As a reminder, the free energy of MMV-2 is identical to that of MMV-1 in the large system limit.

[^10]: The MMV with $J=1$ becomes an SMV. The MMSE results in Fig. \[fig:MMV\_J1\] match with the SMV MMSE in Krzakala et. al. [@Krzakala2012probabilistic; @krzakala2012statistical] and Zhu and Baron [@ZhuBaronCISS2013].

[^11]: When the assumptions about the measurement matrix and signal [@DMM2009; @Montanari2012; @Bayati2011; @Krzakala2012probabilistic; @krzakala2012statistical; @Barbier2015] are violated, AMP might suffer from convergence issues.

[^12]: Note that Algorithm \[algo:AMP\_MMV\] is a straightforward simplification of the AMP algorithm in Barbier and Krzakala [@Barbier2015].

[^13]: We simulated both $J$ different measurement matrices $\underline{F}^j$ and $J$ identical $\underline{F}^j$. Both results match the predicted BP MSE, which support our conclusion that the MMSE’s of both settings are the same. Fig. \[fig:AMPoverMSE\] is with $J$ different $\underline{F}^j$.
