Bayesian Approach with Extended Support Estimation for Sparse Regression

TL;DR

This paper introduces BMMP, a Bayesian greedy algorithm that improves sparse signal support estimation by considering multiple candidates and extended support sets, outperforming existing methods.

Contribution

The paper proposes a novel Bayesian greedy algorithm, BMMP, which uses extended support estimation and multiple candidates to enhance sparse regression support detection.

Findings

BMMP outperforms state-of-the-art methods in simulations.

BMMP approaches the ideal sparsity limit.

Extended support estimation improves support detection accuracy.

Abstract

A greedy algorithm called Bayesian multiple matching pursuit (BMMP) is proposed to estimate a sparse signal vector and its support given $m$ linear measurements. Unlike the maximum a posteriori (MAP) support detection, which was proposed by Lee to estimate the support by selecting an index with the maximum likelihood ratio of the correlation given by a normalized version of the orthogonal matching pursuit (OMP), the proposed method uses the correlation given by the matching pursuit proposed by Davies and Eldar. BMMP exploits the diversity gain to estimate the support by considering multiple support candidates, each of which is obtained by iteratively selecting an index set with a size different for each candidate. In particular, BMMP considers an extended support estimate whose maximal size is $m$ in the process to obtain each of the support candidates. It is observed that BMMP outperforms other state-of-the-art methods and approaches the ideal limit of the signal sparsity in our simulation setting.

Equations

$$y=\Phi x^{\ast }+w,$$

$$\mathbb{E}[{\Vert r(\Delta )\Vert }^2]$$

$$\stackrel{(b)}{=}(|\Omega \setminus \Delta |\cdot {\sigma }^2\cdot v_x^2+{\sigma }_w^2)\cdot (m-d),$$

$$\Psi (\Delta ):=\frac{\max ({\Vert r(\Delta )\Vert }^2/(m-d)-{\sigma }_w^2,0)}{{\sigma }^2\cdot v_x^2}.$$

$$g_i(\Delta )$$

$$=\Vert {\bar{\varphi }}_i(\Delta )\Vert x_i+\sum _{l\in \Omega \setminus (\Delta \cup \{i\})}{\dot{\varphi }}_i(\Delta {)}^{\top }{\bar{\varphi }}_l(\Delta )x_l+{\dot{\varphi }}_i(\Delta {)}^{\top }\bar{w}(\Delta ),$$

$${\mathcal{T}}_0:g_i(\Delta )$$

$${\mathcal{T}}_1:g_i(\Delta )$$

$$\mathbb{E}[z_i|x_i=0]$$

$$=\sum _{l\in \Omega \setminus \Delta }\mathbb{E}[({\dot{\varphi }}_i(\Delta {)}^{\top }{\varphi }_l)]m_x+\mathbb{E}[({\dot{\varphi }}_i(\Delta {)}^{\top }w)]\stackrel{(a)}{=}0$$

$$\mathbb{E}[z_i^2|x_i=0]$$

$$=\sum _{l\in \Omega \setminus \Delta }\mathbb{E}[({\dot{\varphi }}_i(\Delta {)}^{\top }{\varphi }_l{)}^2]v_x^2+\mathbb{E}[({\dot{\varphi }}_i(\Delta {)}^{\top }w{)}^2]$$

$$\stackrel{(a)}{=}|\Omega \setminus \Delta |\cdot {\sigma }^2{v}_x^2+{\sigma }_w^2$$

$$\stackrel{(b)}{\simeq }\Psi (\Delta )\cdot {\sigma }^2{v}_x^2+{\sigma }_w^2=:{\tilde{\sigma }}_0^2,$$

$$\mathbb{E}[z_i|x_i=u]=\mathbb{E}[{\Vert {\bar{\varphi }}_i(\Delta )\Vert }_{2}u]=\sigma \tau u$$

$$\mathbb{E}[(z_i-\mathbb{E}[z_i]{)}^2|x_i=u]$$

$$\simeq (\Psi (\Delta )-1)\cdot {\sigma }^2{v}_x^2+{\sigma }_w^2=:{\tilde{\sigma }}_1^2,$$

$$\Theta (z_i)$$

