# Unitary Approximate Message Passing for Sparse Bayesian Learning

Man Luo, Qinghua Guo, , Ming Jin, Yonina C. Eldar, , Defeng (David) Huang, ,
and Xiangming Meng
Part of this work was presented in the 16th IEEE APWCS 2019. Corresponding
author: Qinghua Guo.M. Luo and Q. Guo are with the School of Electrical
Computer and Telecommunications Engineering, University of Wollongong,
Australia, (e-mail<EMAIL_ADDRESS>qguo@uow.edu.au).M. Jin is with the
Faculty of Electrical Engineering and Computer Science, Ningbo University,
Ningbo, China (e-mail: jinming@nbu.edu.cn).Y. C. Eldar is with the Faculty of
Math and CS, Weizmann Institute of Science, Rehovot, Israel (email:
yonina.eldar@weizmann.ac.il).D. Huang is with the School of Engineering,
University of Western Australia, Perth Australia (e-mail:
david.huang@uwa.edu.au).X. Meng is with the Institute for Physics of
Intelligence, The University of Tokyo, Hongo, Tokyo 113-0033, Japan (email:
meng@g.ecc.u-tokyo.ac.jp).

###### Abstract

Sparse Bayesian learning (SBL) can be implemented with low complexity based on
the approximate message passing (AMP) algorithm. However, it does not work
well for a generic measurement matrix, which may cause AMP to diverge. Damped
AMP has been used for SBL to alleviate the problem at the cost of reducing
convergence speed. In this work, we propose a new SBL algorithm based on
structured variational inference, leveraging AMP with a unitary transformation
(UAMP). Both single measurement vector and multiple measurement vector
problems are investigated. It is shown that, compared to state-of-the-art AMP-
based SBL algorithms, the proposed UAMP-SBL is more robust and efficient,
leading to remarkably better performance.

###### Index Terms:

Sparse Bayesian learning, structured variational inference, approximate
message passing.

##  I Introduction

We consider the problem of recovering a sparse signal $\mathbf{x}$ from noisy
measurements $\mathbf{y}=\mathbf{Ax}+\mathbf{w}$, where $\mathbf{A}$ is a
known measurement matrix [1]. This problem finds numerous applications in
various areas of signal processing, statistics and computer science. One
approach to recovering $\mathbf{x}$ is to use sparse Bayesian learning (SBL),
where $\mathbf{x}$ is assumed to have a sparsity-promoting prior [2].
Conventional implementation of SBL involves matrix inversion in each
iteration, resulting in prohibitive computational complexity for large scale
problems.

The approximate message passing (AMP) algorithm [3], [4] has been proposed for
low-complexity implementation of SBL [5, 6]. AMP was originally developed for
compressive sensing based on loopy belief propagation (BP) [4]. Compared to
convex optimization based algorithms such as LASSO [7] and greedy algorithms
such as iterative hard-thresholding [8], AMP has low complexity and its
performance can be rigorously characterized by a scalar state evolution (SE)
in the case of a large independent and identically distributed (i.i.d.)
(sub-)Gaussian matrix $\mathbf{A}$ [9]. AMP was later extended in [10] to
solve general estimation problems with a generalized linear observation model
[11]. By implementing the E-step using AMP in the expectation maximization
(EM) based SBL method, matrix inversion can be avoided, leading to a
significant reduction in computational complexity. However, AMP does not work
well for a generic matrix such as non-zero mean, rank-deficient, correlated,
or ill-conditioned matrix $\mathbf{A}$ [12], resulting in divergence and poor
performance.

Many variants to AMP have been proposed to address the divergence issue and
achieve better robustness to a generic $\mathbf{A}$, such as the damped AMP
[12], swept AMP [13], generalized approximate message passing algorithm (GAMP)
with adaptive damping [14], vector AMP [15], orthogonal AMP [16], memory AMP
[17], convolutional AMP [18] and more. In [19], by incorporating damped
Gaussian generalized AMP (GGAMP) to the EM-based SBL method, a GGAMP-SBL
algorithm was proposed. Although the robustness of the approach is
significantly improved, it comes at the cost of slowing the convergence. In
addition, the algorithm still exhibits significant performance gap from the
support-oracle bound when the measurement matrix has relatively high
correlation, large condition number or non-zero mean.

For a general linear inverse problem, [20, 21] proposed to apply AMP to a
unitary transform of the original model, where the unitary matrix for the
transformation can be obtained by the singular value decomposition (SVD) of
$\mathbf{A}$. In the case of a circulant $\mathbf{A}$, the normalized discrete
Fourier transform matrix can be used for the unitary transformation, enabling
highly efficient implementation with the fast Fourier transform (FFT)
algorithm [22]. This leads to an AMP variant named unitary AMP (UAMP), which
was formerly called AMP with unitary transformation (UTAMP). 111 SVD plays an
important role in both UAMP and VAMP. In UAMP, SVD is used to obtain the
unitary transformed model. In VAMP, the linear minimum mean squared error
(LMMSE) estimator results in cubic complexity in each iteration, and VAMP
relies on SVD to implement the LMMSE estimator with low complexity. In this
work, we apply this concept to SBL, resulting in a new SBL algorithm called
UAMP-SBL. UAMP-SBL achieves more efficient sparse signal recovery with
significantly enhanced robustness, compared to the state-of-the-art AMP-based
SBL algorithm GGAMP-SBL [19].

To develop UAMP-SBL, we apply structured variational inference (SVI) [23],
[24], [25]. In particular, the formulated problem is represented by a factor
graph model, based on which approximate inference is implemented in terms of
structured variational message passing (SVMP) [24], [25], [26]. The use of
SVMP allows the incorporation of UAMP to the message passing algorithm to
handle the most computational intensive part of message computations with high
robustness and low complexity. In UAMP-SBL, a Gamma distribution is used as
the hyperprior for the precisions of the elements of $\mathbf{x}$. We propose
to tune the shape parameter of the Gamma distribution automatically during
iterations. We show by simulations that, in many cases  with a generic
measurement matrix, UAMP-SBL can still approach the support-oracle bound
closely. We also investigate SE-based performance prediction for UAMP-SBL and
analyze the impact of the shape parameter on SBL. In addition, the UAMP-SBL
algorithm is extended from single measurement vector (SMV) problems to
multiple measurement vector (MMV) problems [27], [28], [29]. Based on our
preliminary results in [30]222Compared to [30], we present a new derivation of
UAMP-SBL, extend it from SMV to MMV, and provide theoretical analyses and
comprehensive comparisons., UAMP-SBL was applied to inverse synthetic aperture
radar (ISAR) [31], where the measurement matrix can be highly correlated in
order to achieve high Doppler resolution. Real data experiments in [31]
demonstrate its superiority in terms of both recovery performance and speed.

The rest of the paper is organized as follows. We briefly introduce SBL and
(U)AMP in Section II. In Section III, UAMP-SBL is derived for SMV problems and
SE-based performance prediction for UAMP-SBL is also discussed. The impact of
the shape parameter is analyzed in Section IV. UAMP-SBL is extended to the MMV
setting in Section V. Numerical results are provided in Section VI, followed
by conclusions in Section VII.

Throughput the paper, we use boldface lowercase and uppercase letters to
represent column vectors and matrices, respectively. The superscript
$(\cdot)^{H}$ represents the conjugate transpose for a complex matrix, and the
transpose for a real matrix. We use $\mathbf{1}$ and $\mathbf{0}$ to denote an
all-one vector and an all-zero vector with proper sizes, respectively. The
notation $\mathcal{N}(\mathbf{x}|\bm{\mu},\bm{\Sigma})$ denotes a Gaussian
distribution of $\mathbf{x}$ with mean $\bm{\mu}$ and covariance
$\bm{\Sigma}$, and ${\mathrm{Ga}(\bm{\gamma}|\epsilon,\eta)}$ is a Gamma
distribution with shape parameter $\epsilon$ and rate parameter $\eta$. We use
$\left|\cdot\right|^{.2}$ to denote the element-wise magnitude squared
operation, and $\left\|\cdot\right\|$ for the $l_{2}$ norm. The notation
$\left\langle f(\mathbf{x})\right\rangle_{q(\mathbf{x})}$ denotes the
expectation of $f(\mathbf{x})$ with respect to probability density function
$q(\mathbf{x})$, and $E[\cdot]$ is the expectation over all random variables
involved in the brackets. We use $Diag(\mathbf{a})$ to represent a diagonal
matrix with elements of $\mathbf{a}$ on its diagonal, $Z_{m,n}$ is the
$(m,n)$th element of $\mathbf{Z}$, and $a_{n}$ is the $n$th element of vector
$\mathbf{a}$. The element-wise product and division of two vectors
$\mathbf{a}$ and $\mathbf{b}$ are written as $\mathbf{a}\cdot\mathbf{b}$ and
$\mathbf{a}./\mathbf{b}$, respectively. The superscript $\mathbf{a}^{t}$ is
the $t$th iteration in an iterative algorithm.

##  II Background

###  II-A Sparse Bayesian Learning

Consider recovering a length-$N$ sparse vector $\mathbf{x}$ from measurements

$\mathbf{y}=\mathbf{Ax}+\mathbf{w},\\\ $ (1)

where $\mathbf{y}$ is a measurement vector of length $M$, the measurement
matrix $\mathbf{A}$ has size $M\times N$, $\mathbf{w}$ denotes a Gaussian
noise vector with mean zero and covariance matrix $\beta^{-1}\mathbf{I}$, and
$\beta$ is the precision of the noise. It is assumed that the elements in
$\mathbf{x}$ are independent and the following two-layer sparsity-promoting
prior is used

$\displaystyle p(\mathbf{x}|\bm{\gamma})$ $\displaystyle=$
$\displaystyle\prod_{n}p(x_{n}|\gamma_{n})=\prod_{n}\mathcal{N}({x_{n}}|0,\gamma_{n}^{-1}),$
(2) $\displaystyle p(\bm{\gamma})$ $\displaystyle=$
$\displaystyle\prod_{n}p(\gamma_{n})=\prod_{n}\mathrm{Ga}(\gamma_{n}|\epsilon,\eta),$
(3)

i.e., the prior of $x_{n}$ is a Gaussian mixture

$p(x_{n})=\int\mathcal{N}({x_{n}}|0,\gamma_{n}^{-1})p(\gamma_{n})d\gamma_{n},$
(4)

where the precision vector
$\mathbf{\gamma}=[\gamma_{1},\gamma_{2},...,\gamma_{N}]^{H}$.

In the conventional SBL algorithm by Tipping [2], the precision vector
$\bm{\gamma}$ is learned by maximizing the a posteriori probability

$p(\bm{\gamma}|\mathbf{y})\propto p(\mathbf{y}|\bm{\gamma})p(\bm{\gamma}),$
(5)

where the marginal likelihood function

$\displaystyle p(\mathbf{y}|\bm{\gamma})=\int
p(\mathbf{y}|\mathbf{x})p(\mathbf{x}|\bm{\gamma})d\mathbf{x}.$ (6)

It can be shown that [2]

$\displaystyle\log{p(\mathbf{y}|\bm{\gamma})}\\!\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!-\frac{1}{2}\left(\log{|\mathbf{B}|}+\mathbf{y}^{H}\mathbf{B}^{-1}\mathbf{y}\right)+const$
(8) $\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}\log{|\mathbf{\Sigma}|}+\log{|Diag(\bm{\gamma})|}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\bm{\zeta}^{H}Diag(\bm{\gamma})\bm{\zeta}\big{)}+const1,$

where $const$ and $const1$ represent terms independent of $\bm{\gamma}$, and

$\displaystyle\mathbf{B}$ $\displaystyle=$
$\displaystyle\beta^{-1}\mathbf{I}+\mathbf{A}Diag(\bm{\gamma})^{-1}\mathbf{A}^{H},$
(9) $\displaystyle\mathbf{\Sigma}$ $\displaystyle=$
$\displaystyle\left(\beta\mathbf{A}^{H}\mathbf{A}+Diag({\bm{\gamma}})\right)^{-1},$
(10) $\displaystyle\bm{\zeta}$ $\displaystyle=$
$\displaystyle\beta\mathbf{\Sigma A}^{H}\mathbf{y}.$ (11)

The a posteriori probability of $\mathbf{x}$

$p(\mathbf{x}|\mathbf{y},\bm{\gamma})=\mathcal{N}(\mathbf{x}|\bm{\zeta},\bm{\Sigma}).$
(12)

By taking the logarithm of $p(\bm{\gamma}|\mathbf{y})$ and ignoring terms
independent of $\bm{\gamma}$, the learning of $\bm{\gamma}$ is to maximize the
following objective function [2]

$\displaystyle\mathcal{L}(\bm{\gamma})=\log{p(\mathbf{y}|\bm{\gamma})}+\sum_{n=1}^{N}(\epsilon\log{\gamma_{n}}-\eta\gamma_{n}).$
(13)

As the value of $\bm{\gamma}$ that maximizes $\mathcal{L}(\bm{\gamma})$ cannot
be obtained in a closed form, iterative re-estimation is employed by taking
advantage of (8), i.e., with a learned ${\bm{\gamma}}$ in the last iteration,
compute $\bm{\Sigma}$ and $\bm{\zeta}$ with (10) and (11), then update
$\bm{\gamma}$ by maximizing $\mathcal{L}(\bm{\gamma})$ with (8) used, which
leads to a closed form to update $\gamma_{n}$

$\gamma_{n}=(2{\epsilon}+1)/(2\eta+|\zeta_{n}|^{2}+\Sigma_{n,n}),n=1,...,N.$
(14)

In summary, Tipping’s SBL algorithm (which is called SBL hereafter) executes
the following iteration [2]:

$\displaystyle\\!\\!\\!\\!\\!\\!Repeat$
$\displaystyle~{}\mathbf{Z}=\left(\beta\mathbf{A}^{H}\mathbf{A}+Diag(\hat{\bm{\gamma}})\right)^{-1}$
(15) $\displaystyle~{}\mathbf{\hat{x}}=\beta\mathbf{ZA}^{H}\mathbf{y}$ (16)
$\displaystyle~{}\hat{\gamma}_{n}=(2{\epsilon}+1)/(2\eta+|\hat{x}_{n}|^{2}+Z_{n,n}),n=1,...,N.$
(17) $\displaystyle\\!\\!\\!\\!\\!\\!Until~{}terminated$

If the noise precision $\beta$ is unknown, its estimation can be incorporated
as well. The SBL algorithm can also be derived based on the EM algorithm [2],
[19]. The SBL algorithm requires a matrix inverse in (15) in each iteration.
This results in cubic complexity in each iteration, which can be prohibitive
for large scale problems. To address this issue, the implementation of the
E-step using AMP has been investigated. GGAMP was used in [19] to implement
the E-Step, where sufficient damping is used to enhance the robustness of the
algorithm against a generic measurement matrix. This leads to the GGAMP-SBL
algorithm with complexity significantly lower than that of SBL.

