---
abstract: 'Real-world applications such as magnetic resonance imaging with multiple coils, multi-user communication, and diffuse optical tomography often assume a linear model where several sparse signals sharing common sparse supports are acquired by several measurement matrices and then contaminated by noise. Multi-measurement vector (MMV) problems consider the estimation or reconstruction of such signals. In different applications, the estimation error that we want to minimize could be the mean squared error or other metrics such as the mean absolute error and the support set error. Seeing that minimizing different error metrics is useful in MMV problems, we study information-theoretic performance limits for MMV signal estimation with arbitrary additive error metrics. We also propose a message passing algorithmic framework that achieves the optimal performance, and rigorously prove the optimality of our algorithm for a special case. We further conjecture the optimality of our algorithm for some general cases, and back it up through numerical examples. As an application of our MMV algorithm, we propose a novel setup for active user detection in multi-user communication and demonstrate the promise of our proposed setup.'
author:
- 'Junan Zhu,  and Dror Baron,  [^1] [^2]'
title: 'Performance Limits with Additive Error Metrics in Noisy Multi-Measurement Vector Problems'
---

[*Keywords*]{}:Active user detection, error metric, message passing, multi-measurement vector problem.

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

Many systems in science and engineering can be approximated by a linear model, where a signal $\x \in \mathbb{R}^N$ is recorded via a measurement matrix $\A\in\mathbb{R}^{M \times N}$, and then contaminated by a measurement channel, $$\label{eq:SMV}
\w = \A \x,\ y_m = \mathcal{Z}(w_m), \forall m\in \{1,\ldots,M\},$$ where $y_m,m\in\{1,\ldots,M\}$, are the entries of the measurements $\y\in\mathbb{R}^M$, and the measurement channel $\mathcal{Z}(\cdot)$ is characterized by a probability density function (pdf), $f(y_m|w_m)$. The goal is to estimate $\x$ from the measurements $\y$ given knowledge of $\A$ and a model for the measurement channel $f(y_m|w_m),\forall m$. We call such a system the [*single measurement vector (SMV)*]{} problem.

In many applications, the signal acquisition systems are distributed, where $J$ measurement matrices measure $J$ different signals individually. The key difference between such a system and $J$ individual SMV’s, is that these $J$ signals are somewhat dependent. An example of a model containing such dependencies is the multi-measurement vector (MMV) problem [@chen2006trs; @cotter2005ssl; @Mishali08rembo; @Berg09jrmm; @LeeKimBreslerYe2011; @LeeBreslerJunge2012; @YeKimBresler2015]. The MMV problem considers the estimation of a set of dependent signals, and has applications such as magnetic resonance imaging with multiple coils [@JuYeKi07; @JuSuNaKiYe09], active user detection in multi-user communication [@FletcherRanganGoyal2009; @Boljanovic2017], and diffuse optical tomography using multiple illumination patterns [@LeeKimBreslerYe2011]. In MMV, thanks to the dependencies among different signals, 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], if the underlying signals share the same sparse supports. A non-rigorous replica analysis of MMV with measurement noise also shows the benefit of having more signal vectors [@ZhuBaronKrzakala2017IEEE; @ZhuDissertation2017].

[**Related work:**]{} There are many estimation approaches for MMV problems. These include 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]. However, these algorithms cannot handle the case of $J$ [*different*]{} measurement matrices. Statistical approaches [@ZinielSchniter2011] often achieve the oracle minimum mean squared error (MMSE). However, when running estimation algorithms for MMV problems, one might be interested in minimizing some other error. For example, if estimating the underlying signal is important, one could use the mean squared error (MSE) metric; when there might be outliers in the estimated signal, using the mean absolute error (MAE) metric might be more appropriate. Seeing that there is no prior work discussing the optimal performance with user-defined error metrics, we study the optimal performance with user-defined additive error metrics in MMV problems where the signals share common sparse supports, and each entry of the measurements is contaminated by parallel measurement channels (i.e., the channel in  satisfies $f(\y|\w)=\prod_{m=1}^M f(y_m|w_m)$). Note that a specific error metric, the MSE, has been studied in Zhu et al. [@ZhuBaronKrzakala2017IEEE], which focuses on the MSE performance limits (i.e., the MMSE) of MMV signal estimation. In contrast, this work explores performance limits and designs an algorithm that can minimize [*arbitrary additive error metrics*]{} beyond MSE.

[**Contributions:**]{} This paper combines insights from Zhu et al. [@ZhuBaronKrzakala2017IEEE] and Tan and coauthors [@Tan2012SSP; @Tan2014], thus yielding a stronger understanding of the MMV problem, which could be extended in future work to other distributed signal acquisition settings, beyond MMV. To be more specific, we make several contributions in this paper. First, by extending Tan and coauthors [@Tan2012SSP; @Tan2014] we provide an algorithm based on a message passing (MP) framework [@RanganGAMP2011ISIT; @Krzakala2012probabilistic] that can be adapted to minimize the expected error for arbitrary additive error metrics. Our algorithm first runs MP until it converges or reaches some stopping criteria, and we then denoise MP’s output using a denoiser that minimizes the given additive error metric. Second, we prove rigorously that our algorithm is optimal in a specific SMV case,[^3] and further conjecture the optimality of our algorithm in MMV. Third, as an example, we derive performance limits for MAE and mean weighted support set error (MWSE) by designing the corresponding optimal denoisers, based on the scalar channel noise variance (Section \[sec:MPA\]) derived from replica analysis (Appendix \[app:inverMMSE\]) [@ZhuBaronKrzakala2017IEEE]. Simulation results show the superiority and optimality of our algorithm. We note in passing that having more signal vectors in MMV helps reduce the estimation error. Finally, as an application of MMV and our algorithm, we propose a novel setup for active user detection in multi-user communications (details in Section \[sec:app\]) and demonstrate the promise of our proposed setup through simulation.

[**Organization:**]{} The remainder of the paper is organized as follows. We introduce our problem setting and MP algorithms in Section \[sec:background\]. Our algorithmic framework, which can minimize arbitrary additive error metrics, is proposed in Section \[sec:achievable\]; we rigorously prove the optimality of our algorithm for an SMV case and conjecture the optimality of our algorithm for MMV cases in Section \[sec:converse\]. For some example error metrics, we derive the corresponding optimal algorithms, together with the theoretical limits for these errors, in Section \[sec:example\]. Synthetic simulation results are discussed in Section \[sec:numeric\_synth\], followed by an application of our metric-optimal algorithm to a real-world problem in Section \[sec:app\]. We conclude in Section \[sec:conclusion\].

[**Notations:**]{} In this paper, bold capital letters represent matrices, bold lower case letters represent vectors, and normal font letters represent scalars. The $m$-th entry (scalar) of a vector $\z$ is denoted by $z_m$.