$$=\ln (\frac{{\int }_{-\infty }^{\infty }\mathbb{P}(z_i|x_i=u)f_{x^{\ast }}(u)du}{\mathbb{P}(z_i|x_i=0)})$$

$$\stackrel{(a)}{\simeq }\ln (\frac{{\int }_{-\infty }^{\infty }\frac{1}{{\tilde{\sigma }}_1\sqrt{2\pi }}\exp (-\frac{|z_i-\sigma \tau u{|}^2}{2{\tilde{\sigma }}_1^2})f_{x^{\ast }}(u)du}{\frac{1}{{\tilde{\sigma }}_0\sqrt{2\pi }}\exp (-\frac{|z_i{|}^2}{2{\tilde{\sigma }}_0^2})}),$$

$$\frac{(z_i{)}^2}{2{\tilde{\sigma }}_0^2}+\ln (\text{erf}(\frac{\sigma \tau b-z_i}{{\tilde{\sigma }}_1\sqrt{2}})-\text{erf}(\frac{\sigma \tau a-z_i}{{\tilde{\sigma }}_1\sqrt{2}})),$$

$$(\gamma ,\eta )=\underset{\bar{\gamma }\in {\mathbb{R}}^a,\bar{\eta }\in \mathbb{R}}{\mathrm{argmin}}L(\bar{\Phi },y):=\log |\Sigma |+y^{\top }{\Sigma }^{-1}y,$$

$$\Vert P_{\mathcal{R}(A)}^{\perp }b\Vert =\sqrt{\sum _{j=1}^{m-d}(b^{\top }{e}_j{)}^2},$$

$$E:={\cup }_{\Gamma \subseteq \{1:n\}\text{ s.t. }|\Gamma \cap \Omega |<|\Omega |\text{ and }|\Gamma |<m}\mathcal{R}({\Phi }_{\Gamma }).$$

Full text

I Introduction

Sparse linear regression, refered to as compressed sensing (CS) [1, 2], has been widely studied in many fields ranging from denoising [3] to super-resolution restoration [4, 5]. A signal vector $x^{*}\in\mathbb{R}^{n}$ is called $k$-sparse if it has at most $k$ non-zero elements. Then, the objective of sparse linear regression is to recover a $k$-sparse signal vector $x^{*}\in\mathbb{R}^{n}$ from the given measurement vector $y\in\mathbb{R}^{m}$, such that

[TABLE]

where $\Phi\in\mathbb{R}^{m\times n}$ is a known sensing matrix with $m<n$ and a noise vector $w\in\mathbb{R}^{m}$.

To improve recovery performance of the sparse signal $x^{*}$ and its support $\Omega$, the set of non-zero indices in $x^{*}$, a method called MAP support detection [6] to recover a partial support of $x^{*}$ was proposed by exploiting the maximum a posteriori (MAP) estimation. This method determines whether the index $i$ belongs to $\Omega$ by using a likelihood ratio test under the true and null hypotheses given a correlation between the residual vector and its $l_{2}$-normalized column vector $\phi_{i}/\left\|\phi_{i}\right\|$ in $\Phi:=[\phi_{1},...,\phi_{n}]$. By applying this method to existing algorithms such as generalized OMP (gOMP) [7], compressive sampling matching pursuit (CoSaMP) [8], and subspace pursuit (SP) [9], some greedy methods were proposed and shown to outperform other existing methods[6].

We propose a new MAP-based method to estimate partial support by modifying the MAP support detection to use the correlation term given in rank aware order recursive matching pursuit (RA-ORMP) [10]. This proposed scheme improves the existing MAP support detection by additionally utilizing an orthogonal complement onto each column vector in $\Phi$ to calculate the correlation, since the correlation can be better represented through this consideration.

Based on this method to estimate a partial support, we develop an algorithm called Bayesian multiple matching pursuit (BMMP) to recover $x^{*}$ and its support $\Omega$. BMMP also uses the following four techniques to enhance the performance further. Simulation results show that BMMP outperforms other state-of-the-art methods in both noiseless and noisy cases.