###  II-B (U)AMP

AMP was derived based on the loopy BP with Gaussian and Taylor-series
approximations [4], [10], which can be used to efficiently solve linear
inverse problems due to its low complexity. An issue with AMP is that it can
easily diverge in the case of a generic measurement matrix, such as
correlated, ill-conditioned, non-zero mean or rank-deficient $\mathbf{A}$
[12]. Inspired by the work in [22], it was shown in [20] that the robustness
of AMP is remarkably improved through simple pre-processing, i.e., performing
a unitary transformation to the original linear model [20], [21]. As any
matrix $\mathbf{A}$ has an SVD $\mathbf{A=U\Lambda V}$ with $\mathbf{U}$ and
$\mathbf{V}$ being two unitary matrices, performing a unitary transformation
with $\mathbf{U}^{H}$ leads to the following model

$\mathbf{r=\Phi x}+\bm{\omega},$ (18)

where $\mathbf{r=U}^{H}\mathbf{y}$,
$\mathbf{\Phi}=\mathbf{U}^{H}\mathbf{A}=\mathbf{\Lambda}\mathbf{V}$,
$\mathbf{\Lambda}$ is an $M\times N$ rectangular diagonal matrix, and
$\bm{\omega}=\mathbf{U}^{H}\mathbf{w}$ remains a zero-mean Gaussian noise
vector with the same covariance matrix $\beta^{-1}\mathbf{I}$. Applying the
vector step size AMP [10] with model (18) leads to the first version of UAMP
(called UAMPv1) shown in Algorithm 1.333By replacing $\mathbf{r}$ and
$\mathbf{\Phi}$ with $\mathbf{y}$ and $\mathbf{A}$ in Algorithm 1
respectively, the original AMP algorithm is recovered.

Algorithm 1 UAMP (UAMPv2 executes operations in [ ])

Initialize $\bm{\tau}_{x}^{(0)}(\mathrm{or}~{}{\tau}_{x}^{(0)})>0$ and
${{\mathbf{x}}^{(0)}}$. Set $\mathbf{s}^{(-1)}=\mathbf{0}$ and $t=0$. Define
vector $\bm{\lambda}=\mathbf{\Lambda\Lambda}^{H}\textbf{1}$.
Repeat

1: $\bm{\tau}_{p}$ = $\mathbf{|\Phi|}^{.2}\bm{\tau}^{t}_{x}$
$\left[\mathrm{or}~{}\bm{\tau}_{p}=\tau^{t}_{x}\bm{\lambda}\right]$

2:
$\mathbf{p}=\mathbf{\Phi}{{\mathbf{x}}^{t}}-\bm{\tau}_{p}\cdot\mathbf{s}^{t-1}$

3: $\bm{\tau}_{s}=\mathbf{1}./(\bm{\tau}_{p}+\beta^{-1}\mathbf{1})$

4: $\mathbf{s}^{t}=\bm{\tau}_{s}\cdot(\mathbf{r}-\mathbf{p})$

5: $\mathbf{1}./\bm{\tau}_{q}$ = $|\mathbf{\Phi}^{H}|^{.2}\bm{\tau}_{s}$
$\left[\mathrm{or}~{}\mathbf{1}./\bm{\tau}_{q}=(\frac{1}{N}\bm{\lambda}^{H}\bm{\tau}_{s})\mathbf{1}\right]$

6:
$\mathbf{q}={{\mathbf{x}}^{t}}+\bm{\tau}_{q}\cdot(\mathbf{\Phi}^{H}\mathbf{s}^{t})$

7: $\bm{\tau}_{x}^{t+1}$ = $\bm{\tau}_{q}\cdot
g_{x}^{\prime}(\mathbf{q},\bm{\tau}_{q})$
$\left[\mathrm{or}~{}\tau_{x}^{t+1}\\!=\\!\frac{1}{N}\mathbf{1}^{H}\left(\bm{\tau}_{q}\cdot
g_{x}^{\prime}(\mathbf{q},\tau_{q})\right)\right]$

8: ${{\mathbf{x}}^{t+1}}=g_{x}(\mathbf{q},\bm{\tau}_{q})$

9: $t=t+1$

Until terminated

We can apply an average operation to two vectors: $\bm{\tau}_{x}$ in Line 7
and $|\mathbf{\Phi}^{H}|^{.2}\bm{\tau}_{s}$ in Line 5 of UAMPv1 in Algorithm
1, leading to the second version of UAMP [20] (called UAMPv2), where the
operations in the brackets of Lines 1, 5 and 7 are executed (refer to [21] for
the derivation). Compared to AMP and UAMPv1, UAMPv2 does not require matrix-
vector products in Lines 1 and 5, so that the number of matrix-vector products
is reduced from 4 to 2 per iteration. This is a significant reduction in
computational complexity because the complexity of AMP-like algorithms is
dominated by matrix-vector products.

In the (U)AMP algorithms, $g_{x}(\mathbf{q},\bm{\tau}_{q})$ is related to the
prior of $\mathbf{x}$ and returns a column vector with the $n$th element
$[g_{x}(\mathbf{q},\bm{\tau}_{q})]_{n}$ given by

$[g_{x}(\mathbf{q},\bm{\tau}_{q})]_{n}=\frac{\int
x_{n}p(x_{n})\mathcal{N}(x_{n};q_{n},\tau_{q_{n}})dx_{n}}{\int
p(x_{n})\mathcal{N}(x_{n};q_{n},\tau_{q_{n}})dx_{n}},$ (19)

where we note that $p(x_{n})$ represents a general known prior for $x_{n}$.
The function $g_{x}^{\prime}(\mathbf{q},\bm{\tau}_{q})$ returns a column
vector and the $n$th element is denoted by
$[g_{x}^{\prime}(\mathbf{q},\bm{\tau}_{q})]_{n}$, where the derivative is
taken with respect to $q_{n}$.

##  III Sparse Bayesian Learning Using UAMP

###  III-A Problem Formulation and Approximate Inference

To enable the use of UAMP, we employ the unitary transformed model
$\mathbf{r=\Phi x}+\bm{\omega}$ in (18). As in many applications the noise
precision $\beta$ is unknown, its estimation is also considered. The joint
conditional distribution of ${\mathbf{x}}$, $\bm{\gamma}$ and $\beta$ can be
expressed as

$\displaystyle p(\mathbf{x},\bm{\gamma},\beta|\mathbf{r})\propto
p(\mathbf{r}|\mathbf{x},\beta)p(\mathbf{x}|\bm{\gamma})p(\bm{\gamma})p(\beta),$
(20)

where $p(\mathbf{x}|\bm{\gamma})$ and $p(\bm{\gamma})$ are given by (2) and
(3), respectively. We assume an improper prior $p(\beta)\propto 1/\beta$ for
the noise precision [2]. According to the transformed model (18),
$p(\mathbf{r}|\mathbf{x},\beta)=\mathcal{N}(\mathbf{r}|\mathbf{\Phi
x},\beta^{-1}\mathbf{I})$. Our aim is to find the marginal distribution
$p(\mathbf{x}|\mathbf{r})$. The a posteriori mean is then used as an estimate
of $\mathbf{x}$ in the sense of minimum mean squared error (MSE). However,
exact inference is intractable due to the high dimensional integration
involved, so we resort to approximate inference techniques.

Variational inference is a machine learning method for approximate inference,
and it has been widely used to approximate posterior densities for Bayesian
models [23], [24], [25]. In variational inference, a trial density function is
chosen and optimized by minimizing the Kullback-Leibler (KL) divergence
between the trial function and the true a posteriori function. Instead of
using fully factorized trial functions where all variables are assumed to be
independent (thereby likely resulting in poor approximations), more structured
factorizations can be used, leading to SVI algorithms. With graphical models,
SVI can be formulated as message-passing algorithms [24], [25], [26], which is
termed SVMP. In this work, SVMP is adopted because the use of SVMP facilitates
the incorporation of the message passing algorithm UAMP into SVMP. We will
show how UAMP can be used to handle the most computational intensive part of
message computations, enabling us to achieve low complexity and high
robustness. With SVMP, we can find an approximation to the marginal
distribution $p(\mathbf{x}|\mathbf{r})$, where an approximation to
$p(\bm{\gamma}|\mathbf{r})$ is also involved (the approximate inference for
$\mathbf{x}$ and $\bm{\gamma}$ is performed alternately).

We introduce an auxiliary variable $\mathbf{h}=\mathbf{\Phi x}$ to facilitate
the incorporation of UAMP, which is crucial to an efficient realization of
SBL. Then the conditional joint distribution is

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!p(\mathbf{x},\mathbf{h},\bm{\gamma},\beta|\mathbf{r})$
$\displaystyle\\!\\!\\!\\!\\!\propto
p(\mathbf{r}|\mathbf{h},\beta)p(\mathbf{h}|\mathbf{x})p(\mathbf{x}|\bm{\gamma})p(\bm{\gamma}|{\epsilon})p(\beta)$
$\displaystyle\\!\\!\\!\\!\\!=\prod_{m=1}^{M}\mathcal{N}(r_{m}|h_{m},{\beta^{-1}})\prod_{m=1}^{M}\delta(h_{m}-[\mathbf{\Phi}]_{m}{\mathbf{x}})$
$\displaystyle~{}\prod_{n=1}^{N}\mathcal{N}(x_{n}|0,\gamma_{n}^{-1})\prod_{n=1}^{N}\mathrm{Ga}(\gamma_{n}|\epsilon,\eta)p(\beta).$
(21)

To facilitate the derivation of the message passing algorithm, a factor graph
representation of the factorization in (21) is shown in Fig.1, where the local
functions $f_{\beta}(\beta)\propto 1/\beta$,
$f_{r_{m}}({r_{m}},h_{m},\beta)=\mathcal{N}(r_{m}|h_{m},\beta^{-1})$,
$f_{\delta_{m}}(h_{m},{\mathbf{x}})=\delta(h_{m}-[\mathbf{\Phi}]_{m}{\mathbf{x}})$,
$f_{x_{n}}({x_{n}},\gamma_{n})=\mathcal{N}(x_{n}|0,\gamma_{n}^{-1})$,
$f_{\gamma_{n}}(\gamma_{n})=\mathrm{Ga}(\gamma_{n}|\epsilon,\eta)$ and
$[\mathbf{\Phi}]_{m}$ is the $m$th row of matrix $\mathbf{\Phi}$.

Following SVI, we define the following structured trial function

$\tilde{q}(\mathbf{x},\mathbf{h},\bm{\gamma},\beta)=\tilde{q}(\beta)\tilde{q}(\mathbf{x},\mathbf{h})\tilde{q}(\bm{\gamma}).$
(22)

In terms of SVMP, the use of the above trial function corresponds to a
partition of the factor graph shown by the dotted boxes in Fig. 1, where
$\tilde{q}(\beta)$, $\tilde{q}(\mathbf{x},\mathbf{h})$ and
$\tilde{q}(\bm{\gamma})$ are associated with Subgraphs 1, 2 and 3,
respectively.

Figure 1: Factor graph of (21) for deriving UAMP-SBL.

As the KL divergence

$\displaystyle\mathcal{KL}\big{(}\tilde{q}(\beta)\tilde{q}(\mathbf{x},\mathbf{h})\tilde{q}(\bm{\gamma})||p(\mathbf{x},\mathbf{h},\bm{\gamma},\beta|\mathbf{r})\big{)},$
(23)

is minimized, it is expected that

$\displaystyle\tilde{q}(\mathbf{x},\mathbf{h})$ $\displaystyle\approx$
$\displaystyle p(\mathbf{x},\mathbf{h}|\mathbf{r}),$ (24)
$\displaystyle\tilde{q}(\bm{\gamma})$ $\displaystyle\approx$ $\displaystyle
p(\bm{\gamma}|\mathbf{r}),$ (25) $\displaystyle\tilde{q}(\beta)$
$\displaystyle\approx$ $\displaystyle p(\beta|\mathbf{r}).$ (26)

Integrating out $\mathbf{h}$ in (24), which corresponds to running BP in
Subgraph 2 (except the factor nodes connecting external variable nodes), we
have $\tilde{q}(\mathbf{x})\approx p(\mathbf{x}|\mathbf{r})$. Running BP in
Subgraph 2 involves the most intensive computations; fortunately it can be
handled efficiently and with high robustness using UAMP. The derivation of
UAMP-SBL is shown in Appendix A, and the algorithm is summarized in Algorithm
2.

Algorithm 2 UAMP-SBL

Unitary transform: $\mathbf{r=U}^{H}\mathbf{y}=\mathbf{\Phi x}+\bm{\omega}$,
where $\mathbf{\Phi=U}^{H}\mathbf{A}=\mathbf{\Lambda V}$, and $\mathbf{A}$ has
SVD $\mathbf{A=U\Lambda V}$.
Define vector $\bm{\lambda}=\mathbf{\Lambda\Lambda}^{H}\mathbf{1}$.
Initialization: ${\tau_{x}^{(0)}=1}$, $\hat{\mathbf{x}}^{(0)}=\textbf{0}$,
${{\epsilon}=0.001}$, $\hat{{\bm{\gamma}}}=\textbf{1}$, $\hat{\beta}=1$,
$\mathbf{s}=\mathbf{0}$, and $t=0$.

Do

1: $\bm{\tau}_{p}$ = ${\tau^{t}_{x}}\bm{\lambda}$

2: $\mathbf{p}=\mathbf{\Phi\hat{x}}^{t}-\bm{\tau}_{p}\cdot\mathbf{s}$

3: $\mathbf{v}_{h}=\bm{\tau}_{p}./(\bm{1}+\hat{\beta}\bm{\tau}_{p})$

4:
$\mathbf{\hat{h}}=(\hat{\beta}\bm{\tau}_{p}\cdot\mathbf{r}+\mathbf{p})./(\bm{1}+\hat{\beta}\bm{\tau}_{p})$

5:
$\hat{\beta}={M}/({||\mathbf{r}-\mathbf{\hat{h}}||^{2}+\mathbf{1}^{H}\mathbf{v}_{h}})$

6: $\bm{\tau}_{s}=\mathbf{1}./(\bm{\tau}_{p}+\hat{\beta}^{-1}\mathbf{1})$

7: $\mathbf{s}=\bm{\tau}_{s}\cdot(\mathbf{r}-\mathbf{p})$

8: $1/\tau_{q}=({1}/{N})\bm{\lambda}^{H}\bm{\tau}_{s}$

9: $\mathbf{q}=\mathbf{\hat{x}}^{t}+\tau_{q}\mathbf{\Phi}^{H}\mathbf{s}$