Problem Setting and Background {#sec:background}
==============================

Problem setting
---------------

[**Signal model**]{}: We consider an ensemble of $J$ signal vectors, $\x^{(j)}\in\mathbb{R}^N,\ j\in\{1,\ldots,J\}$, where $j$ is the index of the signal. Consider a [*super-symbol*]{} $\x_n=[x_n^{(1)},\ldots,x_n^{(J)}],\ n\in\{1,\ldots,N\}$; all super-symbols in this paper are row vectors. The super symbol $\x_n$ follows a $J$-dimensional distribution, $$\label{eq:jsm}
f(\x_n)=\rho \phi(\x_n)+(1-\rho)\delta(\x_n),$$ where $\rho\in (0,1)$ determines the percentage of non-zeros in the signal and is called the sparsity rate, $\phi(\x_n)$ is a $J$-dimensional pdf, and $\delta(\x_n)=\left\{
                \begin{array}{ll}
                 1, & \x_n = \textbf{0},\\
                 0, & \text{else}.
                \end{array}
            \right.\\
$

\[def:jointly\_sparse\] [*Ensembles of signals are called jointly sparse signals with common sparse supports when they obey .*]{}

Note that there are other types of joint sparsity [@BaronDCStech] that fit into the MMV framework. For example, an MMV model with signal vectors that have slowly changing supports is discussed in Ziniel and Schniter [@ZinielSchniter2013MMV]. Since this paper only focuses on the MMV problem with signals sharing common sparse supports , we refer to joint sparsity with common sparse supports as joint sparsity for brevity.

[**Measurement models**]{}: Each signal $\x^{(j)}$ is measured by a measurement matrix $\A^{(j)}\in\mathbb{R}^{M\times N}$ before being corrupted by a random measurement channel, $$\label{eq:linearMixing}
\begin{split}
\w^{(j)}&=\A^{(j)}\x^{(j)},\ y_m^{(j)}=\mathcal{Z}(w_m^{(j)}),\\
 &\quad m\in\{1,\ldots,M\},\ j\in\{1,\ldots,J\},
\end{split}$$ where $y_m^{(j)}, m\in \{1,\ldots,M\}$ are the entries of the measurements $\y^{(j)}$, and the measurement channel $\mathcal{Z}(\cdot)$ is characterized by the pdf $f(y_m^{(j)}|w_m^{(j)})$. In this paper, we only focus on independent and identically distributed (i.i.d.) parallel measurement channels, i.e., the pdf’s $f(y_m^{(j)}|w_m^{(j)}), \forall m,j$, are identical and there is no cross-talk among different channels; our proposed algorithm is readily extended to parallel channels with different $f(y_m^{(j)}|w_m^{(j)}), \forall m,j$. When the number of signal vectors $J=1$, we call this MMV model an SMV problem .

\[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$,*]{} $$\lim_{N\rightarrow\infty} \frac{M(N)}{N} = R>0.$$

We call $R$ the measurement rate.

For MMV problems, we are given the matrices $\A^{(j)}$ and measurements $\y^{(j)},\ \forall j$, as well as knowledge of the measurement channel . Our task is to estimate the underlying signal vectors $\x^{(j)},\ \forall j$. Suppose that the estimate is $\widehat{\x}^{(j)}$. Define $\X=[\x^{(1)},\cdots,\x^{(J)}]$ and $\widehat{\X}=[\widehat{\x}^{(1)},\cdots,\widehat{\x}^{(J)}]$. Therefore, $\X=[\x_1^T,\cdots,\x_N^T]^T$, where $\{\cdot\}^T$ denotes transpose. The estimation quality is quantified by a user-defined error metric $D_{\text{UD}}(\X,\widehat{\X})$, where the subscript UD denotes “user-defined." We define this [*additive error metric $D_{\text{UD}}(\cdot,\cdot)$ as*]{} $$D_{\text{UD}}(\X, \widehat{\X})={\sum_{n=1}^{N}}d_{\text{UD}}(\x_n,\widehat{\x}_n),$$ where $d_{\text{UD}}(\cdot,\cdot): \mathbb{R}^J\times \mathbb{R}^J \rightarrow \mathbb{R}$ is an arbitrary user-defined error metric on each super-symbol. The smaller the $D_{\text{UD}}(\cdot,\cdot)$ is, the better the estimation is.

Message passing algorithms {#sec:MPA}
--------------------------

Message passing (MP) algorithms consider a factor graph [@RanganGAMP2011ISIT; @Krzakala2012probabilistic; @BarbierKrzakala2017IT], which expresses the relation between the signals $\x$ and measurements $\y$. We begin by discussing the factor graph for SMV, followed by that of MMV.

![Factor graph for SMV (left) and MMV (right).[]{data-label="fig:factorGraph"}](factorGraph.png){width="8cm"}

[**Factor graph for SMV:**]{} The left panel of Fig. \[fig:factorGraph\] illustrates the factor graph concept for an SMV problem  with i.i.d. entries in the signal $\x$. The round circles are the variable nodes (representing the distribution of the signal), and the squares denote the factor nodes (representing the measurement channel). The variables $x_n,\ \forall n$, are driven by each factor node $f(x_n)$ individually, because the signal has i.i.d. entries $x_n$. There are two types of messages passed in the factor graph shown on the left panel of Fig. \[fig:factorGraph\]: the message passed by the variable node $x_n$ to the factor node $y_m$, $\mathcal{M}_{n\rightarrow m}(x_n)$, and the message passed by the factor node $y_m$ to the variable node $x_n$, $\widehat{\mathcal{M}}_{m\rightarrow n}(x_n)$. According to the literature on MP algorithms [@RanganGAMP2011ISIT; @Krzakala2012probabilistic; @BarbierKrzakala2017IT], we have the following relation: $$\begin{split}
\mathcal{M}_{n\rightarrow m}(x_n) &=\frac{1}{Z_{n\rightarrow m}}f(x_n) \prod_{\widehat{m} \neq m} \widehat{\mathcal{M}}_{\widehat{m}\rightarrow n}(x_n), \\
\widehat{\mathcal{M}}_{m\rightarrow n}(x_n) &=\frac{1}{Z_{m\rightarrow n}} \int f(y_m | \x) \prod_{\widehat{n}\neq n} \mathcal{M}_{\widehat{n}\rightarrow m}(x_{\widehat{n}})\prod_{\widehat{n}\neq n}dx_{\widehat{n}},
\end{split}$$ where we use a single integral sign to denote a multi-dimensional integration for brevity, and $Z_{n\rightarrow m}$ and $Z_{m\rightarrow n}$ are normalization factors. The aim of this paper is not the detailed derivation of MP algorithms. Instead, the key property of MP utilized in this paper is that MP converts  into the following equivalent scalar channel, $$\q=\x+\v,$$ where $\q$ is the noisy [*pseudo data*]{}, and $\v$ is the equivalent scalar channel additive white Gaussian noise (AWGN) whose variance $\Delta_v$ can be approximated. After obtaining $\q$ and $\Delta_v$, each variable node $x_n$ updates the estimate $\widehat{x}_n$ by denoising $q_n$. If the signal entries are i.i.d. and the matrix is either i.i.d. or sparse and locally tree-like, then MP algorithms yield a density function $f(x_n|q_n)$ that is statistically equivalent to $f(x_n|\y)$ [@RanganGAMP2011ISIT].

[**Factor graph for MMV:**]{} An MMV problem with jointly sparse signals can be expressed as the factor graph shown in the right panel of Fig. \[fig:factorGraph\]. We can see that the MP in each channel is similar to an SMV problem. The only difference is that the $J$ variable nodes $x_n^{(j)},\ j\in\{1,\ldots,J\}$, for fixed $n$ are driven by one factor node $f(\x_n)$. By grouping the entries from different signal vectors together into super-symbols as in , we have i.i.d. super-symbols. The noisy super-symbol pseudo data $\q_n=[q_n^{(1)},\ldots,q_n^{(J)}]$ is denoised, in order to update the estimate for $\x_n=[x_n^{(1)},\ldots,x_n^{(J)}]$.

Main Results {#sec:main}
============

We first present our metric-optimal algorithm in Section \[sec:achievable\] and then in Section \[sec:converse\] we rigorously prove that our metric-optimal algorithm is optimal in the SMV case under certain conditions.[^4] Based on our proof for the SMV case, we conjecture that the proposed algorithm is optimal for arbitrary additive error metrics in the MMV case.

Achievable part: Metric-optimal estimation algorithm {#sec:achievable}
----------------------------------------------------

The metric-optimal algorithm consists of two parts, as illustrated in Algorithm \[algo:metric\_opt\_MMV\]. We first run an MP algorithm (Algorithm \[algo:AMP\_MMV\] provides an implementation of an MP algorithm) to get the noisy pseudo data $\q^{(j)}, \forall j$, and the noise variance $\Delta_v$ (details below). Next, we denoise $\q^{(j)}, \forall j$, using an optimal denoiser tailored to minimize the given error metric. The following discusses both parts in detail.

\
[**Inputs:**]{} Measurements $\y^{(j)}$ and matrices $\A^{(j)}, \forall j$\
[**Part 1 (Algorithm \[algo:AMP\_MMV\]):**]{} Obtain pseudo data $\q^{(j)},\ \forall j$, and scalar channel noise variance $\Delta_v$ from MP($\y^{(j)}, \A^{(j)}, \forall j$)\
[**Part 2 (examples in Section \[sec:example\]):**]{} Obtain optimal estimate $\widetilde{\x}^{(j)}$ from denoiser using $\q^{(j)},\Delta_v, \forall j$\
[**Outputs:**]{} $\widetilde{\x}^{(j)},\ \forall n$

\
[**Inputs:**]{} Maximum number of iterations $t_{\text{max}}$, threshold $\epsilon$, sparsity rate $\rho$, noise variance $\Delta_z$, measurements $\y^{(j)}$, and measurement matrices $\A^{(j)}, \forall j$\
[**Initialize:**]{} $t=1,\delta=\infty,\k^{(j)}=\y^{(j)},\Theta_m^{(j)}=0,s^{(j)}_n=\rho\Delta_z,\widehat{x}_n^{(j)}=0, h_m^{(j)}=0, \forall m,n,j$ \[line:beginForLoop\]\
$\boldsymbol{\Theta}^{(j)}=(\A^{(j)})^2 \s^{(j)}$\[line:Theta\]\
$\k^{(j)}=\A^{(j)} \widehat{\x}^{(j)}-\text{diag}(\boldsymbol{\Theta}^{(j)}) \h^{(j)}$\
$h_m^{(j)}=g_{\text{out}}\l(k_m^{(j)},y_m^{(j)},\Theta_m^{(j)}\r)$\[line:g\_out\]\
$r_m^{(j)} = -\frac{\partial}{\partial k_m^{(j)}} g_{\text{out}}\l(k_m^{(j)},y_m^{(j)},\Theta_m^{(j)}\r)$\[line:deriv\_g\_out\]\
// Scalar channel noise variance\
$\Delta_v^{(j)}=\l\{\frac{1}{N} \mathbf{1}^T \l[(\A^{(j)})^T\r]^2 \mathbf{r}^{(j)}\r\}^{-1}$ \[line:scalarNoiseVar\]\
$\q^{(j)}=\widehat{\x}^{(j)}+\Delta_v^{(j)} (\A^{(j)})^T \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,
$$ $$f_{v_n}(\Delta_v,\q_n)\!=\!-[f_{a_n}(\Delta_v,\q_n)]^2+\frac{\rho}{C(\Delta_v+1)}\!\l[\frac{\q_n^2}{\Delta_v+1}\!+\!\Delta_v\r],$$ where $\q_n^2=\l[\l(q_n^{(1)}\r)^2,\ldots,\l(q_n^{(J)}\r)^2\r]$ and $$C=\rho+(1-\rho)\l(1+\frac{1}{\Delta_v}\r)^{\frac{J}{2}}\exp\l[-\frac{\q_n\q_n^T}{2\Delta_v(\Delta_v+1)}\r].$$ Notice that in Line \[line:scalarNoiseVar\] we take the mean of a vector to obtain a scalar $\Delta_v^{(j)}$, which is the average of the variances for the estimates of signal entries $x_n^{(j)}$. This is because the super-symbols $\x_n, \forall n$, of the signals are i.i.d and the $J$ measurement channels are i.i.d. For the same reason, $\Delta_v^{(j)},\ j\in\{1,\ldots,J\}$, should be close to each other; hence, Line \[line:mean\_delta\]. Note that Algorithm \[algo:AMP\_MMV\] assumes that the entries of $\A^{(j)}$ scale with $\frac{1}{\sqrt{N}}$, and is a more generic form of an algorithm from our prior work with Krzakala [@ZhuBaronKrzakala2017IEEE].

[**Metric-optimal denoiser:**]{} The second part of our metric-optimal algorithm takes as inputs the noisy pseudo data $\q_n=\x_n+\v_n$ and the estimated variance of $\v_n,\ \Delta_v$, from the MP algorithm. Using Bayes’ rule, we can derive the posterior $f(\x_n|\q_n)$ and use $f(\x_n|\q_n)$ to formulate the optimal estimator in the sense of the user-defined error metric. The optimal estimate is $$\label{eq:estimator}
\widetilde{\x}_n=\arg\min_{\widehat{\x}_n} \int d_{\text{UD}}(\x_n,\widehat{\x}_n) f(\x_n|\q_n) d\x_n.
$$

Converse part: The optimal estimate {#sec:converse}
-----------------------------------

The reason why  is optimal is based on the insight from Rangan [@RanganGAMP2011ISIT] that in SMV the density function $f(x_n|q_n)$ converges to the posterior $f(x_n|\y)$. A rigorous proof for a certain SMV case is provided below, followed by our conjecture that  is optimal in the MMV case.

\[lemma:optSMV\] Consider an SMV  in the large system limit (Definition \[def:largeSystemLimit\]) with an AWGN measurement channel, $f(y_m|w_m)=\frac{1}{\sqrt{2\pi\sigma_Z^2}}\exp\l[\frac{(y_m-w_m)^2}{2\sigma_Z^2}\r], \forall m$. The estimate $\widetilde{\x}_n$  is optimal in the sense that $$\label{eq:optSMV}
\lim_{N\rightarrow\infty} \frac{1}{N}\sum_{n=1}^N d_{\text{UD}}(\widetilde{x}_n,x_n) = \text{MUDE},$$ where MUDE denotes the minimum user-defined error, if all the conditions below hold.

1.  The entries of the measurement matrix are i.i.d. Gaussian, $A_{mn}\sim \mathcal{N}(0,\frac{1}{N})$,

2.  the signal entries are i.i.d. with bounded fourth moment $\mathbb{E}[X^4]<B$, where $B$ is some constant,

3.  the free energy given by replica analysis has one fixed point [@ZhuBaronCISS2013; @ZhuBaronKrzakala2017IEEE; @Krzakala2012probabilistic],[^6]

4.  the user-defined error metric $d_{\text{UD}}(\cdot,\cdot)$ is pseudo-Lipschitz [@Bayati2011],[^7]

5.  the optimal estimator  as a function of $q_n,\ \widetilde{x}_n(q_n): \mathbb{R}\rightarrow \mathbb{R}$ , is Lipschitz continuous, and

6.  Part 1 converges before entering Part 2 in Algorithm \[algo:metric\_opt\_MMV\].

According to Theorem 1 in Bayati and Montanari [@Bayati2011], we know that $\lim_{N\rightarrow\infty} \frac{1}{N}\sum_{n=1}^N d_{\text{UD}}(\widetilde{x}_n(q_n),x_n) = \mathbb{E}\l[d_{\text{UD}}(\widetilde{x}_{n}(X+\Delta_v Z),X)\r]$, where $X$ is a random variable following the same distribution as the signal entries, $Z\sim \mathcal{N}(0,1)$, and $\Delta_v$ is given in Line \[line:mean\_delta\] of Algorithm \[algo:AMP\_MMV\]. The question boils down to showing that $\text{MUDE} = \mathbb{E}\l[d_{\text{UD}}(\widetilde{x}_{n}(X+\Delta_v Z),X)\r]$. In fact, when Algorithm \[algo:AMP\_MMV\] converges, $\Delta_v$ corresponds to the scalar channel noise variance associated with the MMSE given by replica analysis [@Bayati2011; @RushVenkataramanan2016]. In addition, the MMSE provided by replica analysis is proved to be exact under the conditions asserted in Lemma \[lemma:optSMV\] [@ReevesPfister2016]. Therefore, $\mathbb{E}\l[d_{\text{UD}}(\widetilde{x}_{n}(X+\Delta_v Z),X)\r]$ is indeed the MUDE. Hence, we proved , and the estimator  is optimal.

**Remark 1:** The optimality of the metric-optimal estimator $\widetilde{\x}_n$  is stated in the sense that the ensemble mean user-defined error converges almost surely to the MUDE. It is possible that there exist other estimators that achieve the MUDE.

**Remark 2:** The proof is made possible by linking three rigorous proofs [@ReevesPfister2016; @Bayati2011; @RushVenkataramanan2016] from the prior art. That said, numerical examples in Section \[sec:numeric\_synth\] demonstrate that Algorithm \[algo:metric\_opt\_MMV\] yields promising results even if the conditions required by Lemma \[lemma:optSMV\] are not met.

After proving the optimality of our metric-optimal algorithm in the SMV scenario with certain conditions, we state the following conjecture.

\[conj:2\] In the large system limit, for the MMV model with the signal in  and the user-defined additive error metric $d_{\text{UD}}(\x_n,\widehat{\x}_n)$, the optimal estimate of the signal vectors is $$\widetilde{\x}_n=\arg\min_{\widehat{\x}_n} \mathbb{E}[d_{\text{UD}}(\x_n,\widehat{\x}_n)|\q_n].$$

As the reader can see from the proof of Lemma \[lemma:optSMV\], in order to rigorously prove Conjecture \[conj:2\], we need to show ([*i*]{}) the MMSE given by the replica analysis for the MMV case [@ZhuBaronKrzakala2017IEEE] is exact, ([*ii*]{}) the scalar channel noise variance $\Delta_v$ in Algorithm \[algo:AMP\_MMV\] corresponds to the MMSE given by replica analysis, and ([*iii*]{}) a result similar to Theorem 1 in Bayati and Montanari [@Bayati2011] holds. None of these three results exists in the prior art, so we do not foresee what exact conditions are needed for Conjecture \[conj:2\]. Proving these three results (and hence Conjecture \[conj:2\]) is beyond the scope of this work. Instead, we provide some intuition that explains why we believe Algorithm \[algo:metric\_opt\_MMV\] is optimal in the MMV scenario.

In the SMV case, $f(x_n|q_n)$ converges to the posterior $f(x_n|\y)$ in relaxed BP [@RanganGAMP2011ISIT]. As an extension to the MMV case, we intuitively think that $f(\x_n|\q_n)$ would converge to the posterior $f(\x_n|\{\y^{(j)}\}_{j=1}^J)$, based on two observations: ([*i*]{}) the measurement channels are i.i.d., so that it suffices to update the estimate of the channel by passing the messages $\mathcal{M}_{n\rightarrow m}(x_n^{(j)})$ and $\widehat{\mathcal{M}}_{m\rightarrow n}(x_n^{(j)})$ for different $j$ individually, and ([*ii*]{}) the super-symbols $\x_n$ are i.i.d., so that denoising each super-symbol individually accounts for all the information needed to update the estimate.

The optimal estimate of each super-symbol $\x_n$ in the signal vectors is $$\label{eq:estimatorTrue}
\widetilde{\x}_{\text{true},n}=\arg\min_{\widehat{\x}_n} \int d_{\text{UD}}(\x_n,\widehat{\x}_n) f(\x_n|\{\y^{(j)}\}_{j=1}^J) d\x_n.$$ Comparing  and , provided that $f(\x_n|\q_n)$ converges to the posterior $f(\x_n|\{\y^{(j)}\}_{j=1}^J)$, we have $$\widetilde{\x}_n=\arg\min_{\widehat{\x}_n} \mathbb{E}[d_{\text{UD}}(\x_n,\widehat{\x}_n)|\q_n]\approx \widetilde{\x}_{\text{true},n},$$ which results in Conjecture \[conj:2\].

Example Metric-optimal Estimators and Performance Limits {#sec:example}
========================================================

In order to derive metric-optimal estimators, we need to know the scalar channel noise variance, $\Delta_v$. Below we obtain the optimal estimator, which is then used as part of Algorithm \[algo:metric\_opt\_MMV\] and to evaluate the performance limits of our metric-optimal algorithm.

For i.i.d. matrices and AWGN channels , replica analysis in our previous work with Krzakala [@ZhuBaronKrzakala2017IEEE] yields the information-theoretic scalar channel noise variance $\Delta_v$ for message passing algorithms, which characterizes the posterior $f(\x_n|\q_n)$.[^8] Hence, we can obtain the information-theoretic optimal performance with arbitrary additive error metrics. For other types of matrices, our replica analysis [@ZhuBaronKrzakala2017IEEE] does not hold. Nevertheless, MP algorithms still yield a posterior $f(\x_n|\q_n)$, which we conjecture converges to the true posterior $f(\x_n|\{\y^{(j)}\}_{j=1}^J)$ (Conjecture \[conj:2\]). Hence, we can still assume that the scalar channel noise variance $\Delta_v$ is known. Based on the known variance $\Delta_v$, we build metric-optimal estimators  for two examples, mean weighted support set error (Section \[sec:weightedSupportSet\]) and mean absolute error (Section \[sec:MAE\]).

Mean weighted support set error (MWSE) {#sec:weightedSupportSet}
--------------------------------------

[**MWSE-optimal estimator:**]{} In support set estimation, the goal is to estimate the support of the signal, which is 1 if the corresponding entry in the signal is non-zero and 0 if it is zero. There are two types of errors in support set estimation: false alarms (support is 0, but estimated to be 1) and misses (support is 1, but estimated to be 0). In some applications such as medical imaging and radar detection, a miss may mean that the doctor misses an illness of the patient, or the radar misses an incoming missile. Hence, the cost paid for a miss could be tremendous compared to a false alarm. There are other applications where a false alarm is more costly than a miss. For example, in court, if an innocent person is mistakenly judged guilty, he/she will likely suffer a great deal. Therefore, we should weight these two errors differently in different applications. Let $b_n$ and $\widehat{b}_n$ be the true support and the estimated support of the $n$-th entry of the signal, respectively, and $\beta\in[0,1]$ is an application-dependent weight, which reflects the trade-off between the false alarms and misses. Hence, the MWSE given the pseudo data $\q_n$ is $$\label{eq:MWSE}
\begin{split}
    &\text{MWSE}|\q_n=\mathbb{E}[d_{\text{WSE}}(b_n,\widehat{b}_n)|\q_n]=\\
    &\left\{
                \begin{array}{ll}
                 (1-\beta)\Pr(b_n=1|\q_n),&\   \widehat{b}_n=0\ \text{and}\ b_n=1,\\
                 \beta\Pr(b_n=0|\q_n),&\   \widehat{b}_n=1\ \text{and}\ b_n=0,\\
                 0,\quad\quad\quad&\ \widehat{b}_n=b_n,
                \end{array}
            \right.\\
\end{split}$$ where $\Pr(\cdot)$ denotes probability. The optimal estimate $\widetilde{b}_n$ minimizes $\mathbb{E}[d_{\text{WSE}}(b_n,\widehat{b}_n)|\q_n]$ , which implies $$\label{eq:weighted_optB}
    \widetilde{b}_n\!=\!\left\{
                \begin{array}{ll}
                 \!0,&\!   (1-\beta)\Pr(b_n=1|\q_n)\!\leq\! \beta\Pr(b_n=0|\q_n),\\
                 \!1,&\! (1-\beta)\Pr(b_n=1|\q_n)\!>\! \beta \Pr(b_n=0|\q_n).
                \end{array}
            \right.\\$$

Since $f(\q_n|b_n=0)=(2\pi\Delta_v)^{-\frac{J}{2}}\exp\l(-\frac{\q_n\q_n^T}{2\Delta_v}\r)$ and $f(\q_n|b_n=1)=[2\pi(\Delta_v+1)]^{-\frac{J}{2}}\exp\l[-\frac{\q_n\q_n^T}{2(\Delta_v+1)}\r]$, we have $$\Pr(b_n=1|\q_n)=\frac{\frac{\rho}{[2\pi (1+\Delta_v)]^{\frac{J}{2}}}\e^{-\frac{\q_n\q_n^T}{2(1+\Delta_v)}}}{\frac{\rho}{[2\pi (1+\Delta_v)]^{\frac{J}{2}}}\e^{-\frac{\q_n\q_n^T}{2(1+\Delta_v)}}+\frac{(1-\rho)}{(2\pi \Delta_v)^{\frac{J}{2}}}\e^{-\frac{\q_n\q_n^T}{2\Delta_v}}},$$ and $\Pr(b_n=0|\q_n)=1-\Pr(b_n=1|\q_n)$. Plugging $\Pr(b_n=1|\q_n)$ and $\Pr(b_n=0|\q_n)$ into , we have $$\label{eq:b_opt}
    \widetilde{b}_n=\left\{
                \begin{array}{ll}
                 0,&\quad   \q_n\q_n^T\leq\tau,\\
                 1,&\quad   \q_n\q_n^T>\tau,
                \end{array}
            \right.\\$$ where $\tau=2\Delta_v(1+\Delta_v)\ln\l[\frac{\beta(1-\rho)}{(1-\beta)\rho}\l(\frac{1+\Delta_v}{\Delta_v}\r)^{\frac{J}{2}}\r]$, and we remind the reader that $\rho$ is the sparsity rate.

[**Performance limits:**]{} Utilizing  and taking expectation over the pseudo data $\q_n$ for $\mathbb{E}[d_{\text{WSE}}(b_n,\widehat{b}_n)|\q_n]$ , we obtain the minimum MWSE (MMWSE), $$\label{eq:MMWSE_1}
\begin{split}
&\text{MMWSE}=\mathbb{E}[d_{\text{WSE}}(b_n,\widetilde{b}_n)]\\
&\quad =\int_{\q_n\q_n^T>\tau} \beta\Pr(b_n=0|\q_n)f(\q_n)d\q_n +\\
&\quad \quad \int_{\q_n\q_n^T\leq\tau} (1-\beta)\Pr(b_n=1|\q_n)f(\q_n)d\q_n.
\end{split}$$ We have two integrals to simplify in , where we show the first below, and the second can be obtained similarly. To derive the first integral, note that $$\label{eq:oneIntegral}
\small \int_{\q_n\q_n^T>\tau} \Pr(b_n\!=\!0|\q_n)f(\q_n)d\q_n\!=\!\Pr(\q_n\q_n^T\!>\!\tau, b_n\!=\!0).$$ Next, we calculate the pdf of the random variable (RV) $g_n=\frac{\q_n\q_n^T}{\Delta_v}$ given $b_n=0$. Because the entries of $\q_n$ are i.i.d. $\mathcal{N}(0,\Delta_v)$ given $b_n=0$, $g_n$ follows the Chi-square distribution, $f_G(g_n)=\frac{g_n^{\frac{J}{2}-1}\exp(-\frac{g_n}{2})}{2^{\frac{J}{2}}\Gamma(\frac{J}{2})}$, where $\Gamma(\cdot)$ is the Gamma function. Let $r_n=\Delta_v g_n=\q_n\q_n^T$, then we obtain $$f(r_n)=\frac{1}{\Delta_v}f_G(\frac{r_n}{\Delta_v})=\frac{r_n^{\frac{J}{2}-1}\exp\l[-\frac{r_n}{2\Delta_v}\r]}{(2\Delta_v)^{J/2}\Gamma(J/2)},$$ which helps to simplify . Therefore,  can be simplified, $$\label{eq:MMWSE_final}
\begin{split}
& \text{MMWSE}=\beta(1-\rho)\underbrace{\int_{r_n=\tau}^{\infty}\frac{r_n^{\frac{J}{2}-1}\exp\l[-\frac{r_n}{2\Delta_v}\r]}{(2\Delta_v)^{J/2}\Gamma(J/2)}dr_n}_{\Pr(\text{false alarm})}\\
& \quad +(1-\beta)\rho\underbrace{\int_{r_n=0}^{\tau}\frac{r_n^{\frac{J}{2}-1}\exp\l[-\frac{r_n}{2(1+\Delta_v)}\r]}{[2(1+\Delta_v)]^{J/2}\Gamma(J/2)}dr_n}_{\Pr(\text{miss})}.
\end{split}$$

[**Hamming distance**]{}: In digital wireless communication systems, the signal only takes discrete values. A useful error metric is the (per-entry) Hamming distance [@Cover06], which equals 1 if the estimate of an entry of the signal differs from the true value. (Section \[sec:app\] will present an example in wireless communication that minimizes the Hamming distance.) The reader can verify that Hamming distance can be interpreted as a special case of weighted support set error, where $\beta=0.5$ provides equal weight to both errors . That said, we provide more insights about this particular case, which is ubiquitous in communication systems.

For the jointly sparse model in , we define the Hamming distance as $$\label{eq:Hamming}
d_{\text{HD}}(\x_n,\widehat{\x}_n)=\mathbbm{1}_{\x_n\neq \widehat{\x}_n},$$ where $\mathbbm{1}_{\mathcal{A}}$ is the indicator function. If ([*i*]{}) the pdf $\phi(\x_n)$ in  is a $J$-dimensional Dirac-delta function $\delta(\x_n-\mathbf{1})$, where $\mathbf{1}$ is an all-one row vector, and ([*ii*]{}) the estimate satisfies $\widehat{x}_n^{\{1\}}=\cdots =\widehat{x}_n^{\{J\}}, \forall n\in \{1,\ldots,N\}$, then the weighted support set error  with weight $\beta=0.5$ is equal to half of the Hamming distance  for super symbols. In –, we briefly derive the Hamming distance-optimal estimator when $\x_n\in \{\mathbf{1},\mathbf{0}\}$, where $\mathbf{0}$ is an all-zero row vector. The mean Hamming distance (MHD) given the pseudo data $\q_n$ is $$\label{eq:MHD}
\begin{split}
    &\text{MHD}|\q_n=\mathbb{E}[d_{\text{MHD}}(\x_n,\widehat{\x}_n)|\q_n]=\\
    &\left\{
                \begin{array}{ll}
                 \Pr(\x_n=\mathbf{1}|\q_n),&\   \widehat{\x}_n=\mathbf{0}\ \text{and}\ \x_n=\mathbf{1},\\
                 \Pr(\x_n=\mathbf{0}|\q_n),&\   \widehat{\x}_n=\mathbf{1}\ \text{and}\ \x_n=\mathbf{0},\\
                 0,\quad\quad\quad&\ \widehat{\x}_n= \x_n.
                \end{array}
            \right.\\
\end{split}$$ Following the steps in –, the minimum MHD (MMHD) estimator is $$\label{eq:MHD_optTh}
\widetilde{\x}_n=\left\{
                \begin{array}{ll}
                 \mathbf{1},&   \sum_{j=1}^J q_n^{(j)} \geq\frac{J}{2}+\Delta_v \ln\l[\frac{1-\rho}{\rho}\r],\\
                 \mathbf{0},&   \sum_{j=1}^J q_n^{(j)} <\frac{J}{2}+\Delta_v \ln\l[\frac{1-\rho}{\rho}\r].
                \end{array}
            \right.\\$$

Mean absolute error (MAE) {#sec:MAE}
-------------------------

[**MAE-optimal estimator:**]{} The element-wise absolute error (AE) is $$\label{eq:AE}
  d_{\text{AE}}(x_n^{(j)},\widehat{x}_n^{(j)})=|x_n^{(j)}-\widehat{x}_n^{(j)}|.$$ In order to find the minimum mean absolute error (MMAE) estimate, $\widetilde{\x}_n$, we need to find the stationary point of , $$\label{eq:optimalCondition}
\frac{d\mathbb{E}[|x_n^{(j)}-\widehat{x}_n^{(j)}|\ | \q_n]}{d\widehat{x}_n^{(j)}}\bigg|_{\widehat{x}_n^{(j)}=\widetilde{x}_n^{(j)}}=0,$$ $\forall j\in\{1,\ldots,J\},n\in\{1,\ldots,N\}$. It can be proved that $\mathbb{E}[X]=\int_{0}^{\infty} \Pr(X>x)dx$, if $X\geq 0$. Therefore, $$\label{eq:expectation}
\begin{split}
  &\mathbb{E}[|x_n^{(j)}-\widehat{x}_n^{(j)}|\ | \q_n]=\int_{0}^{\infty} \Pr(|x_n^{(j)}-\widehat{x}_n^{(j)}|>t| \q_n)dt\\
  &=\int_{-\infty}^{\widehat{x}_n^{(j)}} \Pr(x_n^{(j)}<t| \q_n)dt+\int_{\widehat{x}_n^{(j)}}^{\infty} \Pr(x_n^{(j)}>t| \q_n)dt.
\end{split}$$ Using  and , we obtain $$\Pr(x_n^{(j)}<\widetilde{x}_n^{(j)} | \q_n)=\Pr(x_n^{(j)}>\widetilde{x}_n^{(j)}| \q_n).$$ That is, $$\label{eq:optimalEst}
  \int_{-\infty}^{\widetilde{x}_n^{(j)}} f(x_n^{(j)} | \q_n)dx_n^{(j)}=\int_{\widetilde{x}_n^{(j)}}^{\infty} f(x_n^{(j)} | \q_n)dx_n^{(j)}=\frac{1}{2},$$ through which we solve for the optimal estimator $\widetilde{x}_n^{(j)}$ numerically.

[**Performance limits:**]{} We calculate the MMAE as follows, $$\label{eq:MMAE}
\begin{split}
  \text{MMAE}\!&=\!\mathbb{E}[|\widetilde{x}_n^{(j)}-x_n^{(j)}|]\!=\!\int_{-\infty}^{\infty}\! \mathbb{E}[|\widetilde{x}_n^{(j)}-x_n^{(j)}|\ | \q_n] f(\q_n)d \q_n\\
  &=\int_{-\infty}^{\infty}\Bigg[\int_{-\infty}^{\widetilde{x}_n^{(j)}} -x_n^{(j)} f(x_n^{(j)}| \q_n)d x_n^{(j)}+\\
  &\quad\quad\quad\quad \int_{\widetilde{x}_n^{(j)}}^{\infty} x_n^{(j)}f(x_n^{(j)}| \q_n)d x_n^{(j)}\Bigg] f(\q_n)d \q_n,
\end{split}$$ which has to be numerically approximated.

![Comparison of simulation results to theoretic MMWSE for weighted support set estimation under different number of channels $J$ and measurement rates $R$ ($\beta=0.2$, noise variances $\Delta_z=0.01$).[]{data-label="fig:SimWSE"}](WeightSuppEstSimVsTheory.pdf){width="8cm"}

Synthetic Simulations {#sec:numeric_synth}
=====================

After deriving the minimum mean weighted support set error (MMWSE) and minimum mean absolute error (MMAE) estimators, this section provides numerical results for Algorithm \[algo:metric\_opt\_MMV\]. In the case of i.i.d. random matrices and AWGN channels , replica analysis yields the MMSE of the MMV problem [@ZhuBaronKrzakala2017IEEE]. By inverting the MMSE (details in Appendix \[app:inverMMSE\]), we obtain the scalar channel noise variance $\Delta_v$,[^9] which characterizes the posterior $f(\x_n|\q_n)$. Given $\Delta_v$, we characterize the MMWSE and MMAE theoretically. In the following simulations, we use i.i.d. Gaussian matrices with unit-norm rows, i.i.d. $J$-dimensional Bernoulli-Gaussian signals  with $J=1,3$, and 5, and sparsity rate $\rho=0.1$. The signal length is $N=10000$, and the measurement rate $R=\frac{M}{N}$ varies from $0.3$ to $0.7$. For each setting, the simulation results are averaged over 50 realizations of the problem.

**Mean weighted support set error in AWGN channels:** We simulate AWGN channels  in this case with the noise variance being $\Delta_z\in \{0.01,0.001\}$. Fig. \[fig:SimWSE\] shows the weighted support set estimation results using our metric-optimal algorithm compared to the MMWSE . The red dashed curve, the blue dashed-dotted curve, and the black solid curve correspond to the MMWSE of $J=1,\ 3$, and 5, respectively. The red circles, blue crosses, and black triangles represent the simulation results. We can see that our simulation results match the theoretically optimal performance.

**Remark:** The optimal weighted support set estimator  is not Lipschitz continuous. Hence, for $J=1$,  is not guaranteed to yield the MMWSE, according to Lemma \[lemma:optSMV\]. Nevertheless, numerical results for $J=1$ show that the MWSE given by  is close to the MMWSE.

For weighted support set estimation, we further study the information-theoretic optimal receiver operating characteristic (ROC) curves, which are plotted in Fig. \[fig:ROC\] for different $J$’s. The red curves and black curves are plotted for noise variances $\Delta_z=0.001$ and $0.01$, respectively. The solid curves, the dashed curves, and the dotted curves represent $J=1,\ 3$, and 5, respectively. The true positive rate (TPR) and false positive rate (FPR) in both axes are defined as TPR$=\frac{\text{\# accurately predicted positives}}{\text{\# all positives in truth}}$ and FPR$=\frac{\text{\# wrongly predicted negatives}}{\text{\# all negatives in truth}}$. We can see that having more signal vectors $J$ leads to a larger area under the ROC curve, which indicates better trade-offs between true positives and false positives.

![Receiver operating characteristic curves for weighted support set estimation under different channel noise variances $\Delta_z$ and number of channels $J$ ($\beta=0.2$, measurement rate $R=0.3$).[]{data-label="fig:ROC"}](ROCweightSuppEst.pdf){width="8cm"}

**Mean absolute error in logistic channels:** We simulate i.i.d. logistic channels  with parameters $a=10$ and $30$; the smaller $a$ is, the noisier the channel becomes. Fig. \[fig:MAEvsMMAE\] plots the simulated MAE (crosses) and the theoretic MMAE  (curves) for various settings.[^10] Different colors and line shapes refer to different $a$’s and $J$’s, respectively. We can see that our simulation results match the theoretically optimal performance.

![Comparison of simulation results to theoretic MMAE under different logistic channels , number of channels $J$, and measurement rates $R$.[]{data-label="fig:MAEvsMMAE"}](logitMAE.pdf){width="8cm"}

**Remark:** Both simulations yield better performance for larger $J$. This is intuitive, because more signal vectors that share the same support should make the estimation process easier due to more information being available.

![Simulation results of OMP and Algorithm \[algo:metric\_opt\_MMV\] for active user detection in multi-user communication (SMV).[]{data-label="fig:OMP"}](OMP_vs_metric_opt.pdf){width="8.5cm"}

Application {#sec:app}
===========

In this section, we discuss active user detection (AUD) in a multi-user communication setting that can be viewed as compressed sensing [@DonohoCS; @CandesRUP; @BaraniukCS2007] (CS, which is closely related to SMV) in some scenarios. Next, we simulate AUD using our metric-optimal algorithm. Finally, we discuss how to solve the AUD problem using an MMV setting.

[**CS based active user detection:**]{} One application of MMV is AUD in a massive random access (MRA) scenario for multi-user communication [@FletcherRanganGoyal2009; @Boljanovic2017]. In the MRA scenario, multiple end users (EUs) are requesting access to the network simultaneously by sending their unique identification codewords, $\a_n\in \{-1,+1\}^{M\times 1},\ n\in \{1,\ldots,N\}$, to the base station, where $n$ denotes the user id, and each user’s codeword is known by the base station. The base station needs to determine which EUs are requesting access to the network (active) and which are not (inactive), so that it can allocate resources to the active EUs. Fletcher et al. [@FletcherRanganGoyal2009] proposed a CS based AUD scheme for MRA, which was recently revisited by Boljanovi[ć]{} et al. [@Boljanovic2017]. In their setup [@FletcherRanganGoyal2009; @Boljanovic2017], all EUs are synchronized, i.e., all active EUs send each entry of their codewords to the base station simultaneously in one time slot. Denote the status of the $n$-th EU by $x_n\in \{0,1\}$, where $x_n=1$ means active and $x_n=0$ is inactive. Denote the received signal vector at the base station by $\y\in \mathbb{R}^{M\times 1}$. Since all EUs are synchronized, we can express the received signal $\y$ by $$\label{eq:MRA_SMV}
\y=\A\x+\z,$$ where $\A=[\a_1,\ldots,\a_N]\in \{0,1\}^{M\times N}$, and $\z$ is AWGN. The base station estimates $\x$ to determine which EUs are active.

In the following, we first apply a mean Hamming distance (MHD) optimal algorithm to estimate $\x$ from , and compare to the algorithm of Boljanovi[ć]{} et al. [@Boljanovic2017], which is orthogonal match pursuit (OMP) [@Pati1993]. Next, we propose an MMV based scheme for AUD in MRA.

[**Simulation with MHD-optimal algorithm:**]{} As the reader can see, the CS based active user detection [@FletcherRanganGoyal2009; @Boljanovic2017] is an SMV problem , which is MMV for $J=1$ . Later in this section we will extend the scheme by Boljanovi[ć]{} et al. to an MMV setting, and so we keep using MMV notations. Note that the entries of the measurement matrix in this AUD problem take values of $\pm 1$. Because the derivation of Algorithm \[algo:AMP\_MMV\] assumes that entries of $\A^{(j)}$ scale with $\frac{1}{\sqrt{N}}$,[^11] we scale $\A^{(j)}$ down by $\sqrt{N}$ using a modified $\widetilde{\y}^{(j)}=\frac{\y^{(j)}}{\sqrt{N}}$ .

Following the discussion above, we simulate the settings of measurement rate $R\in \{0.2,0.25,\ldots,0.6\}$ and noise variance $\Delta_z\in \{10^{-1},10^{-1.5},10^{-2}\}$. For each setting, we randomly generate 50 realizations of the Bernoulli signal $\x^{(1)}\in \{0,1\}^{N\times 1}$ with sparsity rate $\rho=0.1$, and measurement matrix $\A^{(1)}\in \l\{-\frac{1}{\sqrt{N}},+\frac{1}{\sqrt{N}}\r\}^{M\times N}$, where $N=10000$. We run Algorithm \[algo:metric\_opt\_MMV\] with $J=1$ to estimate the underlying signal $\x^{(1)}$. Note that $f_{a_n}(\Delta_v,\q_n)$ and $f_{v_n}(\Delta_v,\q_n)$ in Lines \[line:mean\]–\[line:var\] of Algorithm \[algo:AMP\_MMV\], which consists of Part 2 of Algorithm \[algo:metric\_opt\_MMV\], are given by the following, $$\begin{split}
f_{a_n}(\Delta_v,\q_n)&=\frac{\rho}{\rho+(1-\rho)\exp\l[-\frac{\sum_{j=1}^J q_n^{(j)} -\frac{J}{2}}{\Delta_v}\r]}\mathbf{1},\\
f_{v_n}(\Delta_v,\q_n)&=f_{a_n}(\Delta_v,\q_n)-f_{a_n}(\Delta_v,\q_n)^2,
\end{split}$$ where the power-of-two in the last term of $f_{v_n}(\cdot,\cdot)$ is applied element-wise.

Our results are compared to OMP in Fig. \[fig:OMP\]. The solid, dashed, and dash-dotted curves represent noise variance $\Delta_z=10^{-1},\ 10^{-1.5},\ 10^{-2}$, respectively. The black curves are the OMP results and the red curves with circle markers are the results of Algorithm \[algo:metric\_opt\_MMV\] when optimizing for MHD. We can see that our algorithm consistently outperforms OMP.[^12]

[**MMV scheme for active user detection:**]{} Reminiscing on Section \[sec:numeric\_synth\] and our previous work with Krzakala on MMV [@ZhuBaronKrzakala2017IEEE], more measurement vectors (larger $J$) lead to better estimation quality. We propose to convert the SMV style of the AUD problem into an MMV style by having each EU send $J$ different identification codewords, $\a_n^{(j)}\in \{-1,+1\}^{M\times 1}, \forall j\in \{1,\ldots,J\}$, to the base station. However, since the underlying signal $\x$  is Bernoulli and does not change during the AUD period, the resulting MMV is equivalent to an SMV with $J$ times more measurements; the column $\a_n$ of the measurement matrix in the equivalent SMV scheme is $\a_n = \l[\l(\a_n^{(1)}\r)^T,\cdots, \l(\a_n^{(J)}\r)^T\r]^T$. Hence, Algorithm \[algo:metric\_opt\_MMV\] should yield the same detection accuracy for the MMV scheme with $J$ channels and the SMV scheme with $J$ times larger measurement rate. Nevertheless, there is one advantage in adopting the MMV scheme: Lines \[line:beginForLoop\]–\[line:endForLoop\] of Algorithm \[algo:AMP\_MMV\][^13] can be parallelized with $J$ processing units (for example using a general purpose graphics processing unit or multicore computing system). After parallelizing Algorithm \[algo:AMP\_MMV\], the base station can perform the detection procedure with less runtime.

$$\label{eq:Ez}
\mathbb{E}[w|k,y,\Theta]=\left\{
                \begin{array}{ll}
                \displaystyle \frac{1}{\widetilde{Z}}\int  \frac{w}{1+\operatorname{e}^{-aw}} \frac{\operatorname{e}^{-\frac{1}{2\Theta}(w-k)^2}}{\sqrt{2\pi\Theta}}  dz=\frac{\sum_{u=1}^{u_{\text{max}}} \alpha_u T_1(u)}{\sum_{u=1}^{u_{\text{max}}} \alpha_u T_0(u)}=k+\frac{\Theta\sum_{u=1}^{u_{\text{max}}} \frac{\alpha_u \phi(\eta_u)}{\sqrt{\l(\frac{\sigma_u}{a}\r)^2+\Theta}}}{\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u)},\ &y=1,\\
               \displaystyle  \frac{1}{\widetilde{Z}}\int \frac{w\operatorname{e}^{-aw}}{1+\operatorname{e}^{-aw}} \frac{\operatorname{e}^{-\frac{1}{2\Theta}(w-k)^2}}{\sqrt{2\pi\Theta}}  dz=\frac{k-\sum_{u=1}^{u_{\text{max}}} \alpha_u T_1(u)}{1-\sum_{u=1}^{u_{\text{max}}} \alpha_u T_0(u)}=k-\frac{\Theta\sum_{u=1}^{u_{\text{max}}} \frac{\alpha_u \phi(\eta_u)}{\sqrt{\l(\frac{\sigma_u}{a}\r)^2+\Theta}}}{1-\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u)},\ &y=0.
                \end{array}
            \right.\\$$

$$\label{eq:mmse_p1}
\begin{split}
\mathbb{E}[(\mathbb{E}[\x_n|\q_n])^2]&=\int_{\q_n} f(\q_n) (\mathbb{E}[\x_n|\q_n])^2 d\q_n\\
&=\frac{\l(\frac{\rho}{1+\Delta_v}\r)^2}{\l[2\pi (1+\Delta_v)\r]^{J/2}}\int_{\q_n} \frac{\q_n \q_n^T}{\rho \exp\l[\frac{\q_n\q_n^T}{2(1+\Delta_v)}\r]+(1-\rho)\l(1+\frac{1}{\Delta_v}\r)^{J/2} \exp\l[\frac{\q_n\q_n^T (\Delta_v-1)}{2\Delta_v(1+\Delta_v)}\r]}  d\q_n.
\end{split}$$

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

In this paper, we studied the MMV signal estimation problem with user-defined additive error metrics on the estimate. We proposed an algorithmic framework that is optimal under arbitrary additive error metrics. We showed the optimality of our algorithm under certain conditions for SMV and conjectured its optimality for MMV. As examples, we derived algorithms that yield the optimal estimates in the sense of mean weighted support set error and mean absolute error, respectively. Numerical results not only verified the theoretic performance but also verified the intuition that having more signal vectors in MMV problems is beneficial to the estimation algorithm. We further provided simulation results for active user detection problem in multi-user communication systems, which is a real-world application of MMV models with the goal of minimizing the Hamming distance. Simulation results demonstrated the promise of our algorithm.

Derivation of $g_{\text{out}}$ for logistic channels  {#app:logit}
-----------------------------------------------------

Byrne and Schniter [@ByrneSchniter2015ArXiv] provide a method to derive $g_{\text{out}}$ for logistic channels , but the actual formula for $g_{\text{out}}$ is not given in their paper. To make our paper self-contained, we outline the derivation of $g_{\text{out}}$ for logistic channels. In order to calculate $g_{out}(k,y,\Theta)$ , we need to find $f(w|k,y,\Theta)$  and calculate $\mathbb{E}[w|k,y,\Theta]$. For logistic channels , $$\label{eq:logit_cond}
\begin{split}
& f(w|k,y,\Theta)=\\
&\frac{1}{\widetilde{Z}}\times \l[\frac{\delta(y-1)}{1+\operatorname{e}^{-w}}+\delta(y)\frac{\operatorname{e}^{-w}}{1+\operatorname{e}^{-w}}\r]\frac{1}{\sqrt{2\pi\Theta}} \operatorname{e}^{-\frac{1}{2\Theta}(w-k)^2},
\end{split}$$ where $\widetilde{Z}$ is a normalization factor. Therefore, it is difficult to calculate $\mathbb{E}[w|k,y,\Theta]$. Instead of calculating $\mathbb{E}[w|k,y,\Theta]$ by brute force, Byrne and Schniter [@ByrneSchniter2015ArXiv] use a mixture of Guassian cumulative distribution functions (CDF’s) to approximate the sigmoid function $\frac{1}{1+\text{exp}(-aw)} \approx\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\frac{w}{\sigma_u/a})$ [@Stefanski1991], where $u_{\text{max}}$ is the maximum number of Gaussian CDF’s one wants to use, $\Phi(\frac{w}{\sigma_u/a})$ denotes the Gaussian CDF whose standard deviation is $\frac{\sigma_u}{a}$, and $\alpha_u$ is the weight.

Following Byrne and Schniter [@ByrneSchniter2015ArXiv], we define the $i$-th moment $$\int w^i \mathcal{N}(w;k,\Theta) \Phi\l(\frac{w}{\sigma_u/a}\r)dz=T_i(u),$$ where $\mathcal{N}(w;k,\Theta)$ is the pdf of an RV $w$ with mean $k$ and variance $\Theta$, and $\Phi(\cdot)$ is the CDF of a standard Gaussian RV. Defining $\eta_u=\frac{k}{\sqrt{\l(\frac{\sigma_u}{a}\r)^2+\Theta}}$, we obtain $$\begin{aligned}
T_0(u)&=&\Phi(\eta_u),\nonumber\\
T_1(u)&= &k\Phi(\eta_u)+\frac{\Theta\phi(\eta_u)}{\sqrt{\l(\frac{\sigma_u}{a}\r)^2+\Theta}},\nonumber\end{aligned}$$ $$\begin{aligned}
T_2(u)=\frac{(T_1(u))^2}{\Phi(\eta_u)}+\Theta \Phi(\eta_u)-\frac{\Theta^2\phi(\eta_u)}{\l(\frac{\sigma_u}{a}\r)^2+\Theta}\l(\eta_u+\frac{\phi(\eta_u)}{\Phi(\eta_u)}\r),\nonumber\end{aligned}$$ where $\phi(\eta_u)$ is the pdf of a standard Gaussian RV at $\eta_u$. Hence, the normalization factor $\widetilde{Z}$ in  can be derived, $$\begin{split}
\widetilde{Z}=& \int \l[\frac{\delta(y-1)}{1+\operatorname{e}^{-aw}}+\delta(y)\frac{\operatorname{e}^{-aw}}{1+\operatorname{e}^{-aw}}\r]\frac{\operatorname{e}^{-\frac{1}{2\Theta}(w-k)^2}}{\sqrt{2\pi\Theta}}  dz\\
 = &\left\{
                \begin{array}{ll}
                 \sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u),\ &y=1,\\
                 1-\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u),\ &y=0.
                \end{array}
            \right.\\
\end{split}$$ We can further obtain the expression for $\mathbb{E}[w|k,y,\Theta]$ in . Hence, we can calculate $g_{\text{out}}$ .

Apart from $g_{\text{out}}$, we also need to find the partial derivative of $g_{\text{out}}$ , which according to Rangan [@RanganGAMP2011ISIT] satisfies $$\label{eq:deri_g_out}
-\frac{\partial}{\partial k}g_{\text{out}}(k,y,\Theta)=\frac{1}{\Theta}\l(1-\frac{\text{var}(w|k,y,\Theta)}{\Theta}\r),$$ where $\text{var}(w|k,y,\Theta)=\mathbb{E}[w^2|k,y,\Theta]-[\mathbb{E}[w|k,y,\Theta]]^2$. Note that $\mathbb{E}[w^2|k,y,\Theta]$ can be derived in the same way as  and the result is given below, $$\mathbb{E}[w^2|k,y,\Theta]=\left\{
                \begin{array}{ll}
                  \displaystyle k^2+\Theta+\frac{ \displaystyle\sum_{u=1}^{u_{\text{max}}} \alpha_u \xi_u}{\displaystyle\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u)},\ &y=1,\\
                 \displaystyle k^2+\Theta-\frac{\displaystyle\sum_{u=1}^{u_{\text{max}}} \alpha_u \xi_u}{\displaystyle 1-\sum_{u=1}^{u_{\text{max}}} \alpha_u \Phi(\eta_u)},\ &y=0,
                \end{array}
            \right.\\$$ where $$\xi_u=\frac{2k\Theta\phi(\eta_u)}{\sqrt{\l(\frac{\sigma_u}{a}\r)^2+\Theta}}-\frac{\Theta^2 \eta_u\phi(\eta_u)}{\l(\frac{\sigma_u}{a}\r)^2+\Theta}.$$

Inverting the MMSE {#app:inverMMSE}
------------------

For the MMV problem with i.i.d. matrices and joint Bernoulli-Gaussian signals, Zhu et al. provide an information-theoretic characterization of the MMSE by using replica analysis [@ZhuBaronKrzakala2017IEEE]. Suppose that we have already obtained the MMSE for an MMV problem. This appendix briefly shows how to invert the MMSE expression in order to obtain the equivalent scalar channel noise variance $\Delta_v$.

The optimal denoiser for the pseudo data is $\widetilde{\x}_n=\mathbb{E}[\x_n|\q_n]=f_{a_n}(\Delta_v,\q_n)$, where $f_{a_n}(\Delta_v,\q_n)$ is given in . We then express MMSE expression using $\mathbb{E}[\x_n|\q_n]$ as follows, $$\text{MMSE}=\mathbb{E}[(\widetilde{\x}_n-\x_n)^2]=J\rho-\mathbb{E}[(\mathbb{E}[\x_n|\q_n])^2].\label{eq:mmse1}$$ We calculate $\mathbb{E}[(\mathbb{E}[\x_n|\q_n])^2]$ in , where the $J$-dimensional integral can be simplified by a change of coordinates. Then, we plug  into , and express the MMSE as a function of $\Delta_v$. Finally, we numerically solve $\Delta_v$ for any given MMSE.

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

We would like to thank Jin Tan for providing valuable insights into achieving metric-optimal performance in signal estimation. Jong Chul Ye, Yavuz Yapici, and Ismail Guvenc helped us identify some real-world applications for minimizing error metrics that are different from the MSE. Finally, we are very grateful to the reviewers and Associate Editor Prof. Lops. In addition to their excellent suggestions, they were unusually flexible with us during the review process.

[10]{}

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: [R]{}ecovering 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.

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.

A. K. Fletcher, S. Rangan, and V. K. Goyal, “On-off random access channels: [A]{} compressed sensing framework,” , Mar. 2009.

V. Boljanovi[ć]{}, D. Vukobratovi[ć]{}, P. Popovski, and C. Stefanovi[ć]{}, “User activity detection in massive random access: Compressed sensing vs. coded slotted [ALOHA]{},” , June 2017.

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.

J. Zhu, D. Baron, and F. Krzakala, “Performance limits for noisy multimeasurement vector problems,” , vol. 65, no. 9, pp. 2444–2454, May 2017.

J. Zhu, , Ph.D. thesis, North Carolina State University, Raleigh, NC, Jan. 2017.

J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. [P]{}art [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. [P]{}art [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 Asilomar Conf. Signals, Syst., and Comput.*]{}, Nov. 2011, pp. 1447–1451.

J. Tan, D. Carmon, and D. Baron, “Optimal estimation with arbitrary error metric in compressed sensing,” in [*Proc. IEEE Stat. Signal Process. Workshop (SSP)*]{}, Aug. 2012, pp. 588–591.

J. Tan, D. Carmon, and D. Baron, “Signal estimation with additive error metrics in compressed sensing,” , vol. 60, no. 1, pp. 150–158, Jan. 2014.

S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” in [*Proc. IEEE Int. Symp. Inf. Theory (ISIT)*]{}, St. Petersburg, Russia, July 2011, pp. 2168–2172.

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.

S. Rangan, “Estimation with random linear mixing, belief propagation and compressed sensing,” in [*Proc. IEEE Conf. Inf. Sci. Syst. (CISS)*]{}, Princeton, NJ, Mar. 2010.

D. Baron, M. B. Wakin, M. F. Duarte, S. Sarvotham, and R. G. Baraniuk, “Distributed compressed sensing,” Tech. [R]{}ep. ECE-0612, Rice University, Dec. 2006.

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

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

J. Barbier and F. Krzakala, “Approximate message-passing decoder and capacity-achieving sparse superposition codes,” , vol. 63, no. 8, pp. 4894–4927, Aug. 2017.

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

E. Byrne and P. Schniter, “Sparse multinomial logistic regression via approximate message passing,” , Sept. 2015.

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

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.

C. Rush and R. Venkataramanan, “Finite sample analysis of approximate message passing algorithms,” , June 2016.

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

G. Hannak, A. Perelli, N. Goertz, G. Matz, and M. E. Davies, “Performance analysis of approximate message passing for distributed compressed sensing,” , Dec. 2017.

T. M. Cover and J. A. Thomas, , New York, NY, USA: Wiley-Interscience, 2006.

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

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.

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

Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “[Orthogonal matching pursuit: [R]{}ecursive function approximation with applications to wavelet decomposition]{},” , pp. 40–44, Nov. 1993.

L. Stefanski, “A normal scale mixture representation of the logistic distribution,” , vol. 11, no. 1, pp. 67–70, 1991.

[^1]: The work was supported in part by the National Science Foundation under the Grants CCF-1217749 and ECCS-1611112.

[^2]: Junan Zhu is with Bloomberg L.P., New York, NY 10017, and Dror Baron is with the Department of Electrical and Computer Engineering, NC State University, Raleigh, NC 27695. E-mail: {jzhu9, barondror}@ncsu.edu.

[^3]: In SMV, our algorithm closely resembles the one proposed by Tan and coauthors [@Tan2012SSP; @Tan2014], with the difference being Tan and coauthors rely on relaxed belief propagation (relaxed BP, an MP algorithm) [@Rangan2010CISS] and do not provide rigorous proofs for the optimality of their algorithm.

[^4]: Our proof is based on the approximate message passing algorithm [@DMM2009], while Tan and coauthors [@Tan2012SSP; @Tan2014] rely on relaxed BP and do not provide rigorous proofs.

[^5]: Byrne and Schniter [@ByrneSchniter2015ArXiv] describe without detail how to derive $g_{\text{out}}(\cdot)$  for i.i.d. parallel logistic channels ; we present a detailed derivation for completeness in Appendix \[app:logit\], and do not claim it as a contribution.

[^6]: Free energy is a term brought from statistical physics and is used to describe the interaction between the signals and the measurement matrices in linear models [@ZhuBaronCISS2013; @ZhuBaronKrzakala2017IEEE; @Krzakala2012probabilistic]. When the free energy has two fixed points, (G)AMP is not optimal with i.i.d. Gaussian matrices, and neither is Algorithm \[algo:metric\_opt\_MMV\]. We refer interested readers to the literature for detailed discussions [@ZhuBaronCISS2013; @ZhuBaronKrzakala2017IEEE; @Krzakala2012probabilistic].

[^7]: Pseudo-Lipschitz is a concept discussed in Bayati and Montanari [@Bayati2011]: For $k\geq 1$, we say a function $\phi: \mathbb{R}^m \rightarrow \mathbb{R}$ is [*pseudo-Lipschitz*]{} of order $k$ if there exists a constant $L>0$ such that for all $\x,\y\in \mathbb{R}^m$: $|\phi(\x)-\phi(\y)| \leq L(1+\|\x\|^{k-1}+\|\y\|^{k-1})\|\x-\y\|$.

[^8]: Our work with Krzakala [@ZhuBaronKrzakala2017IEEE] focuses on a diagonal covariance matrix for $\x_n$ in . A recent work by Hannak et al. [@Hannak2017] extends our work [@ZhuBaronKrzakala2017IEEE] to non-diagonal covariance matrices for $\x_n$. Following Hannak et al. [@Hannak2017], we can extend the performance limits analysis in this paper to non-diagonal covariance matrices for $\x_n$.

[^9]: Algorithm \[algo:metric\_opt\_MMV\] also applies to problems with other types of matrices, as long as the entries in the measurement matrices scale with $\frac{1}{\sqrt{N}}$. However, when the matrices are not i.i.d., there is no easy way to [*theoretically*]{} characterize the MMSE, the equivalent scalar channel noise variance $\Delta_v$, and the metric-optimal error. Such a theoretic characterization is sometimes necessary, because the MMSE behaves differently under different noise variances $\Delta_z$  and measurement rates $R$ [@ZhuBaronCISS2013; @ZhuBaronKrzakala2017IEEE].

[^10]: We do not have a replica analysis for logistic channels. In order to compute the theoretic MMAE, we use the average $\Delta_v$ from all the 50 simulations for each setting and calculate the MMAE with .

[^11]: Details can be found in Krzakala et al. [@Krzakala2012probabilistic] and Barbier and Krzakala [@BarbierKrzakala2017IT].

[^12]: Note that the entries of the signal estimated by OMP are not exactly 0’s and 1’s. Hence, in order to provide meaningful results, we threshold the OMP estimates before calculating the Hamming distance.

[^13]: Recall that Algorithm \[algo:metric\_opt\_MMV\] runs Algorithm \[algo:AMP\_MMV\] as a subroutine.