(Generating an extended support estimate of size bounded by $m$) The extended support is defined by an index set whose size is larger than the actual sparsity $k$ of $x^{*}$ and includes $\Omega$. gOMP [7], SP [9], and CoSaMP [8] exploited an extended support estimate in the process to estimate $\Omega$. To improve the support recovery performance, the maximum size of the extended support estimate in BMMP is set to $m$ while the maximum size of the extended support estimate in gOMP, SP, and CoSaMP, is $t\cdot\min(k,\left\lfloor{m/t}\right\rfloor)$ , $2k$, and $3k$, respectively, for a constant $t$ smaller than $k$. The reason for setting this to $m$ in BMMP is explained in the section describing BMMP. 2. 2.

(Generating multiple support candidates) BMMP utilizes the diversity gain by iteratively selecting multiple indices of different size to generate each of multiple support candidates. In other words, BMMP generates the $i$th support candidate by iteratively selecting $i$ indices through the proposed partial support estimation. There have been related studies generating multiple support candidates through a tree structure: multipath matching pursuit [11], multi-branch matching pursuit [12]. However, the complexity of the tree-based methods scales in the sparsity $k$ and is higher than those of other CS algorithms, since the number of support candidates increases exponentially with $k$ (the depth of the tree) due to the structural nature of the tree. In contrast with these methods, BMMP does not use the tree-based approach to reduce complexity so that the number $g$ of support candidates is independent of $k$ and the complexity is linearly increased with $g$ in BMMP. 3. 3.

(Updating extended support estimate by replacing its subset) To improve the support recovery performance, SP [9] and CoSaMP [8] update an extended support estimate $\Delta$ by selecting $k$ indices and replacing the rest of the selected indices in $\Delta$ at each iteration. BMMP modifies this technique by iteratively selecting $k/2$ indices in a given extended support estimate $\Delta$ and updating its complement in $\Delta$. 4. 4.

A ridge regression proposed by Wipf and Rao[13] is used to improve robustness to noise.

II Notation

The set $\{1,2,...,i\}$ is denoted by $\{1:i\}$. The submatrix of a matrix $A:=[a_{1},...,a_{n}]$ $\in\mathbb{R}^{m\times n}$, where $a_{i}$ is its $i$th column, with columns indexed by $J\subseteq\{1:n\}$ is denoted by $A_{J}$. $A^{Q}$ denotes the submatrix of $A$ with rows indexed by $Q\subseteq[m]$. $\mathcal{R}(A)$ denotes the range space spanned by the columns of $A$. $A^{\top}$ denotes the transpose of $A$. $P_{\mathcal{R}(A)}:=A(A^{\top}A)^{-1}A^{\top}$ and $P^{\perp}_{\mathcal{R}(A)}$ denote the orthogonal projection onto $\mathcal{R}(A)$ and its orthogonal complement, respectively. $\left\|\cdot\right\|$ denotes the Frobenius norm. For a set $\Delta\subseteq\{1:n\}$ and a space $\mathcal{R}(A_{\Delta})$ of $\mathbb{R}^{d}$, $\dot{a}_{i}({\Delta}):={P^{\perp}_{\mathcal{R}(A_{\Delta})}a_{i}}/{\left\|P^{\perp}_{\mathcal{R}(A_{\Delta})}a_{i}\right\|}$, and $\dot{a}_{i}:={a_{i}}/{\left\|a_{i}\right\|}$. Similarly, $\bar{a}_{i}({\Delta}):=P^{\perp}_{\mathcal{R}(A_{\Delta})}a_{i}$.

III Partial support estimation by MAP

The proposed MAP-based algorithm for estimating a partial true support is introduced in this section. We assume the following. $\Omega$ is uniformly distributed on $\{1:n\}$. Each nonzero element $u$ of $x^{*}$ is i.i.d and follows an arbitrary distribution $f_{x^{*}}(u)$ whose mean and variance are $m_{x}$ and $\sigma_{x}^{2}$, respectively. Each element in $\Phi$ and $w$ follows the Gaussian distribution whose mean is [math] and variance is $\sigma^{2}$ and $\sigma_{w}^{2}$, respectively.