10:
${\tau^{t+1}_{x}}=(\tau_{q}/N)\mathbf{1}^{H}(\mathbf{1}./(\mathbf{1}+\tau_{q}\hat{{\bm{\gamma}}}))$

11:
$\hat{\mathbf{x}}^{t+1}=\mathbf{q}./(\mathbf{1}+\tau_{q}\hat{{\bm{\gamma}}})$

12:
$\hat{{\gamma}}_{n}=({2{\epsilon}+1})/(|{{\hat{x}}_{n}^{t+1}}|^{2}+{\tau^{t+1}_{x}}),n=1,...,N.$

13:
${\epsilon}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{\hat{\gamma}_{n}})-\frac{1}{N}\sum_{n}{\log{\hat{\gamma}}_{n}}}$

14: $t=t+1$

while
($||\hat{\mathbf{x}}^{t+1}-\hat{\mathbf{x}}^{t}||^{2}/||\hat{\mathbf{x}}^{t+1}||^{2}>\delta_{x}$
and $t<t_{max}$)

Regarding the UAMP-SBL in Algorithm 2, we have the following remarks:

  1. 1.

UAMPv2 is employed in Algorithm 2. Similarly, UAMPv1 can also be used. By
comparing UAMPv1 and UAMPv2, the differences lie in Lines 1, 8, 9 and 10 as
vectors $\bm{\tau}^{t}_{x}$ and $\bm{\tau}_{q}$ need to be used. The UAMP-SBL
algorithms with two version of UAMP deliver comparable performance, but UAMP-
SBL with UAMPv2 has lower complexity.

  2. 2.

In SBL with Gamma hyperprior, the shape parameter $\epsilon$ and the rate
parameter $\eta$ are normally chosen to be very small values [2], and
sometimes the value of the shape parameter $\epsilon$ is chosen empirically,
e.g., $\epsilon=1$ in [32]. In UAMP-SBL, we propose to tune the shape
parameter automatically (as shown in Line 13) with the following empirical
rule

${\epsilon}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{\hat{\gamma}_{n}})-\frac{1}{N}\sum_{n}{\log{\hat{\gamma}}_{n}}},$
(27)

i.e., $\epsilon$ is learned iteratively with the iteration, starting from a
small positive initial value. We note that, as the log function is concave,
the parameter $\epsilon$ in (27) is guaranteed to be non-negative. In Section
IV, we will show that the shape parameter $\epsilon$ in the SBL algorithms
functions as a selective amplifier for $\\{\gamma_{n}\\}$, and a proper
$\epsilon$ plays a significant role in promoting sparsity, leading to
considerable performance improvement.

Figure 2: SE and evolution trajectory of UAMP-SBL with a nonzero mean A
($N=10240$, $M=8192$, SNR = 50dB, sparsity rate $\rho=0.1$ and matrix mean
$\mu=10$).

###  III-B SE-Based Performance Prediction

In this section, leveraging UAMP SE, we study how to predict the performance
of UAMP-SBL empirically. We treat UAMP-SBL as UAMP with a special denoiser,
enabling the use of UAMP SE to predict the performance of UAMP-SBL. The
denoiser in the UAMP-SBL corresponds to Lines 10-13 of the UAMP-SBL algorithm
(Algorithm 2).

As (U)AMP decouples the estimation of vector $\mathbf{x}$, in the $t$th
iteration, we have the following pseudo observation model

$q_{n}^{t}=x_{n}+w^{t}_{n},$ (28)

where $q_{n}^{t}$ is the $n$th element of $\mathbf{q}$ in Line 9 of the UAMP-
SBL algorithm in the $t$th iteration, and $w^{t}_{n}$ denotes a Gaussian noise
with mean 0 and variance $\tau^{t}$, which is given as

$\displaystyle\tau^{t}=\frac{N}{\bm{1}^{H}\big{(}\bm{\lambda}./(v_{x}^{t}\bm{\lambda}+\beta^{-1}\bm{1})\big{)}}.$
(29)

Here $v^{t}_{x}$ is the average MSE of $\\{x_{n}\\}$ after denoising in the
$t$th iteration. As it is difficult to obtain a closed form for the average
MSE, we simulate the denoiser with the additive Gaussian noise model (28) by
varying the variance of noise $\tau^{t}$ (or the SNR), so that we can get a
“function” in terms of a table, with the variance of the noise as the input
and the MSE as the output, i.e.,

$v_{x}=\phi(\tau).$ (30)

The function $\phi(\cdot)$ is independent of the measurement matrix
$\mathbf{A}$. The performance of UAMP-SBL can be predicted using the following
iteration with the initialization of $v_{x}$:

$\displaystyle Repeat$ (31)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\tau=\frac{N}{\bm{1}^{H}\big{(}\bm{\lambda}./(v_{x}\bm{\lambda}+\beta^{-1}\bm{1})\big{)}}\triangleq\psi(v_{x})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}v_{x}=\phi(\tau)$ $\displaystyle
Until~{}terminated$

We show the predicted performance, simulated performance in terms of
normalized MSE (NMSE, which is defined in (43)) and the evolution trajectory
of UAMP-SBL in Fig. 2 for a non-zero mean  measurement matrix $\mathbf{A}$. It
can be seen that the predicted performance matches well the simulated
performance.

###  III-C Computational Complexity

UAMP-SBL works well with a simple single loop iteration, which is in contrast
to the double loop iterative algorithm GGAMP-SBL [19]. The complexity of UAMP-
SBL (with UAMPv2) is dominated by two matrix-vector product operations in Line
2 and Line 9, i.e., $\mathcal{O}(MN)$ per iteration. The algorithm typically
converges fast and delivers outstanding performance as shown in Section VI.
UAMP-SBL involves an SVD , but it only needs to be computed once and may be
carried out off-line. The complexity of economic SVD is
$\mathcal{O}(\min\\{M^{2}N,MN^{2}\\})$. Note that for the runtime comparison
in Section VI, we do not assume off-line SVD computation, and the time
consumed by SVD is counted for UAMP-SBL.

##  IV Impact of the Shape Parameter $\epsilon$ in SBL

In this section, we analyze the impact of the hyperparameter $\epsilon$ on the
convergence of SBL. We focus on the case of an identity matrix $\mathbf{A}$.
The same results for a general $\mathbf{A}$ are demonstrated numerically.

We consider the conventional SBL algorithm ($\eta$ is set to be zero) [2]. In
the case of identity matrix $\mathbf{A}$, it reduces to

$\displaystyle\\!\\!\\!\\!\\!\\!Repeat$
$\displaystyle~{}~{}~{}~{}~{}Z_{n,n}=\left(\beta+{\gamma}_{n}^{t}\right)^{-1}$
$\displaystyle~{}~{}~{}~{}~{}\hat{x}_{n}=\beta Z_{n,n}y_{n}$ (32)
$\displaystyle~{}~{}~{}~{}~{}\gamma_{n}^{t+1}=(2{\epsilon}+1)/(|\hat{x}_{n}|^{2}+Z_{n,n})$
$\displaystyle\\!\\!\\!\\!\\!\\!Until~{}terminated$

Here note that in the above iteration we initialize $\gamma_{n}^{(0)}>0$. The
iteration in terms of $\gamma_{n}$ has a closed form

$\displaystyle\gamma^{t+1}_{n}$
$\displaystyle=\frac{2\epsilon+1}{(\beta(\beta+\gamma^{t}_{n})^{-1}y_{n})^{2}+(\beta+\gamma^{t}_{n})^{-1}}$
$\displaystyle=(2\epsilon+1)\frac{(\beta+\gamma^{t}_{n})^{2}}{(\beta
y_{n})^{2}+\beta+\gamma^{t}_{n}}$ (33) $\displaystyle\triangleq
g_{\epsilon}(\gamma^{t}_{n}).$

Next, we investigate the impact of $\epsilon$ on the convergence behavior and
fixed point of the iteration (33) when $\epsilon=0$ or $\epsilon$ takes a
positive value.

For the iteration (33) with a small positive initial value $\gamma_{n}^{(0)}$,
we have the following proposition and theorem.

Proposition 1: When $\epsilon=0$, if $\beta y^{2}_{n}>1$, $\gamma_{n}^{t}$
converges to a stable fixed point

$\gamma^{\prime}_{n}=\frac{\beta}{\beta y^{2}_{n}-1};$ (34)

if $\beta y^{2}_{n}\leq 1$, $\gamma_{n}^{t}$ goes to $+\infty$.

###### Proof.

See Appendix B. ∎

Theorem 1:  When $\epsilon>0$, if $\beta
y^{2}_{n}>1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, $\gamma_{n}^{t}$
converges to a stable fixed point

$\gamma_{n(a)}=\frac{2\beta(1+2\epsilon)}{\beta
y_{n}^{2}-4\epsilon-1+\sqrt{\beta^{2}y_{n}^{4}-8\epsilon\beta y_{n}^{2}-2\beta
y_{n}^{2}+1}};$ (35)

if $\beta y^{2}_{n}<1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$,
$\gamma_{n}^{t}$ goes to $+\infty$.

###### Proof.

See Appendix C. ∎

Based on Proposition 1 and Theorem 1, we make the following remarks:

  1. 1.

If $\beta y^{2}_{n}\leq 1$, for both $\epsilon=0$ and $\epsilon>0$,
$\gamma^{t}_{n}$ goes to $+\infty$. However, a positive $\epsilon$ accelerates
the move of $\gamma^{t}_{n}$ towards $+\infty$. This can be shown as follows.
As $\beta>0$ and $\beta y^{2}_{n}\leq 1$, we have $(\beta
y_{n})^{2}\leq\beta$. Hence, from (33)

$\displaystyle\gamma^{t+1}_{n}=g_{\epsilon}(\gamma^{t}_{n})$
$\displaystyle\geq(2\epsilon+1)\frac{(\beta+\gamma^{t}_{n})^{2}}{2\beta+\gamma^{t}_{n}}$
$\displaystyle=(2\epsilon+1)\left(\gamma^{t}_{n}+\frac{\beta^{2}}{2\beta+\gamma^{t}_{n}}\right)$
$\displaystyle>(2\epsilon+1)\gamma^{t}_{n}.$ (36)

From (36), compared to $\epsilon=0$, a positive value of $\epsilon$ moves
$\gamma^{t}_{n}$ towards infinity more quickly. Considering a fixed number of
iterations, a positive value of $\epsilon$ can be significant because the
precision can reach a large value much faster.

  2. 2.

When $\epsilon=0$, $\gamma^{t}_{n}$ converges to a finite fixed point if
$\beta y^{2}_{n}>1$. In contrast, when $\epsilon>0$, $\gamma^{t}_{n}$ goes to
$+\infty$ if $\beta
y^{2}_{n}\in(1,1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2})$. This is an
additional range for $\gamma^{t}_{n}$ to go to infinity. Hence, a positive
$\epsilon$ is stronger in terms of promoting sparsity, compared to
$\epsilon=0$.

  3. 3.

When $\epsilon>0$, if $\beta
y^{2}_{n}=1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, $\gamma^{t}_{n}$ may
converge or diverge because the iteration has a unique neutral fixed point as
shown in Theorem 1.

  4. 4.

When $\beta y^{2}_{n}>1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$,
$\gamma^{t}_{n}$ always converges to a fixed point. Based on (34) and (35),
the ratio of the precisions obtained with $\epsilon>0$ and $\epsilon=0$ is
given by

$\frac{\gamma_{n(a)}}{\gamma^{\prime}_{n}}=\frac{2(1+2\epsilon)}{1-\frac{4\epsilon}{\beta
y^{2}_{n}-1}+\sqrt{\left(1-\frac{4\epsilon}{\beta
y^{2}_{n}-1}\right)^{2}-\frac{8\epsilon(1+2\epsilon)}{({\beta
y^{2}_{n}-1})^{2}}}}.$ (37)

The ratio is a function of $\beta y^{2}_{n}$, and

${\gamma_{n(a)}}/{\gamma^{\prime}_{n}}\approx 1+2\epsilon,$ (38)

if $\beta y^{2}_{n}$ is relatively large.

Figure 3: Ratio of precisions with different $\epsilon$. Figure 4: Precisions
and their ratios ($\mathbf{A}$ is an identity matrix).  Figure 5: Precisions
and their ratios ($\mathbf{A}$ is i.i.d Gaussian).

The ratios of the precisions versus $\beta y^{2}_{n}$ are shown in Fig. 3,
where they are not shown for $\beta
y^{2}_{n}<1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$ as they are infinity
when $1<\beta y^{2}_{n}<1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, and
undefined when $\beta y^{2}_{n}\leq 1$ (see the above remarks). It can be seen
that the precision obtained with $\epsilon=0$ is amplified depending on the
value of $\beta y^{2}_{n}$. The smaller the value of $\beta y^{2}_{n}$, the
larger the amplification for the corresponding precision (in the case of
$\beta y^{2}_{n}\leq 1$, the ratios are undefined. However, considering a
fixed number of iterations, the ratios can be large as $\gamma^{t}_{n}$ with a
positive $\epsilon$ goes to infinity much quicker). Note that
$y_{n}=x_{n}+w_{n}$ and $\beta$ is the noise precision. Hence, if $\beta
y^{2}_{n}$ is a small value, it is highly likely that the corresponding
$x_{n}$ is zero, hence the precision $\gamma_{n}$ should go to infinity. If
$\beta y^{2}_{n}$ is a large value, it is highly likely that the corresponding
$x_{n}$ is non-zero, hence $\gamma_{n}$ should be a finite value. We see that
a positive $\epsilon$ tends to a sparser solution, and a proper value of
$\epsilon$ leads to much better recovery performance, compared to
$\epsilon=0$.

The precisions of the elements of the sparse vector obtained by the SBL
algorithm with $\epsilon=1.5$ and $\epsilon=0$ are shown in Fig. 4, where
$\mathbf{A}$ is an identity matrix with size $10000\times 10000$, the sparsity
rate of the signal is 0.1, and SNR = 50dB. It can be seen that the precisions
with $\epsilon=1.5$ are separated into two groups more clearly, and the ratios
for the small precisions are roughly 4 (i.e., $1+2\epsilon$), while other
precisions are amplified significantly. Although the above analysis is for an
identity matrix $\mathbf{A}$, it is interesting that the same results are
observed for a general matrix $\mathbf{A}$ as demonstrated numerically in Fig.
5, where $\mathbf{A}$ is an i.i.d Gaussian matrix with size $5000\times
10000$, the non-zero shape parameter $\epsilon=1.5$, and the sparsity rate and
the SNR are the same as the case of identity matrix. (Similar observations are
observed for other matrices). We see that the small precisions are also
roughly amplified by 4 times while others are amplified significantly, leading
to two well-separated groups.

Figure 6: Performance of the conventional SBL. (a) Gaussian matrix; (b)
correlated matrix with $c=0.3$; (c) low-rank matrix with $R/N=0.6$.

It is noted that the value of $\epsilon$ should be determined properly. If the
matrix $\mathbf{A}$ and the sparsity rate of $\mathbf{x}$ are given, we can
find a proper value for $\epsilon$ through trial and error. However, this is
inconvenient, and the sparsity rate of the signal may not be available. We
found the empirical equation (27) to determine the value of $\epsilon$. Next,
we examine its effectiveness with the SBL algorithm. Plugging the shape
parameter update rule (27) to the conventional SBL algorithm leads to the
following iterative algorithm (assuming the noise precision $\beta$ is known):

$\displaystyle Repeat$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{Z}=\left(\beta\mathbf{A}^{H}\mathbf{A}+Diag(\hat{\bm{\gamma}})\right)^{-1}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{\hat{x}}=\beta\mathbf{ZA}^{H}\mathbf{y}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\hat{\gamma}_{n}=(2{\epsilon}+1)/(|\hat{x}_{n}|^{2}+Z_{n,n}),n=1,...,N$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}{\epsilon}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{\hat{\gamma}_{n}})-\frac{1}{N}\sum_{n}{\log{\hat{\gamma}}_{n}}}$
$\displaystyle Until~{}terminated$

To demonstrate the effectiveness of the shape parameter update rule (27), we
compare the performance of the conventional SBL algorithm with and without
shape parameter update. The results are shown in Fig. 6, where the SNR is
50dB, the size of the measurement matrix is $800\times 1000$, and the sparsity
rate $\rho=0.1$. In this figure, the support-oracle bound is also shown for
reference. The matrices in (a), (b), and (c) are respectively i.i.d. Gaussian,
correlated and low-rank matrices (refer to Section VI for their generations).
It can be seen that there is a clear gap between the performance of the
conventional SBL and the bounds, and with shape parameter updated with our
rule, the SBL algorithm attains the bound. It is worth mentioning the
empirical finding in [19], i.e., replacing the noise variance $\beta^{-1}$
with $3\beta^{-1}$ can lead to better performance of GGAMP-SBL [19]. We use
this for the conventional SBL algorithm and the performance is also included
in Fig. 6. We see that it also leads to substantial performance improvement,
but its performance is inferior to that of SBL with updated $\epsilon$ using
(27). Moreover, in many cases, the noise variance is unknown, and it may be
hard to determine its value accurately. In contrast, our empirical update of
$\epsilon$ does not require any additional information.

##  V Extension to MMV

In this section, we extend UAMP-SBL to the MMV setting, where the relation
among the sparse vectors is exploited, e.g., common support and temporal
correlation.

###  V-A UAMP-SBL for MMV

The objective on an MMV problem is to recover a collection of length-$N$
sparse vectors
$\mathbf{X}=\left[\mathbf{x}^{(1)},\mathbf{x}^{(2)},...,\mathbf{x}^{(L)}\right]$
from $L$ noisy length-$M$ measurement vectors
$\mathbf{Y}=\left[\mathbf{y}^{(1)},\mathbf{y}^{(2)},...,\mathbf{y}^{(L)}\right]$
with the following model

$\displaystyle\mathbf{Y}=\mathbf{A}\mathbf{X}+\mathbf{W},$ (39)

where we assume that the $L$ vectors $\\{\mathbf{x}^{(l)}\\}$ share a common
support (i.e., joint sparsity), $\mathbf{A}$ is a known measurement matrix
with size $M\times N$, and $\mathbf{W}$ denotes an i.i.d. Gaussian noise
matrix with the elements having mean zero and precision $\beta$.

Figure 7: Factor graph representation of (41).

With the SVD $\mathbf{A=U\Lambda V}$, a unitary transformation with
$\mathbf{U}^{H}$ to (39) can be performed, i.e.,

$\mathbf{R=\Phi X}+\bm{\Omega},$ (40)

where
$\mathbf{R=U}^{H}\mathbf{Y}=\left[\mathbf{r}^{(1)},\mathbf{r}^{(2)},...,\mathbf{r}^{(L)}\right]$,
$\mathbf{\Phi=U}^{H}\mathbf{A}=\mathbf{\Lambda V}$ and
$\bm{\Omega}=\mathbf{U}^{H}\mathbf{W}$ is still white and Gaussian with mean
zero and precision $\beta$. Define $\mathbf{h}^{(l)}=\mathbf{\Phi x}^{(l)}$
and $\mathbf{H}=[\mathbf{h}^{(1)},...,\mathbf{h}^{(L)}]$. Then we have the
following joint distribution

$\displaystyle p({\mathbf{X,H}},\bm{\gamma},\beta|{\mathbf{R}})$
$\displaystyle\propto\prod_{l=1}^{L}p(\mathbf{r}^{(l)}|\mathbf{h}^{(l)},\beta)p(\mathbf{h}^{(l)}|\mathbf{x}^{(l)})p(\mathbf{x}^{(l)}|\bm{\gamma})p(\bm{\gamma})p(\beta)$
$\displaystyle=\prod_{l=1}^{L}\prod_{m=1}^{M}\mathcal{N}({r}_{m}^{(l)}|{h}_{m}^{(l)},{\beta^{-1}})\delta({h}_{m}^{(l)}-[\mathbf{\Phi}]_{m}{\mathbf{x}^{(l)}})$
$\displaystyle\times\prod_{l=1}^{L}\prod_{n=1}^{N}\mathcal{N}({x}_{n}^{(l)}|0,\gamma_{n}^{-1})\prod_{n=1}^{N}\mathrm{Ga}(\gamma_{n}|\epsilon,\eta)p(\beta).$
(41)

Define factors
$f_{\mathbf{r}}^{(l)}({\mathbf{r}^{(l)}},\mathbf{h}^{(l)},\beta)=\prod_{m}\mathcal{N}({r}_{m}^{(l)}|{h}_{m}^{(l)},\beta)$,
$f_{\bm{\delta}}^{(l)}(\mathbf{h}^{(l)},\mathbf{x}^{(l)})=\prod_{m}\delta({h}_{m}^{(l)}-[\mathbf{\Phi}]_{m}\mathbf{x}^{(l)})$,
$f_{\beta}(\beta)\propto 1/\beta$,
$f_{\mathbf{x}}^{(l)}(\mathbf{x}^{(l)},\bm{\gamma})=\prod_{n}\mathcal{N}(x_{n}^{(l)}|0,\gamma_{n}^{-1})$,
and
$f_{\bm{\gamma}}(\bm{\gamma},\epsilon)=\prod_{n}\mathrm{Ga}({\gamma_{n}}|\epsilon,\eta)$
denotes the hyperprior of the hyperparameters $\\{\gamma_{n}\\}$. The factor
graph representation of (41) is shown in Fig. 7444The vector variable node
$\bm{\gamma}$ is used in the factor graph to make it neat. We note that each
entry $x^{(l)}_{n}$ of $\bm{x}^{(l)}$ is connected to $\gamma_{n}$ through the
function node between them., based on which the message passing algorithm can
be derived. The message updates related to $\mathbf{x}^{(l)}$ and
$\mathbf{h}^{(l)}$ are the same as those for the SMV case and can be computed
in parallel. The difference lies in the computations of $\hat{\beta}$ and
$\hat{\bm{\gamma}}$, and the relevant derivations are shown in Appendix D. The
UAMP-SBL for MMV is summarized in Algorithm 3, where UAMPv2 is employed. The
complexity of the algorithm is $\mathcal{O}(MNL)$ per iteration.

Algorithm 3 UAMP-SBL for MMV

Unitary transform:
$\mathbf{R=U}^{H}\mathbf{Y}=\mathbf{\Phi}{\mathbf{X}}+{\mathbf{W}}$, where
$\mathbf{\Phi=U}^{H}\mathbf{A}=\mathbf{\Lambda V}$, and $\mathbf{A}$ has SVD
$\mathbf{A=U\Lambda V}$.
Define vector $\bm{\lambda}=\mathbf{\Lambda\Lambda}^{H}\mathbf{1}$.
Initialization: $\forall l$; ${\tau}_{x}^{l(0)}=1$,
$\hat{\mathbf{x}}^{l(0)}=\textbf{0}$, ${\epsilon^{\prime}}=0.001$,
$\hat{{\bm{\gamma}}}=\textbf{1}$, $\hat{\beta}=1$,
$\mathbf{s}^{l}=\mathbf{0}$, and $t=0$.
Do

1: $\forall l$; ${\bm{\tau}_{p}}^{l}$ = $\tau^{l(t)}_{x}\bm{\lambda}$

2: $\forall l$;
$\mathbf{p}^{l}=\mathbf{\Phi}\hat{\mathbf{x}}^{l(t)}-\bm{\tau}_{p}^{l}\cdot\mathbf{s}^{l}$

3: $\forall l$;
${\mathbf{v}^{l}_{h}}=\bm{\tau}_{p}^{l}./(\bm{1}+\hat{\beta}\bm{\tau}_{p}^{l})$

4: $\forall l$;
$\mathbf{\hat{h}}^{l}=(\hat{\beta}\bm{\tau}_{p}^{l}\cdot\mathbf{r}^{l}+\mathbf{p}^{l})./(\bm{1}+\hat{\beta}\bm{\tau}_{p}^{l})$

5:
$\hat{\beta}={LM}/(\sum_{l}({||\mathbf{r}^{l}-\mathbf{\hat{h}}^{l}||^{2}+\mathbf{1}^{H}\mathbf{v}_{h}^{l}}))$;

6: $\forall l$;
$\bm{\tau}_{s}^{l}=\mathbf{1}./(\bm{\tau}_{p}^{l}+{\hat{\beta}}^{-1}\mathbf{1})$

7: $\forall l$;
$\mathbf{s}^{l}=\bm{\tau}^{l}_{s}\cdot(\mathbf{r}^{l}-\mathbf{p}^{l})$

8: $\forall l$; $1/\tau_{q}^{l}=({1}/{N})\bm{\lambda}^{H}\bm{\tau_{s}}^{l}$

9: $\forall l$;
$\mathbf{q}^{l}=\hat{\mathbf{x}}^{l(t)}+\tau_{q}^{l}(\mathbf{\Phi}^{H}\mathbf{s}^{l})$

10: $\forall l$;
$\tau_{x}^{l(t+1)}=(\tau_{q}^{l}/N)\mathbf{1}^{H}(\mathbf{1}./(\mathbf{1}+\tau_{q}^{l}\hat{{\bm{\gamma}}}))$

11: $\forall l$;
$\hat{\mathbf{x}}^{l(t+1)}=\mathbf{q}^{l}./(\mathbf{1}+\tau_{q}^{l}\hat{{\bm{\gamma}}})$

12:
$\hat{{\gamma}}_{n}=\frac{2{\epsilon^{\prime}}+1}{(1/L)\sum^{L}_{l=1}(|{{\hat{x}}_{n}^{l(t+1)}}|^{2}+\tau_{x}^{l(t+1)})},n=1,...,N.$

13:
${\epsilon^{\prime}}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{{\hat{\gamma}_{n}}})-\frac{1}{N}\sum_{n}{\log{{\hat{\gamma}_{n}}}}}$

14: $t=t+1$

while
$\frac{1}{L}\sum_{l=1}^{L}(||\hat{\mathbf{x}}^{l(t+1)}-\hat{\mathbf{x}}^{l(t)}||^{2}/||\hat{\mathbf{x}}^{l(t+1)}||^{2})>\delta_{x}$
and $t<t_{max}$)

###  V-B UAMP-TSBL

With the assumption of a common sparsity profile shared by all sparse vectors,
we further consider exploiting the temporal correlation that exists between
the non-zero elements. The messages update related to $\mathbf{h}^{(l)}$,
${\epsilon}$ and $\beta$ are the same as those for the MMV case, where no
temporal correlation between non-zero elements is assumed. As the correlation
is considered, the differences from the UAMP-SBL MMV algorithm lie in the
computations of $\hat{\gamma}_{n}$ and $\mathbf{x}^{(l)}$.

As in [19], we use an AR(1) process [33] to model the correlation between
$x^{(l)}_{n}$ and $x^{(l-1)}_{n}$, i.e.,

$\displaystyle x^{(l)}_{n}$ $\displaystyle=\alpha
x^{(l-1)}_{n}+\sqrt{1-\alpha^{2}}\vartheta^{(l)}_{n}$ (42) $\displaystyle
p(x^{(l)}_{n}|x^{(l-1)}_{n})$ $\displaystyle=\mathcal{N}(x^{(l)}_{n}|\alpha
x^{(l-1)}_{n},(1-\alpha^{2}){{\gamma}^{-1}_{n}}),l>1$ $\displaystyle
p(x^{(1)}_{n})$
$\displaystyle=\mathcal{N}(x^{(1)}_{n}|0,{{\gamma}^{-1}_{n}}),$

where $\alpha\in(-1,1)$ is the temporal correlation coefficient and
$\vartheta^{(l)}_{n}\sim\mathcal{N}(0,{{\gamma}^{-1}_{n}})$. Due to the
temporal correlation, the conditional prior distribution for the vector
$\mathbf{x}^{(l)}$ changes. We redefine the factors
$\\{f_{x^{(l)}_{n}}(x^{(l)}_{n},{\gamma_{n}})\\}$, i.e.,
$f_{x^{(l)}_{n}}(x^{(l)}_{n},{\gamma_{n}})=p(x^{(l)}_{n}|x^{(l-1)}_{n})$ for
$l>1$ and $f_{x^{(1)}_{n}}(x^{(1)}_{n},{\gamma_{n}})=p(x^{(1)}_{n})$. Thus,
each $x_{n}^{(l)}$ is connected to the factor nodes
$f^{(l)}_{x_{n}}(x^{(l)}_{n}|{\gamma_{n}})$,
$f^{(l+1)}_{x_{n}}(x^{(l+1)}_{n}|{\gamma_{n}})$ and
$\\{f^{(l)}_{\delta_{m}}(h^{(l)}_{m}|\mathbf{x}^{(l)}),\forall m\\}$. The
factor graph characterizing the temporal correlation is shown in Fig. 8. The
remaining part of the graph is omitted as it is the same as that of the MMV
case without temporal correlation. The derivation of the extra message passing
for the UAMP-TSBL algorithm is shown in Appendix E, and the algorithm is
summarized in Algorithm 4. UAMP-TSBL is an extension of the UAMP-SBL algorithm
for MMV (Algorithm 3). The complexity of the UAMP-TSBL algorithm is also
dominated by matrix-vector multiplications, and it is $\mathcal{O}(MNL)$ per
iteration.