Suppose that $\Delta\subseteq\{1:n\}$ such that $|\Delta|=:d<m$ is a given partial support estimate. The goal of the proposed partial support detection is to find an index set belonging to $\Omega\setminus\Delta$ by using the inner product used in RA-ORMP, $g_{i}(\Delta):=r(\Delta)^{\top}\dot{\phi}_{i}(\Delta)$ for $i\in\{1:n\}\setminus\Delta$, where $r(\Delta):=P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y$ is the residual vector, and $\dot{\phi}_{i}(\Delta):={P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}\phi_{i}/\left\|{P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}\phi_{i}\right\|$. $g_{i}(\Delta)$ represents the correlation between the residual vector and $i$.

Although MAP support detection [6] estimates a partial true support from the likelihood ratio of the inner product $h_{i}(\Delta):=r(\Delta)^{\top}\dot{\phi}_{i}$, where $\dot{\phi}_{i}:=\phi_{i}/\left\|\phi_{i}\right\|$, the proposed method uses the likelihood ratio of $g_{i}(\Delta)$ for $i\in\{1:n\}\setminus\Delta$. Note that $g_{i}(\Delta)$ and $h_{i}(\Delta)$ use $\dot{\phi}_{i}(\Delta)$ and $\dot{\phi}_{i}$, respectively. $h_{i}(\Delta)$ can be interpreted as the correlation term in a normalized version of OMP since the normalized vector $\dot{\phi}_{i}$ is used in $h_{i}(\Delta)$, whereas OMP uses $\phi_{i}$ in its correlation term $({P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y)^{\top}\phi_{i}$. This indicates that the proposed method additionally considers the orthogonal complement ${P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}$ of $\Phi_{\Delta}$ for the column vector $\phi_{i}$ to represent the correlation compared to MAP support detection. Since each vector space $\mathcal{R}(\dot{\phi}_{i}(\Delta))$ for $i\in\Omega\setminus\Delta$ is guaranteed to belong to a space $\cup_{j\in\Omega\setminus\Delta}\mathcal{R}({P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}\phi_{j})$, which is extended from the residual space $\mathcal{R}({P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y)$, the index selection based on $g_{i}(\Delta)$ may more successfully find the index in $\Omega\setminus\Delta$ compared to that based on $h_{i}(\Delta)$.111It is not guaranteed that each vector space $\mathcal{R}(h_{i}(\Delta))$ for $i\in\Omega\setminus\Delta$ belongs to the extended residual space $\cup_{j\in\Omega\setminus\Delta}\mathcal{R}({P}^{\perp}_{\mathcal{R}(\Phi_{\Delta})}\phi_{j})$.

The proposed partial support detection uses an estimate of $|\Omega\setminus\Delta|$. This is obtained by the following derivations.

[TABLE]

where $v_{x}:=\sqrt{(m^{2}_{x}+\sigma_{x}^{2})}$, (a) follows from the independence assumption, and (b) follows from Lemma A.1 (with $A=\Phi_{\Delta}$ and $b=\phi_{l}$ or $w$) in that $\mathbb{E}[\left\|\bar{\phi}_{l}(\Delta)\right\|^{2}]=\sigma^{2}\cdot(m-d)$ and $\mathbb{E}[\left\|\bar{w}(\Delta)\right\|^{2}]=\sigma_{w}^{2}\cdot(m-d)$. Then, the estimate $\Psi(\Delta)$ for $|\Omega\setminus\Delta|$ is given as

[TABLE]

Note that the following equalities hold for $i\in\{1:n\}\setminus\Delta$

[TABLE]

where $x_{i}$ is the $i$th element of $x^{*}$ and (a) follows from $r(\Delta)=\sum\limits_{l\in\Omega\setminus\Delta}\bar{\phi}_{l}(\Delta)x_{l}+\bar{w}(\Delta)$. Our derivation in the rest of this section is based on [6]. From the equalities in (2), $g_{i}(\Delta)$ becomes the right-side terms in (3) and (4) under the null and true hypotheses—$\mathcal{T}_{0}$ and $\mathcal{T}_{1}$, respectively. For $\mathcal{T}_{0}$, $i$ does not belong to $\Omega$, and $x_{i}=0$. For $\mathcal{T}_{1}$, $i$ belongs to $\Omega$, and $x_{i}=u$.

[TABLE]

From (2) and (3), the expectation and variance of $z_{i}:=g_{i}(\Delta)$ under the null hypothesis can be respectively approximated as follows:

[TABLE]

[TABLE]

where (a) and (b) follow from the circular symmetry of the Gaussian distribution and (1), respectively. Similarly, by (4) and Lemma A.1 (with $A=\Phi_{\Delta}$ and $b=\phi_{i}$), the expectation and variance of $z_{i}$ under the true hypothesis are obtained respectively as

[TABLE]

and

[TABLE]

where $\tau:={\sqrt{2}\cdot\Gamma(\frac{1+m-d}{2})}/{\Gamma(\frac{m-d}{2})}$, and $\Gamma(\cdot)$ is the Gamma function.

Then, using (5)–(8), we obtain the log-likelihood ratio $\Theta(z_{i}):=\ln({\mathbb{P}(i\in\Omega|z_{i})}/{\mathbb{P}(i\notin\Omega|z_{i})})$ as follows.

[TABLE]

where (a) follows from the Gaussian approximation of the conditional distributions of $z_{i}$ given $x_{i}=0$ ($\mathcal{T}_{0}$) by using (5) and (6) and given $x_{i}=u$ ($\mathcal{T}_{1}$) by using (7) and (8).

Therefore, the proposed method estimates a partial true support as an index set $\Lambda$ by selecting the $|\Lambda|$ largest ratios $\Theta(z_{i})$ for $i\in\{1:n\}\setminus\Delta$ in $(\ref{loglikelihood})$. Note that $\Theta(z_{i})$ is propotional to (10) when each nonzero element of $x^{*}$ follows the uniform distribution $f_{x^{*}}(u;a,b)=1/(b-a)$ for $a\leq u\leq b$

[TABLE]

where $\textup{erf}(x):=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$.

IV Algorithm description

The proposed method, BMMP (Algorithm 1), for estimating ${x^{*}}$ and $\Omega$ is introduced in this section. BMMP returns ($\hat{x},\hat{\Omega}$) as estimates of (${x^{*}},\Omega$) given the tuple ($y,\Phi,k,g,\epsilon$), where $g$ is the number of support candidates from which the final support is estimated, and $\epsilon$ is a threshold used to terminate BMMP when the residual error is smaller than this value. The estimate $\hat{\Omega}$ is obtained in Steps 23–24 by selecting one of $g$ support candidates ($\dot{\Omega}_{t}$ for $t\in\{1:g\}$) with the minimum residual error. The $t$th support candidate $\dot{\Omega}_{t}$ in Step 23 is obtained by the iteration in Steps $3$–$18$ consisting of the following four stages.

In the first stage (generate an extended support estimate $\Delta$) in Steps 4–10, a set $\Lambda$ of $t$ indices in Step 8 is iteratively obtained through the proposed partial support estimation222We note that $\Theta(z_{i})$ in (10) is calculated by using the following values: $(y,\Delta,\sigma,$ $m_{x},\sigma_{x},\sigma_{w})$. Step 6 in Algorithm 1 specifies these values as inputs of $\Theta(\cdot)$. and added to a set $\Delta$ in Step 9 until its size $|\Delta|$ is equal to $m$ or the residual error obtained from the estimate is smaller than $\epsilon$, i.e. $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|\leq\epsilon$. In the second stage (generate a temporary estimate $\bar{x}$ of ${x^{*}}$) in Step 11, an estimate $\bar{x}$ of the signal vector supported on $\Delta$ is obtained by the ridge regression described in Step 11. In the third stage in Steps 12–15, a temporary support estimate $\bar{\Omega}_{i+1}$ is obtained by selecting the $k$-largest elements in $\bar{x}^{\Delta}$. The fourth stage in Steps 16–17 generates a subset of $\Delta$ by selecting $k/2$ indices from $\Delta$. Then, the extended support estimate $\Delta$ is reinitialized as the subset in the next iteration step. This iteration terminates and the $t$th support candidate $\dot{\Omega}_{t}$ is obtained in Step 19 when the residual error obtained from the temporary support estimate $\bar{\Omega}_{i+1}$ in Step 15 does not decrease further in Step 3.

To enhance the robustness to noise, a ridge regression proposed by Wipf and Rao [13] is used in Step 11 with the following parameters $(\gamma,\eta)$ obtained by minimizing the cost $L(\cdot)$ in (11).

[TABLE]

where $\bar{\Phi}:=\Phi_{\Delta}\in\mathbb{R}^{m\times a}$ is a submatrix of $\Phi$ with columns indexed by $\Delta$, $\Sigma=(\bar{\eta}^{-2}{\bar{\Phi}}^{\top}{\bar{\Phi}}+D(\bar{\gamma}))^{-1}$, and $D(\bar{\gamma})$ is the diagonal matrix whose $i$th diagonal element is $\gamma_{i}$ in $\bar{\gamma}=(\gamma_{1},...,\gamma_{a})^{\top}$. We obtain an appoximated solution of $(\ref{sble})$ by using sparse bayesian learning (SBL) [13] to reduce the complexity.

In the noiseless case, the parameters $(\bar{\gamma},\bar{\eta})$ in the ridge regression are set to zero so that the ridge regression becomes the least-square regression ($\underset{x\in\mathbb{R}^{n}}{\arg\min}\left\|\Phi_{\Delta}x^{\Delta}-y\right\|$). The least-squares method has a lower complexity than SBL and it is well-known that the least-squares method provides a unique solution as ${x^{*}}$ in the noiseless case if $|\Delta|\leq m$, $\Delta\supseteq\Omega$, and $\Phi_{\Delta}$ has the full column rank. If $\Delta\supseteq\Omega$ holds, the inversion problem of CS can be simplified as an easier problem where $\Phi$ is replaced by its submatrix $\Phi_{\Delta}$ compared to the original problem. When $\Delta$ has a larger size, the probability of satisfying $\Delta\supseteq\Omega$ is greater. For this reason, we set the maximal size of the extended support estimate $\Delta$ as $m$. Besides, it is guaranteed from Lemma A.2 that in the noiseless case, $\Omega$ is a subset of $\Delta$ ($\Delta\supseteq\Omega$) if $|\Delta|<m$ and $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|$ is equal to zero. Based on this fact, to minimize the size of $\Delta$ satisfying $\Delta\supseteq\Omega$, we increase the size of $\Delta$ from [math] to $m$ and find $\Delta$ satisfying the criterion $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|\leq\epsilon$ in Steps 4–10. Figure 1 illustrates the procedure of BMMP.

V Numerical experiments

We compare the performance of BMMP and existing state-of-the-art algorithms such as MAP-gOMP333Two indices are selected at each step in MAP-gOMP., MAP-SP, MAP-CoSaMP, MAP-OMP [6], RegL1 ($l_{1}$ norm minimization with a nonnegative constraint).444RegL1 outputs the signal estimate $\hat{x}=(x_{1},...,x_{n})^{\top}$ by minimizing $\sum_{i}|x_{i}|$ such that $x_{i}\geq 0$ for $i\in\{1:n\}$ and $\left\|y-\Phi\hat{x}\right\|\leq\tau$ where $\tau$ is the noise magnitude. We note that MAP-SP, MAP-CoSaMP, MAP-OMP are algorithms based on MAP support detection [6]. Each entry of $\Phi$ is i.i.d. and follows the Gaussian distribution $\mathcal{N}(0,1/m)$ with mean zero and variance $1/m$. The elements of ${(x^{*})}^{\Omega}$ are independently and uniformly sampled from [math] to $1$. Then, $(\sigma,m_{x},\sigma_{x})$ is $(1/\sqrt{m},1/2,1/\sqrt{12})$.

The rate of successful support recovery, i.e., $\hat{\Omega}=\Omega$, and the execution time of each algorithms in the noiseless case are shown in Fig. 2(a) and Fig. 2(b), respectively. ($m,n$) is set to ($128,256$) in Fig. 2. We evaluate BMMP whose input $(g,\epsilon,\lambda)$ is set to $(4,\left\|y\right\|\cdot 10^{-\textup{SNR}_{\textup{dB}}/20},\epsilon^{2}/m)$, where SNR${}_{\textup{dB}}$ is the signal-to-noise ratio (SNR $:=\mathbb{E}\left\|\Phi_{\Omega}x^{\Omega}_{0}\right\|^{2}/\,\mathbb{E}\left\|w\right\|^{2}$) in decibels. It has been shown [14] that in the noiseless case, the maximal sparsity $\bar{k}$ required for an ideal approach to recover any signal vector ${x^{*}}$ such that $|\Omega|\leq\bar{k}$ is equal to $m/2$ given $m$, where $\Omega$ is the support of ${x^{*}}$; the $l_{0}$ bound is denoted by this ideal limit of the signal sparsity ($m/2$). Fig. 2(a) and Fig. 2(b) show that BMMP has a lower complexity than RegL1, outperforms other state-of-the-art methods, and approaches the $l_{0}$ bound, i.e., $m/2=64$.

Fig. 2(c) illustrates the mean squared error $\mathbb{E}\left\|\hat{x}-{x^{*}}\right\|^{2}$ of each algorithm when SNR is varied from $20$ to $40$ dB555In the noisy case, $x^{*}$ is randomly sampled such that each of its non-zero elements is ranged from $0.1$ to $1$., the sparsity $k$ is fixed to $60$. It is observed that the signal reconstruction performance of BMMP is better than those of other methods in the noisy case.

Fig. 2(d) shows how each of the following three techniques used in BMMP contributes to the performance improvement of BMMP: exploiting the extended support estimate with its maximal size equal to $m$ ($\mathcal{E}$), using multiple support estimates ($\mathcal{M}$), and iteratively updating the support estimate by replacing its subset ($\mathcal{U}$). BMMP$\setminus\{\mathcal{U}\}$ refers to BMMP without $\mathcal{U}$ by setting the value $\left\lfloor{k/2}\right\rfloor$ to [math] in step 16 in Algorithm 1. BMMP$\setminus\{\mathcal{U},\mathcal{M}\}$ denotes BMMP without $\{\mathcal{U},\mathcal{M}\}$ by setting $g$ to $1$ in BMMP$\setminus\{\mathcal{U}\}$. Similarly, BMMP$\setminus\{\mathcal{U},\mathcal{M},\mathcal{E}\}$ indicates BMMP without $\{\mathcal{U},\mathcal{M},\mathcal{E}\}$ such that $m$ (the maximum size of the extended support estimate $\Delta$) shown in steps 4-10 in Algorithm 1 is set to $k$ and the remaining settings follow BMMP$\setminus\{\mathcal{U},\mathcal{M}\}$. It is shown that the performance improves in the order of BMMP$\setminus\{\mathcal{U},\mathcal{M},\mathcal{E}\}$, BMMP$\setminus\{\mathcal{U},\mathcal{M}\}$, BMMP$\setminus\{\mathcal{U}\}$, and BMMP. Fig. 2(d) shows that the consideration of $\mathcal{E}$ contributes to a significant performance improvement of BMMP.

To compare performance of the proposed partial support estimation and MAP support detection [6], we compare performance of MAP-gOMP and its variant, which is obtained by replacing the MAP support detection with the proposed partial support estimation in MAP-gOMP. Fig. 3 shows the rate of successful support recovery of these two algorithms in the noiseless case where $m$ is varied from $32$ to $256$ with $k=\left\lfloor{m/1.8}\right\rfloor$ and $n=2m$. It is observed that the proposed partial support estimation outperforms the MAP support detection [6].

Fig. 4 shows the performance for reconstructing a grayscale sparse image with a size of $16\times 16$ pixels in the noisy case when SNR = $25$ dB. This image is compressed by using $\Phi\in\mathbb{R}^{138\times 256}$, whose elements are i.i.d. and sampled from the Gaussian distribution $\mathcal{N}(0,$ $1/138)$ with mean [math] and variance $1/138$. It is observed in Fig. 4 that BMMP recovers the image better than other methods.

VI Conclusion

We presented BMMP, which updates multiple extended support estimates with each size equal to $m$ and performs a likelihood ratio test given the correlation term in RA-ORMP. The numerical results show that BMMP achieves an improvement in performance compared to existing state-of-the-art methods, even with noisy measurements. Future related studies will consider the application of a technique for estimating the moments of $\Phi$, $w$, and $x^{*}$ to BMMP and the development of BMMP in the case where the distribution of $x^{*}$ is designable or the a priori information of $x^{*}$ is available, i.e., a communication system for nonorthogonal multiple access [15, 16].

Appendix A Lemmas

Lemma A.1**.**

Let $b\in\mathbb{R}^{m}$ be a vector such that its elements are i.i.d. and follow the Gaussian distribution $\mathcal{N}(0,\hat{\sigma}^{2})$, and let $A:=[a_{1},...,a_{d}]\in\mathbb{R}^{m\times d}$ be a matrix $A$ with rank $d<m$. Then, it is guaranteed that $\mathbb{E}[\left\|P^{\perp}_{\mathcal{R}(A)}b\right\|]=\hat{\sigma}\cdot\tau$ and $\mathbb{E}[\left\|P^{\perp}_{\mathcal{R}(A)}b\right\|^{2}]=\hat{\sigma}^{2}\cdot(m-d)$, where $\tau:={\sqrt{2}\cdot\Gamma(\frac{1+m-d}{2})}/{\Gamma(\frac{m-d}{2})}$, and $\Gamma(\cdot)$ is the Gamma function.

Proof.

Let $\{e_{1},...,e_{m-d}\}$ be an orthonormal basis of $\mathcal{R}^{\perp}(A)$. Then, $P^{\perp}_{\mathcal{R}(A)}b=\sum_{j=1}^{m-d}(b^{\top}e_{j})e_{j}$. Note that each $b^{\top}e_{j}$ for $j\in\{1:m-d\}$ follows $\mathcal{N}(0,\hat{\sigma}^{2})$ by the isotropic property of the Gaussian vector [17]. Since it follows that

[TABLE]

$\left\|P^{\perp}_{\mathcal{R}(A)}b\right\|/\hat{\sigma}$ follows the chi distribution with $m-d$ degrees of freedom. Then, the proof is completed by averaging the chi and chi-squared distributions. ∎

Lemma A.2**.**

Suppose that every $m$ columns in $\Phi$ exhibit full rank, $y$ is uniformly sampled from $\mathcal{R}(\Phi_{\Omega})$, and the sparsity $k$ is smaller than $m$, i.e., $|\Omega|<m$. Then, in the noiseless case, the following statements (a) and (b) hold for any set $\Delta\subseteq\{1:n\}$ such that $|\Delta|<m$.

1. (a)

$\Delta\nsupseteq\Omega$* holds if $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|>0$.* 2. (b)

$\Delta\supseteq\Omega$* holds almost surely if $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|=0$.*

Proof.

The following subspace $E$ is defined as

[TABLE]

For any set $\Delta\subseteq\{1:n\}$ such that $|\Delta|<m$, $y\notin\mathcal{R}(\Phi_{\Delta})$ holds if $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|>0$. If $\Delta\supseteq\Omega$, $\mathcal{R}(\Phi_{\Delta})$ includes $\mathcal{R}(\Phi_{\Omega})$ so that $y\in\mathcal{R}(\Phi_{\Delta})$. Thus, $y\notin\mathcal{R}(\Phi_{\Delta})$ implies that $\Delta\nsupseteq\Omega$ so that the statement (a) is satisfied.

Given that from the assumption for $\Phi$, the rank of $\mathcal{R}(\Phi_{\Omega})$ is strictly larger than that of $\mathcal{R}(\Phi_{\Gamma})\cap\mathcal{R}(\Phi_{\Omega})$ for any index set $\Gamma\subseteq\{1:n\}$ such that $|\Gamma\cap\Omega|<|\Omega|$ and $|\Gamma|<m$, the event region satisfying $y\in\mathcal{R}(\Phi_{\Omega})\cap E$ has Lebesgue measure zero on the range space $\mathcal{R}(\Phi_{\Omega})$. Thus, the condition $y\in\mathcal{R}(\Phi_{\Omega})\setminus E$ holds almost surely. Given that the condition $\left\|P^{\perp}_{\mathcal{R}(\Phi_{\Delta})}y\right\|=0$ implies that $y\in\mathcal{R}(\Phi_{\Delta})$, $y\in\mathcal{R}(\Phi_{\Delta})\cap(\mathcal{R}(\Phi_{\Omega})\setminus E)$ holds almost surely. And this condition $y\in\mathcal{R}(\Phi_{\Delta})\cap(\mathcal{R}(\Phi_{\Omega})\setminus E)$ implies that ${\Delta}\supseteq\Omega$ from the definition of $E$. Thus, the statement (b) is satisfied. ∎