Figure 8: The additional factor graph for deriving UAMP-TSBL. Algorithm 4
UAMP-TSBL

Unitary transform:
$\mathbf{R=U}^{H}\mathbf{Y}=\mathbf{\Phi}{\mathbf{X}}+\mathbf{W}$, where
$\mathbf{\Phi=U}^{H}\mathbf{A}=\mathbf{\Lambda V}$, and $\mathbf{A}$ has SVD
$\mathbf{A=U\Lambda V}$.
Define vector $\bm{\lambda}=\mathbf{\Lambda\Lambda}^{H}\mathbf{1}$.
Initialization: $\forall l$: ${\tau}_{x}^{l(0)}={1}$,
$\hat{\mathbf{x}}^{l(0)}=\mathbf{0}$,
$\mathbf{q}^{l}={\mathbf{0}}$,${\bm{\tau}_{q}^{{l(0)}}}=\mathbf{1}$,
$\bm{\xi}^{l(0)}=\mathbf{0}$, $\bm{\psi}^{l(0)}={\mathbf{1}}$,
$\bm{\theta}^{l(0)}=\mathbf{0}$, $\bm{\phi}^{l(0)}={\mathbf{1}}$,
$\mathbf{s}^{l(-1)}=\mathbf{0}$, ${\epsilon^{\prime}}=0.001$,
$\hat{{\bm{\gamma}}}^{(0)}=\textbf{1}$, $\hat{\beta}=1$, and $t=0$.
Do

1: $\bm{\xi}^{1}=\mathbf{0}$

2: $\bm{\psi}^{1}=\mathbf{1}./\hat{\bm{\gamma}}^{(t)}$

3: for $l=2,...,L$

4:
$\quad\bm{\xi}^{l}=\alpha\left(\frac{\mathbf{q}^{l-1}}{\bm{\tau}_{q}^{l-1}}+\frac{\bm{\xi}^{l-1}}{\bm{\psi}^{l-1}}\right)\cdot\left(\frac{\bm{\tau}_{q}^{l-1}\cdot\bm{\psi}^{l-1}}{\bm{\tau}_{q}^{l-1}+\bm{\psi}^{l-1}}\right)$

5:
$\quad\bm{\psi}^{l}=\alpha^{2}\left(\frac{\bm{\tau}_{q}^{l-1}\cdot\bm{\psi}^{l-1}}{\bm{\tau}_{q}^{l-1}+\bm{\psi}^{l-1}}\right)+(1-\alpha^{2})/\hat{\bm{\gamma}}^{(t)}$

6: end

7: for $l=1,...,L$

8: $\quad{\bm{\tau}_{p}}^{l}$ = $\tau^{l(t)}_{x}\bm{\lambda}$

9:
$\quad\mathbf{p}^{l}=\mathbf{\Phi}\hat{\mathbf{x}}^{l(t)}-\bm{\tau}_{p}^{l}\cdot\mathbf{s}^{l(t-1)}$

10:
$\quad{\mathbf{v}^{l}_{h}}=\bm{\tau}_{p}^{l}./(\bm{1}+\hat{\beta}\bm{\tau}_{p}^{l})$

11:
$\quad\mathbf{\hat{h}}^{l}=(\hat{\beta}\bm{\tau}_{p}^{l}\cdot\mathbf{r}^{l}+\mathbf{p}^{l})./(\bm{1}+\hat{\beta}\bm{\tau}_{p}^{l})$

12: end

13:
$\hat{\beta}={LM}/(\sum_{l}({||\mathbf{r}^{l}-\mathbf{\hat{h}}^{l}||^{2}+\mathbf{1}^{H}\mathbf{v}_{h}^{l}}))$

14: for $l=1,...,L$

15:
$\quad\bm{\tau}_{s}^{l}=\mathbf{1}./(\bm{\tau}_{p}^{l}+\hat{\beta}^{-1}\mathbf{1})$

16:
$\quad\mathbf{s}^{l(t)}=\bm{\tau}^{l}_{s}\cdot(\mathbf{r}^{l}-\mathbf{p}^{l})$

17: $\quad 1/\tau_{q}^{l}=({1}/{N})\bm{\lambda}^{H}\bm{\tau_{s}}^{l}$

18:
$\quad\mathbf{q}^{l}=\hat{\mathbf{x}}^{l(t)}+\tau_{q}^{l}(\mathbf{\Phi}^{H}\mathbf{s}^{l(t)})$

19:
$\quad{\tau}_{x}^{l(t+1)}=(1/N)\mathbf{1}^{H}(\mathbf{1}./(\mathbf{1}./\bm{\tau}^{l}_{q}+\mathbf{1}./\bm{\phi}^{l}+\mathbf{1}./\bm{\psi}^{l}))$

20:
$\quad\hat{\mathbf{x}}^{l(t+1)}={\tau}_{x}^{l(t+1)}(\mathbf{q}^{l}./\bm{\tau}^{l}_{q}+\bm{\theta}^{l}./\bm{\phi}^{l}+\bm{\xi}^{l}./\bm{\psi}^{l})$

21: end

22: $\bm{\theta}^{L-1}=\frac{1}{\alpha}\mathbf{q}^{L}$

23:
$\bm{\phi}^{L-1}=\frac{1}{\alpha^{2}}\left(\bm{\tau}_{q}^{L}+(1-\alpha^{2})/\hat{\bm{\gamma}}^{(t)}\right)$

24: for $l=L-2,...,1$

25:
$\quad\bm{\theta}^{l}=\frac{1}{\alpha}\left(\frac{\mathbf{q}^{l+1}}{\bm{\tau}_{q}^{l+1}}+\frac{\bm{\theta}^{l+1}}{\bm{\phi}^{l+1}}\right)\cdot\left(\frac{\bm{\tau}_{q}^{l+1}\bm{\phi}^{l+1}}{\bm{\tau}_{q}^{l+1}+\bm{\phi}^{l+1}}\right)$

26:
$\quad\bm{\phi}^{l}=\frac{1}{\alpha^{2}}\left(\frac{\bm{\tau}_{q}^{l+1}\bm{\phi}^{l+1}}{\bm{\tau}_{q}^{l+1}+\bm{\phi}^{l+1}}+(1-\alpha^{2})/\hat{\bm{\gamma}}^{(t)}\right)$

27: end

28:
$\hat{\bm{\gamma}}^{(t+1)}=L({2{\epsilon^{\prime}}+1})/[|{\mathbf{\hat{x}}^{1(t+1)}}|^{2}+{\tau}_{x}^{1(t+1)}\mathbf{1}$$+\frac{1}{1-\alpha^{2}}\sum_{l=2}^{L}(|{\mathbf{\hat{x}}^{l(t+1)}}|^{2}+{\tau}_{x}^{l(t+1)}\mathbf{1})$$+\frac{\alpha^{2}}{1-\alpha^{2}}\sum_{l=1}^{L-1}(|{\mathbf{\hat{x}}^{l(t+1)}}|^{2}+{\tau}_{x}^{l(t+1)}\mathbf{1})-\frac{2\alpha}{1-\alpha^{2}}\sum_{l=2}^{L}({\mathbf{\hat{x}}^{l(t+1)}}\cdot{\mathbf{\hat{x}}^{(l-1)(t+1)}})]$

29:
${\epsilon^{\prime}}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{{\hat{\gamma}^{(t+1)}_{n}}})-\frac{1}{N}\sum_{n}{\log{{\hat{\gamma}_{n}}}^{(t+1)}}}$

30: $t=t+1$

while
$\frac{1}{L}\sum_{l=1}^{L}(||\hat{\mathbf{x}}^{l(t+1)}-\hat{\mathbf{x}}^{l(t)}||^{2}/||\hat{\mathbf{x}}^{l(t+1)}||^{2})>\delta_{x}$
and $t<t_{max}$

##  VI Numerical results

In this section, we compare the proposed UAMP-(T)SBL algorithms with the
conventional SBL and state-of-the-art AMP-based SBL algorithms. We evaluate
the performance of various algorithms using normalized MSE, defined as

NMSE $\displaystyle\triangleq$
$\displaystyle\frac{1}{{\color[rgb]{0,0,0}K}}\sum_{k=1}^{K}{||\hat{\mathbf{x}}_{k}-{\mathbf{x}_{k}}||^{2}}/{||{\mathbf{x}_{k}}||^{2}},$
(43) NMSE $\displaystyle\triangleq$
$\displaystyle\frac{1}{{\color[rgb]{0,0,0}K}{L}}\sum_{k=1}^{J}\sum_{l=1}^{L}{||\hat{\mathbf{x}}^{(l)}_{k}-{\mathbf{x}}^{(l)}_{k}||^{2}}/{||{\mathbf{x}}^{(l)}_{k}||^{2}}$
(44)

for the SMV and MMV cases respectively, where
$\hat{\mathbf{x}}_{k}~{}(\hat{\mathbf{x}}_{k}^{(l)})$ is the estimate of
${\mathbf{x}}_{k}~{}({\mathbf{x}}_{k}^{(l)})$, and ${\color[rgb]{0,0,0}K}$ is
the number of trials. Since different algorithms have different computational
complexity per iteration and they require a different number of iterations to
converge, as in [19], we measure the runtime of the algorithms to indicate
their relative computational complexity. It is noted that the time consumed by
the SVD in UAMP-SBL is counted for the runtime.

Figure 9: Performance comparison (ill-conditioned matrices). Figure 10:
Performance comparison (correlated matrices).

To test the robustness and performance of the algorithms, we use the following
measurement matrices:

  1. 1.

Ill-conditioned Matrix: Matrix $\mathbf{A}$ is constructed based on the SVD
$\mathbf{A}=\mathbf{U\Lambda V}$ where $\bm{\Lambda}$ is a singular value
matrix with $\Lambda_{i,i}/\Lambda_{i+1,i+1}=\kappa^{1/(M-1)}$ for
$i=1,2,...,M-1$ (i.e., the condition number of the matrix is $\kappa$).

  2. 2.

Correlated Matrix: The correlated matrix A is constructed using
$\textbf{A}=\textbf{C}_{L}^{1/2}\textbf{G}\textbf{C}_{R}^{1/2}$, where G is an
i.i.d. Gaussian matrix with mean zero and unit variance, and $\textbf{C}_{L}$
is an $M\times M$ matrix with the $(m,n)$th element given by $c^{|m-n|}$ where
$c\in[0,1]$. Matrix $\textbf{C}_{R}$ is generated in the same way but with a
size of $N\times N$. The parameter $c$ controls the correlation of matrix A.

  3. 3.

Non-zero Mean Matrix: The elements of matrix $\mathbf{A}$ are drawn from a
non-zero mean Gaussian distribution, i.e.,
$a_{m,n}\sim\mathcal{N}(a_{m,n}|\mu,1)$. The mean $\mu$ measures the
derivation from the i.i.d. zero-mean Gaussian matrix.

  4. 4.

Low Rank Matrix: The measurement matrix $\textbf{A}=\textbf{BC}$, where the
size of B and C are $M\times R$ and $R\times N$, respectively, and $R<M$. Both
B and C are i.i.d. Gaussian matrices with mean zero and unit variance. The
rank ratio $R/N$ is used to measure the deviation of matrix A from the i.i.d.
Gaussian matrix.

###  VI-A Numerical Results for SMV

In this section, we compare UAMP-SBL against the conventional SBL [2] and the
state-of-the-art AMP based SBL algorithm GGAMP-SBL [19] with estimated noise
variance and 3 times of the true noise variance. The vector $\mathbf{x}$ is
drawn from a Bernoulli-Gaussian distribution with a non-zero probability
$\rho$. The SNR is defined as $\text{SNR}\triangleq
E\left\|\mathbf{\mathbf{Ax}}\right\|^{2}/E\left\|\mathbf{\mathbf{w}}\right\|^{2}$.
As a performance benchmark, the support-oracle MMSE bound [19] is also
included. We set $M=800$, $N=1000$ and the SNR is set to be 60dB, unless it is
specified. For UAMP-SBL we set the maximum iteration number $t_{max}=300$
(note that there is no inner iteration in UAMP-SBL). GGAMP-SBL is a double
loop algorithm, the maximum numbers of E-step and outer iteration are set to
be 50 and 1000 respectively. The damping factor for GGAMP-SBL is 0.2 to
enhance its robustness against tough measurement matrices. It is noted that
the damping factor can be increased to reduce the runtime of GGAMP-SBL but at
the cost of reduced robustness.

Figure 11: Performance comparison: (a)low rank matrices; (b) non-zero mean
matrices.

In Fig. 9, the performance of various algorithms in terms of NMSE versus the
condition number is shown in (a) for a sparsity rate of $\rho=0.1$ and (b) for
a sparsity rate of $\rho=0.3$. It can be seen from Fig. 9(a) that UAMP-SBL
delivers the best performance (even better than the conventional SBL
algorithm), which closely approaches the support-oracle bound. With a larger
sparsity rate in Fig. 9(b), UAMP-SBL still exhibits excellent performance and
it performs slightly better than SBL and significantly better than GGAMP-SBL
when the condition number is relatively large. In addition, the simulation
performance of UAMP-SBL matches well with the performance predicted with SE.

Fig. 10 shows the performance of various algorithms versus a range of
correlation parameter $c$ from $0.1$ to $0.5$, where the sparsity rate
$\rho=0.1$ in (a) and $\rho=0.3$ in (b). It can be seen that, UAMP-SBL still
delivers exceptional performance, which is better than SBL and significantly
better than GGAMP-SBL when the correlation parameter $c$ is relatively large.
The gap between UAMP-SBL and GGAMP-SBL becomes more notable with a higher
sparsity rate. The performance of UAMP-SBL matches well with SE again.

In Fig. 11, we examine the performance of the algorithms versus rank ratio in
(a), where the sparsity rate $\rho=0.1$, and versus non-zero mean in (b),
where the sparsity rate $\rho=0.3$. It can be seen that UAMP-SBL still
delivers performance which closely matches the support-oracle bound, and is
slightly better than that of SBL. We can also see that GGAMP-SBL diverges when
the mean $\mu$ is relatively large. The performance of UAMP-SBL matches well
with SE as well.

Figure 12: Support recovery rate comparison: (a) low rank matrices; (b) non-
zero mean matrices.

In Fig. 12, we evaluate the support recovery rate of the algorithms versus
correlation parameter $c$ for correlation matrices in (a) and mean value $\mu$
for non-zero mean matrices in (b), where the sparse rate $\rho=0.3$. The
support recovery rate is defined as the percentage of successful trials in the
total trials [34]. In the noiseless case, a successful trial is recorded if
the indexes of estimated non-zero signal elements are the same as the true
indexes. In the noisy case, as the true sparse vector cannot be recovered
exactly, the recovery is regarded to be successful if the indexes of the
estimated elements with the $\mathcal{K}$ largest absolute values are the same
as the true indexes of non-zero elements in the sparse vector $\mathbf{x}$,
where $\mathcal{K}$ is the number of non-zero elements in $\mathbf{x}$. From
the results, we can see that UAMP-SBL and SBL deliver similar performance and
they can significantly outperform GGAMP-SBL when $c$ or $\mu$ is relatively
large.

The average runtime of various algorithms is shown in Figs. 13, where the
sparsity rate $\rho=0.3$, and the measurement matrice are correlated in (a)
and ill-conditioned in (b). It can be seen that UAMP-SBL is much faster than
GGAMP-SBL and SBL. SBL is normally the slowest as it has the highest
complexity due to the matrix inverse in each iteration. It is noted that, for
GGAMP-SBL, we set the damping factor to be relatively small value 0.2 to
enable it to achieve better performance and robustness. If the damping factor
is increased, GGAMP-SBL could become faster but at the cost of offsetting its
performance and robustness.

Figure 13: Runtime comparison: (a) correlated matrix; (b) ill-conditioned
matrix.

We also compare the performance of various algorithms at SNR = 35dB, and the
NMSE performance and runtime of the algorithms are shown in Fig. 14, where (a)
and (b) are for non-zero mean matrices, and (c) and (d) are for ill-
conditioned matrices. The sparsity rate $\rho=0.1$. Again, we can see that,
compared to GGAMP-SBL, UAMP-SBL delivers better performance with considerably
much smaller runtime when the mean or condition number of the matrices are
relatively large.

Figure 14: Performance and runtime comparisons of various algorithms where SNR
= 35dB.

The key difference between AMP and UAMP is that a unitary transformation is
performed in UAMP, which makes UAMP much more robust against a generic
measurement matrix. Inspired by this, we test the impact of the unitary
transformation on the GGAMP-SBL algorithm, where we first perform the unitary
transformation to the original model and then carry out GGAMP-SBL. We call
this algorithm UT-GGAMP-SBL, and compare it with UAMP-SBL. The performance and
the corresponding runtime are shown in Fig. 15, where (a) and (b) for
correlated matrices, and (c) and (d) for non-zero mean matrices. It can be
seen from this figure that, thanks to the unitary transformation, the
stability of GGAMP-SBL is significantly improved as expected. UT-GGAMP-SBL
with 3 times true noise variance achieves almost the same performance as UAMP-
SBL, however, UT-GGAMP-SBL requires the knowledge of noise variance and it is
significantly slower than UAMP-SBL.

Figure 15: Performance and runtime comparisons of UAMP-SBL and UT-GGAMP-SBL.
Figure 16: Performance comparison of various algorithms in the case of MMV.

###  VI-B Numerical Results for MMV

The elements of the sparse vectors $\left\\{\mathbf{x}^{(l)},l=1:L\right\\}$
are drawn from a Bernoulli-Gaussian distribution, and the vectors share a
common support. The number of measurement vectors is $5$. The performance of
the algorithms with ill-conditioned, correlated, non-zero mean and low-rank
measurement matrices is shown in Fig. 16 (a)-(d), respectively. In this
figure, we also include the performance of the direct extension of the
conventional SBL algorithm to the MMV model (MSBL) [35] and support-oracle
bound. It can be seen from this figure that, when the deviation of the
measurement matrices from the i.i.d. zero-mean Gaussian matrix is small,
GGAMP-SBL (with $3\times$ true noise variance) and UAMP-SBL deliver similar
performance, and both of them can approach the bound closely. MSBL works
slightly worse than GGAMP-SBL and UAMP-SBL. However, when the deviation is
relatively large, MSBL delivers slightly better performance but at high
complexity. In most cases, UAMP-SBL and MSBL almost have the same performance,
and can significantly outperform GGAMP-SBL. As an example, we show the average
runtime of different algorithms in the case of ill-conditioned matrices in
Fig. 16(e), where UAMP-SBL converges significantly faster than GGAMP-SBL and
MSBL.

Figure 17: Performance comparison of various algorithms in the case of MMV
with temporal correlation.

Furthermore, we present a numerical study to illustrate the performance of
UAMP-SBL when incorporating the temporal correlation. Besides the temporally
correlated SBL (TMSBL) [34] and GGAMP-SBL, we also compare the recovery
performance with a lower bound: the achievable NMSE by a support-aware Kalman
smoother (SKS) [36] with the knowledge of the support of the sparse vectors
and the true values of $\beta$, $\alpha$ and $\bm{\gamma}$. The SKS is
implemented in a more efficient way by incorporating UAMP. As examples, we use
low rank and non-zero mean measurement matrices to test their performance.
The sparsity rate $\rho=0.1$, SNR = 50dB and the temporal correlation
coefficient $\alpha=0.8$. It can be seen from Fig. 17 that, UAMP-TSBL can
approach the bound closely and outperform other algorithms significantly when
the rank ratio is relatively low and the mean is relatively high. In addition,
UAMP-TSBL is much faster.

##  VII Conclusion

In this paper, leveraging UAMP, we proposed UAMP-SBL for sparse signal
recovery with the framework of structured variational inference, which
inherits the low complexity and robustness of UAMP against a generic
measurement matrix. We demonstrated that, compared to the state-of-the-art AMP
based SBL algorithm, UAMP-SBL can achieve much better performance in terms of
robustness, speed and recovery accuracy. Future work includes rigorous
analyses of the state evolution of UAMP-SBL and the update mechanism of the
shape parameter.

##  Appendix A Derivation of UAMP-SBL with SVMP

We detail the forward and backward message passing in each subgraph of the
factor graph in Fig. 1 according to the principle of SVMP [23], [24], [26].
The notation $\mathcal{M}_{n_{a}\rightarrow n_{b}}(x)$ is used to denote a
message passed from node $n_{a}$ to node $n_{b}$, which is a function of $x$.
Note that, if a forward message computation requires backward messages, we use
the messages in previous iteration by default.

####  A-1 Message Computations in Subgraph 1

In this subgraph, we only need to compute the outgoing (forward) messages
$\\{\mathcal{M}_{\beta\rightarrow f_{r_{m}}}(\beta)\\}$, which are input to
Subgraph 2. The derivation of the message update rule is delayed in the
message computations in Subgraph 2, and is given in (55).

####  A-2 Message Computations in Subgraph 2

According to SVMP, we need to run BP in this subgraph except at the factor
nodes $\\{f_{r_{m}}\\}$ as they connect external variable nodes. Due to the
involvement of $\mathbf{\Phi}$, this is the most computational intensive part,
and we propose to use UAMP to handle it by integrating it to the message
passing process.

According to the derivation of (U)AMP using loopy BP, UAMP provides the
message from variable node $h_{m}$ to function node $f_{r_{m}}$. Due to the
Gaussian approximation in the the derivation of (U)AMP, the message is
Gaussian, i.e.,

$\mathcal{M}_{{h_{m}}\rightarrow
f_{r_{m}}}(h_{m})=\mathcal{M}_{f_{\delta_{m}}\rightarrow
h_{m}}(h_{m})=\mathcal{N}(h_{m}|p_{m},\tau_{p_{m}}),$ (45)

where the mean $p_{m}$ and the variance $\tau_{p_{m}}$ are respectively the
$m$th elements of $\bm{p}$ and $\bm{\tau}_{p}$ given in Line 2 and Line 1 of
the UAMP algorithm (Algorithm 1), which are also Line 2 and Line 1 of the
UAMP-SBL algorithm (Algorithm 2).

Following SVMP [26], the message
$\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)$ from factor node $f_{r_{m}}$
to variable node ${\beta}$ can be expressed as

$\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)\propto\exp\left\\{\left\langle{\log
f_{r_{m}}({{r_{m}|h_{m},\beta^{-1}}})}\right\rangle_{b\left({h_{m}}\right)}\right\\},\\\
$ (46)

where the belief of $h_{m}$ is given as

$b(h_{m})\propto\mathcal{M}_{{h_{m}}\rightarrow
f_{r_{m}}}(h_{m})\mathcal{M}_{f_{r_{m}}\rightarrow h_{m}}(h_{m}).$ (47)

Later we will see that $\mathcal{M}_{f_{r_{m}}\rightarrow
h_{m}}(h_{m})\propto\mathcal{N}(h_{m}|r_{m},\hat{\beta}^{-1})$ where
$\hat{\beta}^{-1}$ is an estimate of $\beta^{-1}$ (in the last iteration), and
its computation is delayed to (56). Hence $b(h_{m})$ is Gaussian according to
the property of the product of Gaussian functions, i.e.,
$b(h_{m})=\mathcal{N}(h_{m}|\hat{h}_{m},v_{h_{m}})$ with

$\displaystyle v_{h_{m}}$ $\displaystyle=$
$\displaystyle({{1}/{\tau_{p_{m}}}+\hat{\beta}})^{-1}$ (48)
$\displaystyle\hat{h}_{m}$ $\displaystyle=$ $\displaystyle
v_{h_{m}}(\hat{\beta}{r_{m}}+{p_{m}}/{\tau_{p_{m}}}).$ (49)

They can be rewritten in vector form as

$\displaystyle\mathbf{v}_{h}$ $\displaystyle=$
$\displaystyle\bm{\tau}_{p}./(\bm{1}+\hat{\beta}\bm{\tau}_{p})$ (50)
$\displaystyle\mathbf{\hat{h}}$ $\displaystyle=$
$\displaystyle(\hat{\beta}\bm{\tau}_{p}\cdot\bm{r}+\bm{p})./(\bm{1}+\hat{\beta}\bm{\tau}_{p}),$
(51)

to avoid numerical problems as $\bm{\tau}_{p}$ may contain zero elements,
which are Lines 3 and 4 of the UAMP-SBL algorithm. Then, from (46) and the
Gaussianity of $b(h_{m})$, the message
$\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)$ is

$\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)\propto\sqrt{\beta}\exp\left\\{{-\frac{\beta}{2}}{(|r_{m}-\hat{h}_{m}|^{2}+v_{h_{m}})}\right\\}.$
(52)

According to SVMP, the message from function node $f_{r_{m}}$ to variable node
$h_{m}$ is

$\displaystyle\mathcal{M}_{f_{r_{m}}\rightarrow h_{m}}(h_{m})$
$\displaystyle\propto\exp\left\\{\left\langle{\log
f_{r_{m}}({{r_{m}|h_{m},\beta^{-1}}})}\right\rangle_{b\left({\beta}\right)}\right\\}$
(53) $\displaystyle\propto\mathcal{N}(h_{m}|r_{m},\hat{\beta}^{-1}),$

where $\hat{\beta}=\left\langle\beta\right\rangle_{b(\beta)}$ with

$\displaystyle b(\beta)$ $\displaystyle\\!\\!\\!\\!\\!=\\!\\!\\!\\!\\!$
$\displaystyle\mathcal{M}_{\beta\rightarrow
f_{r_{m}}}(\beta)\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)$ (54)
$\displaystyle\\!\\!\\!\\!\\!=\\!\\!\\!\\!\\!$ $\displaystyle
f_{\beta}(\beta)\prod_{m}\mathcal{M}_{f_{r_{m}}\rightarrow\beta}(\beta)$
$\displaystyle\\!\\!\\!\\!\\!\propto\\!\\!\\!\\!\\!$
$\displaystyle{\beta}^{\frac{M}{2}-1}\exp\left\\{{-\frac{\beta}{2}}\sum_{m}{\left(|r_{m}-\hat{h}_{m}|^{2}+v_{h_{m}}\right)}\right\\},$

and

$\displaystyle\mathcal{M}_{\beta\rightarrow
f_{r_{m}}}(\beta)=f_{\beta}(\beta)\prod_{m^{\prime}\neq
m}\mathcal{M}_{f_{r_{m^{\prime}}}\rightarrow\beta}(\beta).$ (55)

It is noted that $b(\beta)$ is a Gamma distribution with the rate parameter
$\frac{1}{2}\sum_{m}\left(|r_{m}-\hat{h}_{m}|^{2}+v_{h_{m}}\right)$ and the
shape parameter $M/2$, so
$\hat{\beta}=\left\langle\beta\right\rangle_{b(\beta)}$ can be computed as

$\hat{\beta}={M}/{\sum_{m}{\left(|r_{m}-\hat{h}_{m}|^{2}+v_{h_{m}}\right)}},$
(56)

which can be rewritten in vector form shown in Line 5 of the UAMP-SBL
algorithm.

From (53), the Gaussian form of the message $\mathcal{M}_{f_{r_{m}}\rightarrow
h_{m}}(h_{m})$ suggests the following model

$r_{m}=h_{m}+w_{m},m=1,...,M,$ (57)

where $w_{m}$ is a Gaussian noise with mean 0 and variance $\hat{\beta}^{-1}$.
This fits into the forward recursion of the UAMP algorithm as if the noise
variance is known. Therefore, Lines 3 - 6 of the UAMP algorithm (Algorithm 1)
can be executed, which are Lines 6 - 9 of the UAMP-SBL algorithm. According to
the derivation of (U)AMP, UAMP produces the message
$\mathcal{M}_{{x_{n}}\rightarrow
f_{x_{n}}}(x_{n})\propto\mathcal{N}(x_{n}|q_{n},\tau_{q})$ with mean $q_{n}$
and variance $\tau_{q}$, which are given in Lines 5 and 6 of the UAMP
algorithm or Line 8 and Line 9 of the UAMP-SBL algorithm. We can see that the
UAMP algorithm is integrated.

The function nodes $\left\\{f_{x_{n}}\right\\}$ connect the external variable
node $\gamma_{n}$. According to SVMP, the outgoing message of Subgraph 2
$\mathcal{M}_{f_{x_{n}}\rightarrow{\gamma_{n}}}(\gamma_{n})$ can be expressed
as

$\displaystyle\mathcal{M}_{f_{x_{n}}\rightarrow{\gamma_{n}}}(\gamma_{n})\propto\exp\left\\{\left\langle{{\log
f_{x_{n}}({{x_{n}|0,\gamma_{n}^{-1}}})}}\right\rangle_{b\left({x_{n}}\right)}\right\\},$
(58)

where the belief $b(x_{n})\propto\mathcal{M}_{{x_{n}}\rightarrow
f_{x_{n}}}(x_{n})\mathcal{M}_{f_{x_{n}}\rightarrow x_{n}}(x_{n})$.

The message $\mathcal{M}_{f_{x_{n}}\rightarrow
x_{n}}(x_{n})\propto\mathcal{N}({{x_{n}|0,\hat{\gamma}_{n}^{-1}}})$ will be
computed in (65), where
$\hat{\gamma}_{n}=\left\langle\gamma_{n}\right\rangle_{b(\gamma_{n})}$. Then
$b(x_{n})$ turns out to be Gaussian, i.e.,
$b(x_{n})=\mathcal{N}({x_{n}}|\hat{x}_{n},\tau_{x_{n}})$ with

$\displaystyle\tau_{x_{n}}=\left({1}/{\tau_{q}}+\hat{\gamma}_{n}\right)^{-1}$
(59) $\displaystyle\hat{x}_{n}={q_{n}}/{(1+\tau_{q}\hat{\gamma}_{n})}.$ (60)

Performing the average operations to $\\{\tau_{x_{n}}\\}$ in (59) and
arranging (60) in a vector form lead to Lines 10 and 11 of the UAMP-SBL
algorithm. According to the above,

$\displaystyle\mathcal{M}_{f_{x_{n}}\rightarrow{\gamma_{n}}}(\gamma_{n})\propto\sqrt{\gamma_{n}}\exp\left\\{-\frac{\gamma_{n}}{2}{({|\hat{x}_{n}|^{2}+\tau_{x}})}\right\\},$
(61)

which is passed to Subgraph 3. This is the end of the message update in
Subgraph 2.

####  A-3 Message Computations in Subgraph 3

The message $\mathcal{M}_{f_{\gamma_{n}}\rightarrow{\gamma_{n}}}(\gamma_{n})$
from the factor node $f_{\gamma_{n}}$ to the variable node $\gamma_{n}$ is a
predefined Gamma distribution with shape parameter $\epsilon$ and rate
parameter $\eta$, i.e.,

$\mathcal{M}_{f_{\gamma_{n}}\rightarrow\gamma_{n}}(\gamma_{n})\propto{\gamma_{n}}^{{\epsilon}-1}\exp\left\\{-{\eta}{\gamma_{n}}\right\\}.$
(62)

According to SVMP, the message

$\mathcal{M}_{f_{x_{n}}\rightarrow
x_{n}}(x_{n})\propto\exp\left\\{\left\langle{{\log
f_{x}({{x_{n}|0,\gamma_{n}^{-1}}})}}\right\rangle_{b\left({\gamma_{n}}\right)}\right\\},$
(63)

where the belief of $\gamma_{n}$

$\displaystyle b(\gamma_{n})$
$\displaystyle\propto\mathcal{M}_{f_{\gamma_{n}}\rightarrow\gamma_{n}}(\gamma_{n})\mathcal{M}_{f_{x_{n}}\rightarrow\gamma_{n}}(\gamma_{n})$
(64)
$\displaystyle\propto{\gamma_{n}}^{{\epsilon}-\frac{1}{2}}\exp\left\\{-\frac{\gamma_{n}}{2}{({|\hat{x}_{n}|^{2}+\tau_{x}+2\eta})}\right\\}.$

Hence, the message

$\mathcal{M}_{f_{x_{n}}\rightarrow
x_{n}}(x_{n})\propto\mathcal{N}({x_{n}}|{0},\hat{\gamma}_{n}^{-1}),$ (65)

where

$\hat{\gamma}_{n}=\left\langle{{\gamma_{n}}}\right\rangle_{b({\gamma_{n}})}=\frac{{2{\epsilon}+1}}{{2\eta+{{|\hat{x}_{n}|^{2}+\tau_{x}}}}}.$
(66)

Here we set $\eta=0$, and $\hat{\gamma}_{n}$ is reduced to
$\frac{(2{\epsilon}+1)}{{|\hat{x}_{n}|^{2}+\tau_{x}}}$, which leads to Line 12
of the UAMP-SBL algorithm.

We propose to tune the parameter automatically with the empirical update rule
for $\epsilon$ shown in Line 13 of the UAMP-SBL algorithm. The iteration is
terminated when either the difference between two consecutive estimates of
$\mathbf{x}$ is smaller than a threshold or the iteration number reaches the
pre-set maximum value $t_{max}$.

##  Appendix B Proof of Proposition 1

When $\epsilon=0$, the iteration in terms of $\gamma_{n}$ has a simplified
closed form, i.e.,

$\displaystyle\gamma^{t+1}_{n}=g_{\epsilon_{0}}(\gamma^{t}_{n})=\frac{(\beta+\gamma^{t}_{n})^{2}}{(\beta
y_{n})^{2}+\beta+\gamma^{t}_{n}}.$ (67)

In order to find the fixed point, we need to solve the following equation

$f(\gamma_{n})={g_{\epsilon_{0}}(\gamma_{n})}-\gamma_{n}=0,$ (68)

which leads to the unique root

$\displaystyle\gamma^{\prime}_{n}=\frac{\beta}{\beta y^{2}_{n}-1}.$ (69)

If $\beta y^{2}_{n}>1$, the root $\gamma^{\prime}_{n}=\frac{\beta}{\beta
y^{2}_{n}-1}>0$. Taking the derivative of $g_{\epsilon_{0}}(\gamma_{n})$ in
(67), we have

$\displaystyle\frac{{d}}{{d}\gamma_{n}}g_{\epsilon_{0}}(\gamma_{n})$
$\displaystyle=1-\left(\frac{{\beta}^{2}y^{2}_{n}}{{\beta}^{2}y^{2}_{n}+\beta+\gamma_{n}}\right)^{2}.$
(70)

It is easy to verify that, when $\gamma_{n}>0$,
$0<\frac{{d}}{{d}\gamma_{n}}g_{\epsilon_{0}}(\gamma_{n})<1$. Thus, the unique
root $\gamma^{\prime}_{n}=\frac{\beta}{\beta y^{2}_{n}-1}$ is a stable fixed
point of the iteration. As
$0<\frac{{d}}{{d}\gamma_{n}}g_{\epsilon_{0}}(\gamma_{n})<1$ when
$\gamma_{n}>0$, with an initial value $\gamma^{(0)}_{n}>0$, $\gamma^{t}_{n}$
will converge to the stable fixed point $\gamma^{\prime}_{n}$ [37].

If $\beta y^{2}_{n}\leq 1$, the root $\gamma^{\prime}_{n}=\frac{\beta}{\beta
y^{2}_{n}-1}<0$ or $\gamma^{\prime}_{n}=+\infty$, i.e., there is no cross-
point between $y=g_{\epsilon_{0}}(\gamma_{n})$ and $y=\gamma_{n}$ when
$\gamma^{n}>0$. As $g_{\epsilon_{0}}(0)=\frac{\beta^{2}}{(\beta
y_{n})^{2}+\beta}>0$, $y=g_{\epsilon_{0}}(\gamma_{n})$ is above $y=\gamma_{n}$
for $\gamma_{n}>0$. In addition, $y=g_{\epsilon_{0}}(\gamma_{n})$ is an
increasing function for $\gamma_{n}>0$. Hence $\gamma_{n}^{t}$ goes to
$+\infty$ with the iteration.

##  Appendix C Proof of Theorem 1

With $\epsilon>0$, the derivative of $g_{\epsilon}(\gamma_{n})$ is given as

$\displaystyle\frac{dg_{\epsilon}(\gamma_{n})}{d\gamma_{n}}=(2\epsilon+1)\left(1-\left(\frac{\beta
u_{n}}{\beta u_{n}+\beta+\gamma_{n}}\right)^{2}\right),$ (71)

where $u_{n}=\beta y_{n}^{2}$. To find the fixed points of the iteration, we
let $f(\gamma_{n})={g_{\epsilon}(\gamma_{n})}-\gamma_{n}=0$, leading to

$2\epsilon\gamma^{2}_{n}-\gamma_{n}\beta(\beta
y^{2}_{n}-4\epsilon-1)+\beta^{2}(1+2\epsilon)=0.$ (72)

The two roots of (72) are given by 555An alternative form for the quadratic
formula is used, which can be deduced from the standard quadratic formula by
Vieta’s formulas.

$\gamma_{n(a)}=\frac{2\beta(1+2\epsilon)}{u_{n}-4\epsilon-1+\sqrt{u_{n}^{2}-8\epsilon
u_{n}-2u_{n}+1}},$ (73)

and

$\gamma_{n(b)}=\frac{2\beta(1+2\epsilon)}{u_{n}-4\epsilon-1-\sqrt{u_{n}^{2}-8\epsilon
u_{n}-2u_{n}+1}}.$ (74)

If $u_{n}>1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, it is not hard to
verify that $u_{n}-4\epsilon-1-\sqrt{u_{n}^{2}-8\epsilon u_{n}-2u_{n}+1}>0$,
so both roots are positive. Hence they are two fixed points of the iteration.
Next, we show that $\gamma_{n(a)}$ is a stable fixed point while
$\gamma_{n(b)}$ is an unstable one.

Plugging the root $\gamma_{n(a)}$ into (71), we have

$\displaystyle\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})\Bigg{|}_{\gamma_{n}=\gamma_{n(a)}}$
$\displaystyle=(2\epsilon+1)\left(1-\left(\frac{\beta u_{n}}{\beta
u_{n}+\beta+\gamma_{n(a)}}\right)^{2}\right).$ (75)

It is clear that the derivative is larger than 0. Verifying that
$\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})|_{\gamma_{n}=\gamma_{n(a)}}<1$
is equivalent to showing that

$l(u_{n})=(2\epsilon+1)(\beta u_{n})^{2}-2\epsilon(\beta
u_{n}+\beta+\gamma_{n(a)})^{2}$ (76)

is larger than 0. Inserting (73) into (76),

$\displaystyle{\color[rgb]{0,0,0}\frac{4\epsilon l(u_{n})}{\beta^{2}}}$
$\displaystyle=l_{1}(u_{n})+((4\epsilon+1)u_{n}-1)\sqrt{-l_{1}(u_{n})},$ (77)

where

$l_{1}(u_{n})=-(u_{n}^{2}-8\epsilon u_{n}-2u_{n}+1)<0.$ (78)

Then

$\displaystyle{\color[rgb]{0,0,0}\frac{4\epsilon
l(u_{n})}{\beta^{2}}}\\!=\\!\sqrt{-l_{1}(u_{n})}\big{(}-\sqrt{-l_{1}(u_{n})}\\!+\\!(4\epsilon
u_{n}\\!+\\!u_{n}\\!-\\!1)\big{)}.$ (79)

Because

$\displaystyle(4\epsilon
u_{n}+u_{n}-1)^{2}-{(-{\color[rgb]{0,0,0}l_{1}}(u_{n}))}=16\epsilon^{2}u_{n}^{2}+8\epsilon
u_{n}>0,$ (80)

the term in (79) $-\sqrt{-l_{1}(u_{n})}+(4\epsilon u_{n}+u_{n}-1)>0$ and we
have $l(u_{n})>0$. Therefore,
$\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})|_{\gamma_{n}=\gamma_{n(a)}}<1$,
i.e., $\gamma_{n(a)}$ is a stable fixed point. Similarly, it is not hard to
show that $l(u_{n})<0$ (i.e.,
$\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})>1$) for
$\gamma_{n}=\gamma_{n(b)}$, i.e., $\gamma_{n(b)}$ is an unstable fixed point.

Then we analyze the convergence behavior. As $\gamma_{n}>0$, the derivative
(71) is an increasing function and it is positive. In the above, it is already
shown that
$\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})|_{\gamma_{n}=\gamma_{n(a)}}<1$.
Therefore, for $\gamma_{n}\in[0,\gamma_{n(a)}]$,
$0<\frac{{d}}{{d}\gamma_{n}}g_{\epsilon}(\gamma_{n})<1.$ Thus, with an initial
$\gamma^{(0)}_{n}$ with the range, $\gamma^{t}_{n}$ converges to the stable
fixed point $\gamma_{n(a)}$ [37].

Next we consider $u_{n}<1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$. For
$u_{n}\in(1+4\epsilon-4\sqrt{\epsilon^{2}+\epsilon/2},1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2})$,
it can be verified that $u_{n}^{2}-8\epsilon u_{n}-2u_{n}+1<0$, leading to two
complex roots $\gamma_{n(a)}$ and $\gamma_{n(b)}$. If $u_{n}\leq
1+4\epsilon-4\sqrt{\epsilon^{2}+\epsilon/2}$, it can be shown that
$u_{n}^{2}-8\epsilon u_{n}-2u_{n}+1\geq 0$ and $u_{n}^{2}-8\epsilon
u_{n}-2u_{n}+1<(u_{n}-4\epsilon-1)^{2}$. Thus
$u_{n}-4\epsilon-1<-4\sqrt{\epsilon^{2}+\epsilon/2}<0$ and
$u_{n}-4\epsilon-1\pm\sqrt{u_{n}^{2}-8\epsilon u_{n}-2u_{n}+1}<0$, leading to
negative $\gamma_{n(a)}$ and $\gamma_{n(b)}$. In summary, if
$u_{n}<1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, the two roots are either
complex or negative. Hence, there is no cross-point between
$y=g_{\epsilon}(\gamma_{n})$ and $y=\gamma_{n}$ for $\gamma_{n}>0$. As
$g_{\epsilon}(0)=(2\epsilon+1)\frac{\beta^{2}}{(\beta y_{n})^{2}+\beta}>0$,
$y=g_{\epsilon}(\gamma_{n})$ is above $y=\gamma_{n}$. Meanwhile
$g_{\epsilon}(\gamma^{t}_{n})$ is an increasing function. Hence,
$\gamma^{t}_{n}$ goes to $+\infty$ with the iteration.

When $u_{n}=1+4\epsilon+4\sqrt{\epsilon^{2}+\epsilon/2}$, there is single root
$\gamma^{*}_{n}=\frac{2\beta(1+2\epsilon)}{u_{n}-1-4\epsilon}$. Plugging
$\gamma^{*}_{n}$ into (71), we
have$\frac{dg_{\epsilon}(\gamma_{n})}{d\gamma_{n}}=1$. Thus $\gamma^{*}_{n}$
is neutral fixed point [37]. Depending on the initial value
$\gamma^{(0)}_{n}$, $\gamma^{t}_{n}$ may converge to the fixed point
$\gamma^{*}_{n}$ or diverge.

##  Appendix D Derivation of UAMP-SBL for MMV

The belief $b(\beta)$ can be represented as

$\displaystyle b(\beta)$ $\displaystyle\propto$ $\displaystyle
f_{\beta}(\beta)\prod_{l,m}{\mathcal{M}_{f^{(l)}_{r_{m}}\rightarrow{\beta}}}(\beta)$
(81) $\displaystyle\propto$ $\displaystyle
1/\beta\prod_{l,m}\mathcal{N}({h}^{(l)}_{m}|{r}^{(l)}_{m},\hat{\beta}^{-1}).$

Then according to $\hat{\beta}=<\beta>_{b(\beta)}$, we have

$\displaystyle\hat{\beta}={ML}/{\sum_{m,l}{\left(|r^{(l)}_{m}-\hat{h}^{(l)}_{m}|^{2}+v^{(l)}_{h_{m}}\right)}}.$
(82)

According to the factor graph in Fig. 7, the belief
${b\left({{{\gamma}_{n}}}\right)}$ can be updated as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!b(\gamma^{(l)}_{n})\propto\mathcal{M}_{f^{(l)}_{\gamma_{n}\rightarrow\gamma^{(l)}_{n}}}(\gamma^{(l)}_{n})\mathcal{M}_{f^{(l)}_{{{x}_{n}}\rightarrow\gamma^{(l)}_{n}}}(\gamma^{(l)}_{n})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!=(\gamma^{(l)}_{n})^{{\epsilon}-1+\frac{1}{2}}\exp\left\\{-\frac{\gamma^{(l)}_{n}}{2}(2\eta+({|\hat{x}_{n}^{(l)}|^{2}+{\tau}^{(l)}_{x}}))\right\\}.$
(83)

Here, we still set $\eta=0$ and the expectation of $\gamma_{n}$ leads to

$\displaystyle\hat{\gamma}_{n}=\frac{2{\epsilon^{\prime}}+1}{(1/L)\sum^{L}_{l=1}({|\hat{{x}}_{n}^{(l)}|^{2}+{\tau}^{(l)}_{x}})},$
(84)

where $\epsilon^{\prime}=\epsilon/L$. By comparing (84) with (66), the update
of $\epsilon^{\prime}$ can be expressed as

$\displaystyle{\epsilon^{\prime}}=\frac{1}{2}\sqrt{\log(\frac{1}{N}\sum_{n}{\hat{\gamma}_{n}})-\frac{1}{N}\sum_{n}{\log{\hat{\gamma}}_{n}}}.$
(85)

##  Appendix E Derivation of UAMP-TSBL

We only derive the message passing for the graph shown in Fig. 8. The message
$\mathcal{M}_{f^{(l)}_{x_{n}}\rightarrow x^{(l)}_{n}}(x^{(l)}_{n})$ is
computed by the BP rule with the product of messages
$\\{\mathcal{M}_{f^{(l-1)}_{\delta_{m}}\rightarrow
x^{(l-1)}_{n}}(x^{(l-1)}_{n}),\forall m\\}$ defined in UAMP and message
$\\{\mathcal{M}_{f^{(l-1)}_{\delta_{m}}\rightarrow
x^{(l-1)}_{n}}(x^{(l-1)}_{n})\\}$, i.e.,

$\displaystyle\mathcal{M}_{f^{(l)}_{x_{n}}\rightarrow
x^{(l)}_{n}}(x^{(l)}_{n})$ (86)
$\displaystyle=\left\langle{f_{x^{(l)}_{n}}(x^{(l)})}\right\rangle_{\mathcal{M}_{f^{(l-1)}_{x_{n}}\rightarrow
x^{(l-1)}_{n}}\prod_{m}\mathcal{M}_{f^{(l-1)}_{\delta_{m}}\rightarrow
x^{(l-1)}_{n}}}$
$\displaystyle\propto\mathcal{N}(x^{(l)}_{n}|\xi^{(l)}_{n},\psi^{(l)}_{n}),$

which leads to Lines 1 to 6 of the UAMP-TSBL algorithm. Similarly, the message
$\mathcal{M}_{f^{(l+1)}_{x_{n}}\rightarrow x^{(l)}_{n}}(x^{(l)}_{n})$ from
factor node $f^{(l+1)}_{x_{n}}$ to variable node $x^{(l)}_{n}$ is also updated
by the BP rule

$\displaystyle\mathcal{M}_{f^{(l+1)}_{x_{n}}\rightarrow
x^{(l)}_{n}}(x^{(l)}_{n})$ (87)
$\displaystyle=\left\langle{f_{x^{(l+1)}_{n}}(x^{(l+1))}}\right\rangle_{\mathcal{M}_{f^{(l+2)}_{x_{n}}\rightarrow
x^{(l+1)}_{n}}\prod_{m}\mathcal{M}_{f^{(l+1)}_{\delta_{m}}\rightarrow
x^{(l+1)}_{n}}}$
$\displaystyle\propto\mathcal{N}(x^{(l)}_{n}|\theta^{(l)}_{n},\phi^{(l)}_{n}),$

leading to Lines 22 to 27 of the UAMP-TSBL algorithm. We compute the belief of
variable $x^{(l)}_{n}$ by

$\displaystyle b(x^{(l)}_{n})\propto\mathcal{M}_{f^{(l)}_{x_{n}}\rightarrow
x^{(l)}_{n}}\mathcal{M}_{f^{(l+1)}_{x_{n}}\rightarrow
x^{(l)}_{n}}\prod_{m}\mathcal{M}_{f^{(l)}_{\delta_{m}}\rightarrow
x^{(l)}_{n}}$ (88)
$\displaystyle\propto\mathcal{N}(x^{(l)}_{n}|\hat{x}^{(l)}_{n},\tau_{x}^{(l)})$

leading to Lines 19 to 20 of the UAMP-TSBL algorithm. With the beliefs
$b({x^{(l)}_{n}})$ and $b({x^{(l-1)}_{n}})$, the message
$\mathcal{M}_{f^{(l)}_{x_{n}}\rightarrow\gamma_{n}}(\gamma_{n})$ can be
obtained as

$\displaystyle\begin{aligned}
\mathcal{M}_{f^{(l)}_{x_{n}}\rightarrow\gamma_{n}}(\gamma_{n})=\exp\left\\{\left\langle{f^{(l)}_{x_{n}}(x^{(l)}_{n}|{\gamma_{n}})}\right\rangle_{b({x^{(l)}_{n}})b({x^{(l-1)}_{n}})}\right\\}.\\\
\end{aligned}$ (89)

Then, with the message
$\mathcal{M}_{f_{\gamma_{n}}\rightarrow\gamma_{n}}(\gamma_{n})$ in (62), the
belief
$b(\gamma_{n})\propto\mathcal{M}_{f_{\gamma_{n}}\rightarrow\gamma_{n}}(\gamma_{n})\mathcal{M}_{f_{x_{n}}\rightarrow\gamma_{n}}(\gamma_{n})$.
Then the update of $\hat{\gamma}_{n}$ can be expressed as

$\displaystyle\hat{\gamma}_{n}=L(2\epsilon^{\prime}+1)/(|\hat{\mathbf{x}}_{n}^{(1)}|^{2}+{\tau}^{(1)}_{x}+\frac{1}{\alpha^{2}}\sum^{L}_{l=2}({|\hat{\mathbf{x}}_{n}^{(l)}|^{2}+{\tau}^{(l)}_{x}})$
(90)
$\displaystyle+\frac{\alpha^{2}}{1-\alpha^{2}}\sum^{L-1}_{l=1}({|\hat{\mathbf{x}}_{n}^{(l)}|^{2}+{\tau}^{(l)}_{x}})-\frac{2\alpha}{1-\alpha^{2}}\sum^{L}_{l=2}({\hat{\mathbf{x}}_{n}^{(l)}\hat{\mathbf{x}}_{n}^{(l-1)}})).$

## References

  * [1] Y. C. Eldar and G. Kutyniok, _Compressed sensing: theory and applications_ , Cambridge Univ. Press, 2012.
  * [2] M. E. Tipping, “Sparse Bayesian learning and the relevance vector machine,” _J. Mach. Learn. Res._ , vol. 1, pp. 211–244, June. 2001.
  * [3] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. motivation and construction,” in _Inf. Theory Workshop on Inf. Theory_. IEEE, Jan. 2010, pp. 1–5.
  * [4] ——, “Message-passing algorithms for compressed sensing,” _Proc. Nat. Acad. Sci._ , vol. 106, no. 45, pp. 18 914–18 919, Nov. 2009.
  * [5] M. Al-Shoukairi and B. Rao, “Sparse Bayesian learning using approximate message passing,” in _Proc. 48th Asilomar Conf. Signals, Syst. Comput._ IEEE, Nov. 2014, pp. 1957–1961.
  * [6] J. Zhu, L. Han, and X. Meng, “An AMP-based low complexity generalized sparse Bayesian learning algorithm,” _IEEE Access_ , vol. 7, pp. 7965–7976, Dec. 2018.
  * [7] R. Tibshirani, “Regression shrinkage and selection via the lasso: a retrospective,” _J. Roy. Statist. Soc.: Ser. B (Statist. Methodol.)_ , vol. 73, no. 3, pp. 273–282, Apr. 2011.
  * [8] T. Blumensath and M. E. Davies, “Normalized iterative hard thresholding: Guaranteed stability and performance,” _IEEE J. Sel. Topics Signal Process._ , vol. 4, no. 2, pp. 298–309, Apr. 2010.
  * [9] A. Javanmard and A. Montanari, “State evolution for general approximate message passing algorithms, with applications to spatial coupling,” _Inf. Infer.: J. IMA_ , vol. 2, no. 2, pp. 115–144, Dec. 2013.
  * [10] S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” in _Proc. Int. Symp. Inf. Theory_. IEEE, July. 2011, pp. 2168–2172.
  * [11] X. Meng, S. Wu, and J. Zhu, “A unified Bayesian inference framework for generalized linear models,” _IEEE Signal Process. Lett._ , vol. 25, no. 3, pp. 398–402, Mar. 2018.
  * [12] S. Rangan, P. Schniter, A. K. Fletcher, and S. Sarkar, “On the convergence of approximate message passing with arbitrary matrices,” _IEEE Trans. Inf. Theory_ , vol. 65, no. 9, pp. 5339–5351, Sept. 2019.
  * [13] A. Manoel, F. Krzakala, E. W. Tramel, and L. Zdeborová, “Sparse estimation with the swept approximated message-passing algorithm,” _arXiv preprint arXiv:1406.4311_ , June. 2014.
  * [14] J. Vila, P. Schniter, S. Rangan, F. Krzakala, and L. Zdeborová, “Adaptive damping and mean removal for the generalized approximate message passing algorithm,” in _Proc. 40th IEEE ICASSP_. IEEE, Apr. 2015, pp. 2021–2025.
  * [15] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message passing,” in _Proc. Int. Symp. Inf. Theory_. IEEE, June. 2017, pp. 1588–1592.
  * [16] J. Ma and L. Ping, “Orthogonal AMP,” _IEEE Access_ , vol. 5, pp. 2020–2033, Jan. 2017.
  * [17] L. Liu, S. Huang, and B. M. Kurkoski, “Memory approximate message passing,” _IEEE ISIT, 2021, arXiv preprint arXiv:2012.10861_ , Dec. 2020.
  * [18] K. Takeuchi, “Convolutional approximate message-passing,” _IEEE Signal Process. Lett._ , vol. 27, pp. 416–420, Feb. 2020.
  * [19] M. Al-Shoukairi, P. Schniter, and B. D. Rao, “A GAMP-based low complexity sparse Bayesian learning algorithm,” _IEEE Trans. Signal Process._ , vol. 66, no. 2, pp. 294–308, Jan. 2018.
  * [20] Q. Guo and J. Xi, “Approximate message passing with unitary transformation,” _arXiv preprint arXiv:1504.04799_ , Apr. 2015.
  * [21] Z. Yuan, Q. Guo, and M. Luo, “Approximate message passing with unitary transformation for robust bilinear recovery,” _IEEE Trans. Signal Process._ , vol. 69, pp. 617–630, Dec. 2020.
  * [22] Q. Guo, D. D. Huang, S. Nordholm, J. Xi, and Y. Yu, “Iterative frequency domain equalization with generalized approximate message passing,” _IEEE Signal Process. Lett._ , vol. 20, no. 6, pp. 559–562, June. 2013.
  * [23] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, “An introduction to variational methods for graphical models,” _Mach. Learn._ , vol. 37, no. 2, pp. 183–233, Nov. 1999.
  * [24] J. Winn and C. M. Bishop, “Variational message passing,” _J. Mach. Learn. Res._ , vol. 6, pp. 661–694, Apr. 2005.
  * [25] E. P. Xing, M. I. Jordan, and S. Russell, “A generalized mean field algorithm for variational inference in exponential families,” _arXiv preprint arXiv:1212.2512_ , Oct. 2012.
  * [26] J. Dauwels, “On variational message passing on factor graphs,” in _Proc. Int. Symp. Inf. Theory_. IEEE, June 2007, pp. 2546–2550.
  * [27] Y. Jin and B. D. Rao, “Support recovery of sparse signals in the presence of multiple measurement vectors,” _IEEE Trans. Inf. Theory_ , vol. 59, no. 5, pp. 3139–3157, May. 2013.
  * [28] M. F. Duarte and Y. C. Eldar, “Structured compressed sensing: From theory to applications,” _IEEE Trans. Signal Process._ , vol. 59, no. 9, pp. 4053–4085, Sept. 2011.
  * [29] M. E. Davies and Y. C. Eldar, “Rank awareness in joint sparse recovery,” _IEEE Trans. Inf. Theory_ , vol. 58, no. 2, pp. 1135–1146, Feb. 2012.
  * [30] M. Luo, Q. Guo, D. Huang, and J. Xi, “Sparse Bayesian learning based on approximate message passing with unitary transformation,” in _VTS APWCS_. IEEE, Aug. 2019, pp. 1–5.
  * [31] H. Kang, J. Li, Q. Guo, and M. Martorella, “Pattern coupled sparse Bayesian learning based on UTAMP for robust high resolution ISAR imaging,” _IEEE Sensors J._ , vol. 20, no. 22, pp. 13 734–13 742, Nov. 2020.
  * [32] N. L. Pedersen, C. N. Manchón, D. Shutin, and B. H. Fleury, “Application of Bayesian hierarchical prior modeling to sparse channel estimation,” in _IEEE Int. Conf. Commun._ IEEE, June 2012, pp. 3487–3492.
  * [33] Z. Zhang and B. D. Rao, “Sparse signal recovery in the presence of correlated multiple measurement vectors,” in _Proc. 35th IEEE ICASSP_. IEEE, Mar. 2010, pp. 3986–3989.
  * [34] ——, “Sparse signal recovery with temporally correlated source vectors using sparse Bayesian learning,” _IEEE J. Sel. Topics Signal Process._ , vol. 5, no. 5, pp. 912–926, Sept. 2011.
  * [35] D. P. Wipf and B. D. Rao, “An empirical Bayesian strategy for solving the simultaneous sparse approximation problem,” _IEEE Trans. Signal Process._ , vol. 55, no. 7, pp. 3704–3716, July. 2007.
  * [36] R. L. Eubank, _A Kalman filter primer_. CRC Press, 2005.
  * [37] R. L. Devaney, _A first course in chaotic dynamical systems: theory and experiment_. Chapman and Hall/CRC, 2020\.