monthyeardate[] , exampleExample hypothesisHypothesis claimClaim Complex-
Valued Signal Recovery using the Bayesian LASSOD. Green, J. Lindbloom, and A.
Gelb

# Complex-Valued Signal Recovery using the Bayesian LASSO111This work is
partially supported by NSF DMS grant #1912685, AFOSR grant #F9550-22-1-0411,
DOE ASCR grant #DE-ACO5-000R22725, and DoD ONR MURI grant #N00014-20-1-2595.

Dylan Green Department of Mathematics, Dartmouth College, Hanover, NH 03755,
USA (, orcid.org/0000-0001-6184-0183 and , orcid.org/0000-0002-1789-2629 and ,
orcid.org/0000-0002-9219-4572<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Jonathan
Lindbloom22footnotemark: 2 Anne Gelb22footnotemark: 2

###### Abstract

Recovering complex-valued image recovery from noisy indirect data is important
in applications such as ultrasound imaging and synthetic aperture radar. While
there are many effective algorithms to recover point estimates of the
magnitude, fewer are designed to recover the phase. Quantifying uncertainty in
the estimate can also provide valuable information for real-time decision
making. This investigation therefore proposes a new Bayesian inference method
that recovers point estimates while also quantifying the uncertainty for
complex-valued signals or images given noisy and indirect observation data.
Our method is motivated by the Bayesian LASSO approach for real-valued sparse
signals, and here we demonstrate that the Bayesian LASSO can be effectively
adapted to recover complex-valued images whose magnitude is sparse in some
(e.g. the gradient) domain. Numerical examples demonstrate our algorithm’s
robustness to noise as well as its computational efficiency.

###### keywords:

Bayesian inference, Bayesian LASSO, complex-valued image recovery

62C10, 62F15, 65C20

##  1 Introduction

Recovering complex-valued images or signals222We use the terms image and
signal interchangeably throughout our manuscript. from noisy and/or under-
sampled data is important in coherent imaging applications such as synthetic
aperture radar (SAR) [9, 19], ultrasound [40], and digital holography [39].
The problem is often modeled as

(1) $\displaystyle\bm{y}=F\bm{z}+\bm{\varepsilon}.$

Here $\bm{z}=\bm{g}\odot e^{i{\bm{\phi}}}\in\mathbb{C}^{n}$ with $\odot$
indicating componentwise multiplication is the unknown signal of interest
decomposed into magnitude $\bm{g}\in\left(\mathbb{R^{+}}\right)^{n}$ and phase
$\bm{\phi}\in[-\pi,\pi)^{n}$, $\bm{y}\in\mathbb{C}^{m}$ represents the
observable data, $F\in\mathbb{C}^{m\crossproduct n}$ is the known forward
linear operator, and $\bm{\varepsilon}\in\mathbb{C}^{m}$ is centered Gaussian
noise with covariance $\sigma^{2}\mathbb{I}_{m}$ and probability density

$\displaystyle
f_{\mathcal{E}}(\bm{\varepsilon})=\frac{1}{\pi^{m}\sigma^{2m}}\exp(-\frac{1}{\sigma^{2}}\norm{\bm{\varepsilon}}_{2}^{2}).$

Many techniques seek to exclusively recover the magnitude ${\bm{g}}$, that is,
without considering the phase $\bm{\phi}$, and as such are designed to promote
sparsity in $\bm{g}$ or in some known transformation of $\bm{g}$ (e.g., its
gradient) [23, 31]. However phase information is also useful for some
applications, including coherent change detection [3, 19, 26, 38], high-
resolution imaging [36], and interferometry [14, 19]. In this regard, some
methods have been developed to recover point estimates of complex-valued
signals to promote sparsity of $\bm{g}$ or of some transformation of $\bm{g}$
to a sparse domain, see e.g. [8, 16, 30, 33, 41].

More recent approaches have incorporated uncertainty quantification (UQ) into
complex-valued signal recovery methods (see e.g. [10, 11, 13]), thereby
providing useful information for real-time decision making. These methods are
generally designed to incorporate prior sparse knowledge of the magnitude or
some linear transform of the magnitude of the signal into the recovery of the
posterior. In the case where sparsity is in some linear transform of the
magnitude, however, these methods do not _infer_ information regarding the
phase. Instead the phase is approximated as part of an optimization step,
limiting the uncertainty information available.

The method developed in this investigation addresses some of these issues and
brings a more comprehensive approach to recovering complex-valued images.
Specifically, we build on the Bayesian LASSO technique [29], which was
originally designed to recover sparse real-valued signals. While other
inference methods, such as generalized sparse Bayesian learning [15], recover
uncertainty information for the signal itself, the Bayesian LASSO method also
provides UQ for hyperparameters that describe the structure and overall
sparsity of the image. The idea in Bayesian LASSO is to recast the
$\ell_{1}$-regularized optimization problem (commonly employed in compressed
sensing applications [7, 12]) in a Bayesian framework by treating both the
data and the unknown signal as random variables, yielding the linear system

(2) $\mathcal{Y}=F\mathcal{X}+\mathcal{E}.$

The forward operator $F\in\mathbb{R}^{m\crossproduct n}$ is known, and
$\mathcal{X}\in\mathbb{R}^{n}$, $\mathcal{E},\mathcal{Y}\in\mathbb{R}^{m}$ are
random variables with respective realizations $\bm{x}$, $\bm{\varepsilon}$,
$\bm{y}=F\bm{x}+\bm{\varepsilon}$. An estimate of the full posterior density
function of the real-valued signal can then be recovered using Bayes’ theorem,
expressed as

(3) $f_{\mathcal{X}|\mathcal{Y}}(\bm{x}|\bm{y})\propto
f_{\mathcal{Y}|\mathcal{X}}(\bm{y}|\bm{x})f_{\mathcal{X}|\upeta}(\bm{x}|\eta)f_{\upeta}(\eta).$

Here $f_{\mathcal{Y}|\mathcal{X}}(\bm{y}|\bm{x})$ is the likelihood density
function determined by $F$ and assumptions on $\mathcal{E}$,
$f_{\mathcal{X}|\upeta}(\bm{x}|\eta)$ is the prior density function that
encodes a priori assumptions about the unknown, and $f_{\upeta}(\eta)$ is the
hyperprior density function on the scale parameter $\eta$ of the prior
density. For example, for sparsifying transform matrix
$L\in\mathbb{R}^{k\crossproduct n}$, $f_{\mathcal{X}|\eta}(\bm{x}|\eta)$ may
be defined as the product of Laplace probability densities yielding

(4)
$f_{\mathcal{X}|\upeta}(\bm{x}|\eta)\propto\prod_{i=1}^{k}\frac{1}{2\eta}\exp\left(-\frac{1}{\eta}|L_{i}\bm{x}|\right),$

where $L_{i}$ is the $i$th row of $L$ for $i=1,\dots,k$.

Our new method extends the Bayesian LASSO technique to complex-valued signals
whose magnitude is sparse in some domain by treating $\bm{y}$, $\bm{z}$, and
$\bm{\varepsilon}$ in (1) as realizations of respective random variables
$\mathcal{Y}$, $\mathcal{Z}$, and $\mathcal{E}$, leading to the probabilistic
forward model

(5) $\displaystyle\mathcal{Y}=F\mathcal{Z}+\mathcal{E}.$

The corresponding estimate of the full posterior density function of the
complex-valued signal can then be recovered according to

(6) $f_{\mathcal{Z}|\mathcal{Y}}(\bm{z}|\bm{y})\propto
f_{\mathcal{Y}|\mathcal{Z}}(\bm{y}|\bm{z})f_{\mathcal{Z}|\upeta}(\bm{z}|\eta)f_{\upeta}(\eta),$

where similarly $f_{\mathcal{Y}|\mathcal{Z}}(\bm{y}|\bm{z})$ is the likelihood
density function determined by $F$ and assumptions on $\mathcal{E}$,
$f_{\mathcal{Z}|\upeta}(\bm{z}|\eta)$ is the prior density function encoding a
priori assumptions about the unknown, and $f_{\upeta}(\eta)$ is the hyperprior
density function on $\eta$, the scale parameter of the prior density. We call
our resulting method the complex-valued Bayesian LASSO (CVBL). As a primary
benefit, the CVBL allows us to fully exploit the sparsity of the underlying
complex-valued signal in the recovery without sacrificing the phase, all while
quantifying the uncertainty regarding the entire complex-valued signal and the
hyperparameters that describe its structure and sparsity.

### Our contribution

Given noisy indirect observable data, we introduce a new Bayesian model that
uses _a priori_ assumptions regarding the magnitude of the underlying complex-
valued image of interest to recover its posterior distribution as well as to
quantify the uncertainty for both the magnitude and the phase of the unknown
signal. Adapting the Bayesian LASSO approach allows us to develop efficient
sampling techniques to simulate random draws from these resulting posterior
distributions.

### Paper organization

The rest of this paper is organized as follows. Section 2 details the
construction of the likelihood, prior, and hyperprior densities for our
method. Section 3 discusses the real-valued Bayesian LASSO approach for sparse
signal recovery, which we then expand to include recovery of signals that are
sparse in some transform domain. The corresponding complex-valued Bayesian
LASSO is then proposed in Section 4. Numerical experiments in Section 5
consider different forward operators $F$ in (1) as well as various signal to
noise (SNR) values, demonstrating our method’s utility and robustness to
noise. We also compare our results to those obtained using the more classical
LASSO maximum a posteriori (MAP) estimate approach [35]. We provide some
concluding remarks and ideas for future work in Section 6.

##  2 Bayesian Formulation

Since (1) is easily understood for one-dimensional problems, we develop our
method for $\bm{z}\in\mathbb{C}^{n}$ and note that higher-dimensional signals
can be readily vectorized to fit this form. Section 5 includes both one and
two-dimensional examples. Below we collect the ingredients needed for our new
method, including a description of the likelihood density function in Section
2.1, the construction of sparsity-promoting priors in Section 2.2, and a
review of commonly used sparsifying transform operators in Section 2.3.

###  2.1 The likelihood

The likelihood density function $f_{\mathcal{Y}|\mathcal{Z}}(\bm{y}|\bm{z})$
in (6) is determined from the density function of the noise present in the
system. Following [29] and adjusting accordingly for the complex-valued signal
case, here we assume that $\mathcal{E}$ follows a central complex normal
distribution with covariance $\sigma^{2}\mathbb{I}_{m}$, yielding

$f_{\mathcal{E}}(\bm{\varepsilon})=\mathcal{CN}(\bm{0},\sigma^{2}\mathbb{I}_{m})=\frac{1}{\pi^{m}\sigma^{2m}}\exp(-\frac{1}{\sigma^{2}}\norm{\bm{\varepsilon}}_{2}^{2}).$

The law of total probability provides that

(7)
$f_{\mathcal{Y}|\mathcal{Z}}(\bm{y}|\bm{z})=\int_{\mathbb{C}^{m}}f_{\mathcal{Y}|\mathcal{Z},\mathcal{E}}(\bm{y}|\bm{z},\bm{\varepsilon})f_{\mathcal{E}}(\bm{\varepsilon})\mathrm{d}\bm{\varepsilon}.$

When conditioned on $\bm{z}=\mathcal{Z}$ and $\bm{\varepsilon}=\mathcal{E}$,
Eq. 5 implies that $\mathcal{Y}$ is entirely determined according to the
density function

(8)
$f_{\mathcal{Y}|\mathcal{Z},\mathcal{E}}(\bm{y}|\bm{z},\bm{\varepsilon})=\delta(\bm{y}-F\bm{z}-\bm{\varepsilon}).$

Combining (7) and Eq. 8, the likelihood density function is therefore

(9)
$f_{\mathcal{Y}|\mathcal{Z}}(\bm{y}|\bm{z})=\int_{\mathbb{C}^{m}}\delta(\bm{y}-F\bm{z}-\bm{\varepsilon})f_{\mathcal{E}}(\bm{\varepsilon})\mathrm{d}\bm{\varepsilon}=f_{\mathcal{E}}(\bm{y}-F\bm{z})=\frac{1}{\pi^{m}\sigma^{2m}}\exp(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F\bm{z}}_{2}^{2}).$

######  Remark 1.

We note that an improper hyperprior was placed on $\sigma^{2}$ in [29]. For
simplicity here we assume that $\sigma^{2}$ is either explicitly known or
easily estimated. In general, $\sigma^{2}$ may be provided _a priori_ when the
errors present in the forward model are well understood. In other situations,
an estimate of $\sigma^{2}$ may be acquired when an area of low intensity in
$\mathcal{Z}$ is known or when multiple measurement vectors (MMV) are
available, see e.g. [18, 32, 42].

###  2.2 The prior density function

A primary goal for our new method is to ensure that the prior density function
used in Eq. 6 enables efficient computation for the complex-valued Bayesian
LASSO (CVBL). In this regard we assume that either the magnitude or some
linear transform of the magnitude of the unknown signal is sparse.

####  2.2.1 Sparse magnitude

Let $\mathcal{Z}$ be decomposed as $\mathcal{Z}=\mathcal{A}+i\mathcal{B}$,
where $\mathcal{A}\in\mathbb{R}^{n}$ and $\mathcal{B}\in\mathbb{R}^{n}$ are
mutually independent with respective realizations $\bm{a}$ and $\bm{b}$. When
$|\mathcal{Z}|$ is sparse we can simply adapt the approach from [42] used for
real-valued signals and define the conditional prior density function as

(10) $\displaystyle
f_{\mathcal{Z}|\upeta}(\bm{z}|\eta)=f_{\mathcal{A},\mathcal{B}|\upeta}(\bm{a},\bm{b}|\eta)=\frac{1}{(2\eta)^{n}}\exp(-\frac{1}{\eta}\sum_{j=1}^{n}\sqrt{a_{j}^{2}+b_{j}^{2}}),$

where $\upeta\in\mathbb{R}^{+}$ is a random variable with realization $\eta$.

####  2.2.2 Sparse transform of magnitude

When the magnitude is presumed sparse in a domain other than the imaging
domain, we first reformulate the model as

(11)
$\displaystyle\mathcal{Y}=F\mathcal{Z}+\mathcal{E}=F\left(\mathcal{G}\odot
e^{i{\Phi}}\right)+\mathcal{E}.$

Here $\mathcal{Z}=\mathcal{G}\odot e^{i\Phi}$, where
$\mathcal{G}\in\left(\mathbb{R}^{+}\right)^{n}$ and $\Phi\in[-\pi,\pi)^{n}$,
with $\mathcal{G}$ and $\Phi$ assumed to be mutually independent. The
componentwise decomposition of $\mathcal{Z}$ allows us to rewrite the
likelihood and prior density functions in (9) respectively as

$\displaystyle
f_{\mathcal{Y}|{\mathcal{G},\Phi}}(\bm{y}|\bm{g},\bm{\phi})=\frac{1}{\pi^{n}\sigma^{2n}}\exp\left(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F\left(\bm{g}\odot
e^{i\bm{\phi}}\right)}^{2}_{2}\right),\quad
f_{\mathcal{Z}|\upeta}(\bm{z}|\eta)=f_{\mathcal{G}|\upeta}(\bm{g}|\eta)f_{\Phi}(\bm{\phi}).$

Let $L\in\mathbb{R}^{k\times n}$ be the rank $n$ operator that transforms
$\mathcal{G}$ to the sparse domain.333The rank $n$ requirement is needed for
the particular computational implementation in Algorithm 1 and Algorithm 4. A
rank deficient sparse transform operator is commonly augmented by imposing
boundary constraints [21]. Denoting $\bm{1}_{\mathbb{R}^{+}}$ as the indicator
function for positive real vectors, the conditional prior probability density
for $\mathcal{G}$ is then

(12) $\displaystyle
f_{{\mathcal{G}}|{\upeta}}(\bm{g}|\eta)\propto\frac{1}{(2\eta)^{k}}\exp\left(-\frac{1}{\eta}\norm{L\bm{g}}_{1}\right)\bm{1}_{\mathbb{R}^{+}}(\bm{g}),$

where again $\upeta\in\mathbb{R}^{+}$ is a random variable with realization
$\eta$.

As there is often no prior information regarding the phase $\Phi$ in coherent
imaging systems, it is reasonable to impose the uniform prior as

(13) $f_{\Phi}(\bm{\phi})=\frac{1}{(2\pi)^{n}}\bm{1}_{[-\pi,\pi)}(\bm{\phi}),$

where $\bm{1}_{[-\pi,\pi)}$ denotes the indicator function for vectors whose
elements are in $[-\pi,\pi)$.

####  2.2.3 The hyperprior

There are, of course, various options for the hyperprior on $\upeta$ in (10)
and (12). In some cases it may be appropriate to simply use the delta density
function

(14) $\displaystyle f_{\upeta}(\eta)=\delta(\eta-\hat{\eta}),$

with point estimate $\hat{\eta}$ for $\upeta$. Such an estimate may be
available when considering signals whose sparsity is known or well-
approximated, see e.g. [32, 42]. Here we follow the original Bayesian LASSO
method [29] and place a gamma hyperprior on $\upeta^{-2}$ in order to maintain
conjugacy, yielding the probability density function444See discussion
surrounding Eq. 21 for explanation regarding using $\upeta^{-2}$ in place of
$\upeta$ in (3).

(15) $\displaystyle
f_{\upeta^{-2}}(\eta^{-2})=\frac{\delta^{r}}{\Gamma(r)}(\eta^{-2})^{r-1}\exp(-\delta\eta^{-2}),\quad
r>0,\delta>0.$

Note that choosing _shape_ parameter $r\leq 1$ ensures that the mode of
$f_{\upeta^{-2}}(\eta^{-2})$ is zero while also encouraging sparsity in the
solution. Furthermore, having _rate_ parameter $\delta\ll 1$ gives
$f_{\upeta^{-2}}(\eta^{-2})$ a fatter tail, making the hyperprior relatively
uninformative. To encourage sparsity while still allowing for small values of
$\eta^{-2}$, unless otherwise noted our numerical experiments use
hyperparameters

(16) $\displaystyle r=1,\quad\delta=10^{-3},$

with no additional tuning.

###  2.3 The sparsifying transform operator

This investigation assumes that the underlying signal is either sparse or that
its magnitude is piecewise constant, in which case the first order
differencing (TV) operator is well suited to suppress variation and noise in
smooth regions. The numerical experiment in Section 5.3 verifies that this
assumption is reasonable even when the true signal has piecewise smooth (but
not constant) magnitude. We note, however, that such an assumption is not an
inherent limitation to our new method, which can be straightforwardly adapted
to other sparsifying transform operators, such as higher order TV (HOTV) [2]
or wavelets [1, 22] as appropriate.

######  Remark 2.

When using $L\in\mathbb{R}^{k\crossproduct n}$ with $k>n$ as the sparsifying
transform, for instance when $L=[L_{1}^{T}\ L_{2}^{T}]^{T}$ is the 2D gradient
operator with $L_{1}\in\mathbb{R}^{n\crossproduct n}$ computing vertical
differences and $L_{2}\in\mathbb{R}^{n\crossproduct n}$ computing horizontal
differences, the posterior magnitude may become overly smoothed. We conjecture
that when $k\gg n$, the resulting posterior density is multimodal and has
significant mass concentrated where $\eta^{-2}$ is large and the magnitude is
relatively flat. Future investigations will explore methods to overcome this
undesirable result.

Before introducing our approach for complex-valued signal recovery in Section
4, Section 3 first discusses its real-valued analog. We use what we will call
the real-valued Bayesian LASSO (RVBL) technique [29], which was initially
developed to sample a posterior consisting of a Gaussian likelihood and
Laplace prior on the signal itself, that is, where $L$ is the identity matrix.
The technique was extended to include the anisotropic TV operator in [20]. We
now adapt these ideas to accommodate _any_ sparsifying rank $n$ linear
operator $L$ for real-valued signal recovery, which to our knowledge has not
been previously done.

##  3 Bayesian LASSO for a real-valued signal (RVBL)

The RVBL method is a blocked Gibbs sampling technique. Critical to the
approach is the equivalent representation of each element of the product in
(4) as a scale mixture of Gaussians with an exponential mixing density written
for each component $L_{i}\bm{x}$, $i=1,\dots,k$. In particular,

(17)
$\displaystyle\frac{1}{2\eta}\exp\left(-\frac{1}{\eta}|L_{i}\bm{x}|\right)=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi
s_{i}}}\exp\left(-\frac{(L_{i}\bm{x})^{2}}{2s_{i}}\right)\frac{1}{2\eta^{2}}\exp\left(-\frac{s_{i}}{2\eta^{2}}\right)\mathrm{d}s_{i}.$

Substituting (17) into (4) and denoting
$D(\bm{s})=\operatorname{diag}(\bm{s})$ yields

(18) $\displaystyle
f_{\mathcal{X}|\upeta}(\bm{x}|\eta)=\int_{(\mathbb{R}^{+})^{k}}\frac{1}{(2\pi)^{\frac{k}{2}}|D(\bm{s})|^{\frac{1}{2}}}\exp\left(-\frac{1}{2}(L\bm{x})^{T}[D(\bm{s})]^{-1}(L\bm{x})\right)\frac{1}{2^{k}\eta^{2k}}\exp\left(-\frac{\norm{\bm{s}}_{1}}{2\eta^{2}}\right)\mathrm{d}\bm{s},$

Now consider the marginal density of $\mathcal{X}$ over the joint distribution
of $\mathcal{X}$ and some random variable $\mathcal{T}^{2}\in\mathbb{R}^{n}$,
which we refer to as the scale mixture parameter, given by

(19) $\displaystyle
f_{\mathcal{X}|\upeta}(\bm{x}|\eta)=\int_{(\mathbb{R}^{+})^{k}}f_{\mathcal{X},\mathcal{T}^{2}|\upeta}(\bm{x},\bm{\tau}^{2}|\eta)\mathrm{d}\bm{\tau}^{2}=\int_{(\mathbb{R}^{+})^{k}}f_{\mathcal{X}|\mathcal{T}^{2},\upeta}(\bm{x}|\bm{\tau}^{2},\eta)f_{\mathcal{T}^{2}|\upeta}(\bm{\tau}^{2}|\eta)\mathrm{d}\bm{\tau}^{2}.$

Note that in Eq. 19 $\bm{\tau}^{2}\in\mathbb{R}^{k}$ is defined componentwise,
with

(20)
$\displaystyle\bm{\tau}^{2}=\begin{bmatrix}\tau_{1}^{2}&\tau_{2}^{2}&\cdots&\tau_{k}^{2}\end{bmatrix}^{T}.$

By making the substitution $\bm{\tau}^{2}$ for $\bm{s}$ and directly comparing
the integrand terms in (18) to those in (19), we obtain the Gaussian density
function

(21) $\displaystyle
f_{\mathcal{X}|\mathcal{T}^{2},\upeta^{-2}}(\bm{x}|\bm{\tau}^{2},\eta^{-2})=f_{\mathcal{X}|\mathcal{T}^{2}}(\bm{x}|\bm{\tau}^{2})\propto\frac{1}{|D(\bm{\tau}^{2})|^{\frac{1}{2}}}\exp\left(-\frac{1}{2}(L\bm{x})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{x})\right)$

as the prior along with the product of exponential densities

(22) $\displaystyle
f_{\mathcal{T}^{2}|\upeta^{-2}}(\bm{\tau}^{2}|\eta^{-2})=\prod_{i=1}^{k}\frac{\eta^{-2}}{2}\exp\left(-\frac{\tau_{i}^{2}\eta^{-2}}{2}\right)$

as the hyperprior in (3). Finally with $f_{\upeta^{-2}}(\eta^{-2})$ given in
(15),555For ease of presentation the rest of this manuscript uses the
hyperparameter $\upeta^{-2}$ in place of the original hyperparameter $\upeta$.
we are now able to completely characterize the hierarchical model from (9),
(15), (21), and (22) as

(23a) $\displaystyle\bm{y}|\bm{x}$
$\displaystyle\sim\mathcal{N}(F\bm{x},\sigma^{2}\mathbb{I})$ (23b)
$\displaystyle\bm{x}|\bm{\tau}^{2}$
$\displaystyle\sim\mathcal{N}(0,{(L^{T}[D(\bm{\tau}^{2})]^{-1}L)^{-1}})$ (23c)
$\displaystyle\bm{\tau}^{2}|\eta^{-2}$
$\displaystyle\sim\prod_{j=1}^{k}\frac{\eta^{-2}}{2}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2})\mathrm{d}\tau_{j}^{2}$
(23d) $\displaystyle\eta^{-2}$ $\displaystyle\sim\Gamma(r,\delta),$

which yields the corresponding posterior density function

(24)
$f_{\mathcal{X},\mathcal{T}^{2},\upeta^{-2}|\mathcal{Y}}(\bm{x},\bm{\tau}^{2},\eta^{-2}|\bm{y})\propto
f_{\mathcal{Y}|\mathcal{X}}(\bm{y}|\bm{x})f_{\mathcal{X}|\mathcal{T}^{2}}(\bm{x}|\bm{\tau}^{2})f_{\mathcal{T}^{2}|\upeta}(\bm{\tau}^{2}|\eta)f_{\upeta^{-2}}(\eta^{-2}),$

where

$\displaystyle f_{\mathcal{Y}|\mathcal{X}}(\bm{y}|\bm{x})$
$\displaystyle\propto$
$\displaystyle\frac{1}{(\sigma^{2})^{\frac{n}{2}}}\exp(-\frac{1}{2\sigma^{2}}\norm{\bm{y}-F\bm{x}}_{2}^{2})$
$\displaystyle f_{\mathcal{X}|\mathcal{T}^{2}}(\bm{x}|\bm{\tau}^{2})$
$\displaystyle\propto$
$\displaystyle\frac{1}{\sqrt{|D(\bm{\tau}^{2})|}}\exp\left(-\frac{1}{2}(L\bm{x})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{x})\right)$
$\displaystyle f_{\mathcal{T}^{2}|\upeta}(\bm{\tau}^{2}|\eta)$
$\displaystyle\propto$
$\displaystyle\prod_{j=1}^{k}\eta^{-2}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2})$
$\displaystyle f_{\upeta^{-2}}(\eta^{-2})$ $\displaystyle\propto$
$\displaystyle(\eta^{-2})^{r-1}\exp(-\delta\eta^{-2}).$

By isolating the parts of (24) that depend on $\bm{x}$ and defining

$\displaystyle G$ $\displaystyle\mathrel{\mathop{\ordinarycolon}}=$
$\displaystyle\frac{1}{\sigma^{2}}F^{T}F+L^{T}[D(\bm{\tau}^{2})]^{-1}L,$
$\displaystyle\bar{\bm{x}}$ $\displaystyle\mathrel{\mathop{\ordinarycolon}}=$
$\displaystyle\frac{1}{\sigma^{2}}G^{-1}F^{T}\bm{y},$

we obtain

(25) $\displaystyle
f_{\mathcal{X}|\mathcal{Y},\mathcal{T}^{2}}(\bm{x}|\bm{y},\bm{\tau}^{2})$
$\displaystyle\propto\exp(-\frac{1}{2}(\bm{x}-\bar{\bm{x}})^{T}G(\bm{x}-\bar{\bm{x}})).$

#### Conditionals on $\mathcal{T}_{j}^{2}$

When conditioned on all other variables including
$\mathcal{T}_{-j}^{2}=\\{\mathcal{T}_{j^{\prime}}^{2}\mathrel{\mathop{\ordinarycolon}}j^{\prime}=1,\dots,k,\
j^{\prime}\neq j\\}$, (24) gives

(26) $\displaystyle
f_{\mathcal{T}_{j}^{2}|\mathcal{T}_{-j}^{2},\mathcal{X},\mathcal{Y},\upeta^{-2}}(\tau_{j}^{2}|\tau_{-j}^{2},\bm{x},\bm{y},\eta^{-2})$
$\displaystyle\propto\left(\tau_{j}^{2}\right)^{-\frac{1}{2}}\exp(-\frac{[L\bm{x}]_{j}^{2}}{2\tau_{j}^{2}}-\frac{\tau_{j}^{2}\eta^{-2}}{2}),\quad
j=1,\dots,k.$

Making the change of variables $\nu_{j}^{2}=1/\tau_{j}^{2}$, the conditional
can alternatively be expressed as

(27) $\displaystyle
f_{\mathcal{V}_{j}^{2}|\mathcal{T}_{-j}^{2},\mathcal{X},\upeta^{-2}}(\nu_{j}^{2}|\tau_{-j}^{2},\bm{x},\eta^{-2})$
$\displaystyle\propto\left(\nu_{j}^{2}\right)^{-\frac{3}{2}}\exp(-\frac{[L\bm{x}]_{j}^{2}\nu_{j}^{2}}{2}-\frac{\eta^{-2}}{2\nu_{j}^{2}}).$

By defining mean parameter $\mu^{\prime}=\sqrt{\eta^{-2}/[L\bm{x}]_{j}^{2}}$
and shape parameter $\lambda^{\prime}=\eta^{-2}$, we observe that Eq. 27 fits
the form of the density function of inverse Gaussian distribution given by

$\displaystyle f_{\mu^{\prime},\lambda^{\prime}}(y)\propto
y^{-3/2}\exp(-\frac{\lambda^{\prime}y}{2(\mu^{\prime})^{2}}-\frac{\lambda^{\prime}}{2y})\propto
y^{-3/2}\exp(-\frac{\lambda^{\prime}(y-\mu^{\prime})^{2}}{2(\mu^{\prime})^{2}y})\quad\mu^{\prime},\lambda^{\prime}>0.$

Thus Eq. 27 can be directly sampled using a standard routine for sampling the
inverse Gaussian or Wald distribution (e.g., see [25]).

######  Remark 3.

In the case where $[L\bm{x}]_{j}\approx 0$ (where $L$ could be the identity
matrix $\mathbb{I}_{n}$), we can refrain from performing the change of
variables to obtain (27) and note that the original conditional (26)
simplifies to

(28) $\displaystyle
f_{\mathcal{T}_{j}^{2}|\mathcal{T}_{-j}^{2},\mathcal{X},\mathcal{Y},\upeta^{-2}}(\tau_{j}^{2}|\tau_{-j}^{2},\bm{x},\bm{y},\eta^{-2})$
$\displaystyle\propto\left(\tau_{j}^{2}\right)^{-\frac{1}{2}}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2}),$

corresponding to the Gamma distribution $\Gamma(1/2,\eta^{-2}/2)$. Thus
whenever $|[L\bm{x}]_{j}|<10^{-8}$ in our numerical experiments we instead
sample the conditional density of each $\mathcal{T}_{j}^{2}$ according to Eq.
28.666Due to system noise, nearly all values $|[L\bm{x}]_{j}|$ are greater
than $10^{-8}$ with typical values in $[.001,.1]$.

#### Conditional on $\upeta^{-2}$

Lastly, to sample the conditional on $\upeta^{-2}$, we note that

(29) $\displaystyle
f_{\upeta^{-2}|\mathcal{X},\mathcal{Y},\mathcal{T}^{2}}(\eta^{-2}|\bm{x},\bm{y},\bm{\tau}^{2})$
$\displaystyle\propto(\eta^{-2})^{r+k-1}\exp\left(-\eta^{-2}\left(\delta+\frac{1}{2}\sum_{j=1}^{k}\tau_{j}^{2}\right)\right),$

which due to conjugacy is the density for the Gamma distribution
$\Gamma(r+k,\delta+\frac{1}{2}\sum_{j=1}^{k}\tau_{j}^{2})$. Algorithm 1
outlines the full Gibbs sampling approach.

Algorithm 1 Real-valued Bayesian LASSO (RVBL)

0: Input Data vector $\bm{y}$, parameters $\sigma$ and $\eta$, sparsifying
operator $L\in\mathbb{R}^{k\times n}$, chain length $N_{M}$, and burn-in
length $B$.

0: Output Samples $\bm{x}^{(s-B+1)}$ for $s=B,\dots,N_{M}$.

0: Set $\bm{x}^{(0)}=F^{T}\bm{y}$ and
$\left(\tau_{j}^{2}\right)^{(0)}=\left(\eta^{-2}\right)^{(0)}=1$ for
$j=1,\dots,k$.

0: For $l=1,\dots,N_{M}$ do

  * i.

Sample $\bm{x}^{(l)}$ from
$\mathcal{X}\Big{|}\mathcal{T}^{2}=\left(\bm{\tau}^{2}\right)^{(l-1)}$

  * ii.

Sample each $\left(\tau_{j}^{2}\right)^{(l)}$ using
$\mathcal{T}_{j}^{-2}\Big{|}\mathcal{X}=\bm{x}^{(l)},\upeta^{-2}=\left(\eta^{-2}\right)^{(l-1)}$.

  * iii.

Sample $\left(\eta^{-2}\right)^{(l)}$ using
$\upeta^{-2}|\mathcal{T}^{2}=\left(\bm{\tau}^{2}\right)^{(l-1)}$.

We conclude this section with illustrative numerical example in Fig. 4 of the
RVBL method applied to a signal deblurring problem, where $F$ is a Gaussian
blurring kernel with standard deviation $\sigma_{\text{blur}}=5.0$ and with a
zero boundary condition. We assume that the data $\bm{y}$ are generated by
$\bm{y}=F\bm{x}+\bm{\varepsilon}$, where $\bm{\varepsilon}$ is Gaussian white
noise with SNR $=30\operatorname{dB}$, where SNR is defined in (51) for
complex-valued signals. The problem data and ground truth signal are shown in
Fig. 4 (a). To apply the RVBL method, we take the sparsifying transformation
$L$ to be the discrete gradient operator (50) and adopt the uninformative
hyperprior given by (15). Fig. 4 (b)-(c) show the resulting posterior mean
estimates for $\bm{x}$ and $\bm{\tau}^{2}$, as well as their 95% credible
intervals.

.28  .28  .28

Figure 1: Deblurring problem data Figure 2: Estimate of $\bm{x}$ Figure 3:
Estimate of $\bm{\tau}^{2}$ Figure 4: The RVBL method applied to a signal
deblurring problem. The posterior mean estimate for the parameter
$\upeta^{-2}$ is $\bar{\eta}^{-2}\approx 401$ and SNR = 30. in (51).

##  4 Complex-valued Bayesian LASSO (CVBL)

We now have all of the ingredients needed to compute and sample from the
complex-valued posterior density functions using the likelihood given in (9),
the priors described in (10) and (12), and the hyperprior (15). We will
describe the CVBL for two distinct cases: Section 4.1 considers the magnitude
itself to be sparse, while Section 4.2 examines the case for which the
magnitude is sparse in the some transform domain.

The scale mixture (17) used in the RVBL relies on the likelihood being a real-
valued probability density function. As such, for our algorithmic development
we will make use of the equivalency that for given $\bm{y}\in\mathbb{C}^{m}$,
$F\in\mathbb{C}^{m\crossproduct n}$, and $\bm{x}\in\mathbb{R}^{n}$, we have

(30)
$\displaystyle\norm{\bm{y}-F\bm{x}}_{2}^{2}=\norm{\tilde{\bm{y}}-\tilde{F}\bm{x}}_{2}^{2},$

where $\tilde{\bm{y}}=[\mathrm{Re}(\bm{y})^{T}\ \mathrm{Im}(\bm{y})^{T}]^{T}$
and $\tilde{F}=[\mathrm{Re}(F)^{T}\ \mathrm{Im}(F)^{T}]^{T}$.

###  4.1 CVBL for the Sparse Magnitude Case

We begin by writing $\mathcal{Z}$ as

$\displaystyle\mathcal{Z}=\mathcal{A}+i\mathcal{B},$

where $\mathcal{Z}\in\mathbb{C}^{n}$ and
$\mathcal{A},\mathcal{B}\in\mathbb{R}^{n}$, with realizations
$\bm{z}=\bm{a}+i\bm{b}$. When $|\mathcal{Z}|=|\mathcal{A}+i\mathcal{B}|$ is
assumed to be sparse, we choose to impose a $1-$norm prior on
$\sqrt{\mathcal{A}^{2}+\mathcal{B}^{2}}$ using (10). We now provide the
details for how this is accomplished.

For observations given in (1), an analogous derivation that yielded (23)
results in the hierarchical model for the complex-valued sparse signal
recovery

$\displaystyle\bm{y}|\bm{a},\bm{b}$
$\displaystyle\sim\mathcal{CN}(F(\bm{a}+i\bm{b}),\sigma^{2}\mathbb{I}),$
$\displaystyle\bm{a}|\bm{\tau}^{2}$
$\displaystyle\sim\mathcal{N}(0,D(\bm{\tau}^{-2})),$
$\displaystyle\bm{b}|\bm{\tau}^{2}$
$\displaystyle\sim\mathcal{N}(0,D(\bm{\tau}^{-2})),$
$\displaystyle\bm{\tau}^{2}|\eta^{-2}$
$\displaystyle\sim\prod_{j=1}^{k}\frac{\eta^{-2}}{2}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2})\mathrm{d}\tau_{j}^{2},$
$\displaystyle\eta^{-2}$ $\displaystyle\sim\Gamma(r,\delta),$

where $\bm{\tau}^{2}$ (defined componentwise as in (20)) is a realization of
the scale mixture parameter $\mathcal{T}^{2}$ introduced in (19). We can
choose hyperparameters $r$ and $\delta$ as in (16) to promote sparsity. The
posterior is then written as

$f_{\mathcal{Z},\mathcal{T}^{2},\upeta^{-2}|\mathcal{Y}}(\bm{z},\bm{\tau}^{2},\eta^{-2}|\bm{y})=f_{\mathcal{A},\mathcal{B},\mathcal{T}^{2},\upeta^{-2}|\mathcal{Y}}(\bm{a},\bm{b},\bm{\tau}^{2},\eta^{-2}|\bm{y}),$

with

(31)
$f_{\mathcal{A},\mathcal{B},\mathcal{T}^{2},\upeta^{-2}|\mathcal{Y}}(\bm{a},\bm{b},\bm{\tau}^{2},\eta^{-2}|\bm{y})\propto
f_{\mathcal{Y}|\mathcal{A},\mathcal{B}}(\bm{y}|\bm{a},\bm{b})f_{\upeta^{-2}}(\eta^{-2})\prod_{j=1}^{k}f_{\mathcal{A}_{j},\mathcal{B}_{j}|\mathcal{T}_{j}^{2}}(a_{j},b_{j}|\tau_{j}^{2})f_{\mathcal{T}_{j}^{2}|\upeta^{-2}}(\tau_{j}^{2}|\eta^{-2}).$

Here

$\displaystyle f_{\mathcal{Y}|\mathcal{A},\mathcal{B}}(\bm{y}|\bm{a},\bm{b})$
$\displaystyle\propto$
$\displaystyle\frac{1}{(\sigma^{2})^{n}}\exp(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F\bm{a}-iF\bm{b}}_{2}^{2}),$
$\displaystyle f_{\upeta^{-2}}(\eta^{-2})$ $\displaystyle\propto$
$\displaystyle(\eta^{-2})^{r-1}\exp(-\delta\eta^{-2}),$ $\displaystyle
f_{\mathcal{A}_{j},\mathcal{B}_{j}|\mathcal{T}_{j}^{2}}(a_{j},b_{j}|\tau_{j}^{2})$
$\displaystyle\propto$
$\displaystyle\frac{1}{\sqrt{\tau_{j}^{2}}}\exp(-\frac{a_{j}^{2}+b_{j}^{2}}{2\tau_{j}^{2}}),$
$\displaystyle f_{\mathcal{T}_{j}^{2}|\upeta^{-2}}(\tau_{j}^{2}|\eta^{-2})$
$\displaystyle\propto$
$\displaystyle\eta^{-2}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2}).$

For ease of notation we define

(32) $\tilde{F}=[\mathrm{Re}(F)^{T}\
\mathrm{Im}(F)^{T}]^{T},\quad\tilde{F}^{*}=[-\mathrm{Im}(F)^{T}\
\mathrm{Re}(F)^{T}]^{T},\quad\tilde{\bm{y}}=[\mathrm{Re}(\bm{y})^{T}\
\mathrm{Im}(\bm{y})^{T}]^{T},$

so that the real and imaginary parts of the observations are respectively

(33) $\bm{y}_{1}=\tilde{\bm{y}}-\tilde{F}^{*}\bm{b}\quad\text{ and
}\quad\bm{y}_{2}=\tilde{\bm{y}}-\tilde{F}\bm{a}.$

#### Sampling on $\mathcal{A}$ and $\mathcal{B}$

Analogous to the procedure resulting in (25), we now sample from the
conditionals on $\mathcal{A}$ and $\mathcal{B}$ by first isolating the parts
of (31) that depend on $\bm{a}$ and $\bm{b}$ respectively, yielding

(34a)
$f_{\mathcal{A}|\mathcal{Y},\mathcal{B},\mathcal{T}^{2}}(\bm{a}|\bm{y},\bm{b},\bm{\tau}^{2})\propto\exp(-\frac{1}{2}(\bm{a}-\bar{\bm{a}})^{T}G(\bm{a}-\bar{\bm{a}})),$
(34b)
$f_{\mathcal{B}|\mathcal{Y},\mathcal{A},\mathcal{T}^{2}}(\bm{b}|\bm{y},\bm{a},\bm{\tau}^{2})\propto\exp(-\frac{1}{2}(\bm{b}-\bar{\bm{b}})^{T}G(\bm{b}-\bar{\bm{b}})),$

where $G=\frac{2}{\sigma^{2}}\tilde{F}^{T}\tilde{F}+[D(\bm{\tau}^{2})]^{-1}$,
$\bar{\bm{a}}=\frac{2}{\sigma^{2}}G^{-1}\tilde{F}^{T}\bm{y}_{1}$ and
$\bar{\bm{b}}=\frac{2}{\sigma^{2}}G^{-1}\tilde{F}^{*T}\tilde{\bm{y}}_{2}$.
Observe that both $\mathcal{A}$ and $\mathcal{B}$ are conditionally Gaussian.

######  Remark 4.

If $F$ is unitary, for example when it represents the $n\times n$ normalized
discrete Fourier transform, then (34a) and (34b) can be simplified to
multivariate Gaussian densities with respective means

$\displaystyle\bm{a}=\frac{2}{\sigma^{2}}G^{-1}\mathrm{Re}(F^{H}(\bm{y}-iF\bm{b})),\quad\bm{b}=\frac{2}{\sigma^{2}}G^{-1}\mathrm{Im}(F^{H}(\bm{y}-F\bm{a})),$

where $G=\frac{2}{\sigma^{2}}\mathbb{I}+[D(\bm{\tau}^{2})]^{-1}.$ The
splitting in (32) is needed in applications where the input data are under-
sampled, or when portions of the data must be discarded.

#### Sampling on scale mixing parameter $\mathcal{T}_{j}^{2}$

Obtaining the conditional density of $\\{\mathcal{T}_{j}^{2}\\}_{j=1}^{k}$
from (31) yields

(35) $\displaystyle
f_{\mathcal{T}_{j}^{2}|\mathcal{A},\mathcal{B},\mathcal{T}_{-j}^{2},\upeta^{-2}}(\tau_{j}^{2}|\bm{a},\bm{b},\tau_{-j}^{2},\eta^{-2})$
$\displaystyle\propto\left(\tau_{j}^{2}\right)^{-\frac{1}{2}}\exp(-\frac{a_{j}^{2}+b_{j}^{2}}{2\tau_{j}^{2}}-\frac{\tau_{j}^{2}\eta^{-2}}{2}),\quad
j=1,\dots,k.$

By change of variable $\mathcal{V}_{j}^{2}=(\mathcal{T}_{j}^{2})^{-1}$, with
realization $\nu_{j}^{2}=(\tau_{j}^{2})^{-1}$, we obtain

(36) $\displaystyle
f_{\mathcal{V}_{j}^{2}|\mathcal{A},\mathcal{B},\mathcal{T}_{-j}^{2},\upeta^{-2}}(\nu_{j}^{2}|\bm{a},\bm{b},\tau_{-j}^{2},\eta^{-2})$
$\displaystyle\propto\left(\nu_{j}^{2}\right)^{-\frac{3}{2}}\exp(-\frac{1}{{2\nu_{j}^{2}}}\left(a_{j}^{2}+b_{j}^{2}\right)\left(\nu_{j}^{2}-\sqrt{\frac{\eta^{-2}}{(a_{j}^{2}+b_{j}^{2})}}\right)^{2}).$

Observe that comparable to (27), (36) is the probability density function for
an inverse Gaussian distribution with mean parameter
$\mu^{\prime}=\sqrt{\eta^{-2}/(a_{j}^{2}+b_{j}^{2})}$ and shape parameter
$\lambda^{\prime}=\eta^{-2}$. Moreover, since each
$(\mathcal{T}_{j}^{2})^{-1}$ is mutually independent from
$(\mathcal{T}_{-j}^{2})^{-1}$, $j=1,\dots,k$, we can efficiently sample the
conditionals on each $(\mathcal{T}_{j}^{2})^{-1}$ in parallel.

######  Remark 5.

Similar to the discussion in Remark 3, we observe that the density in (36) is
not well-defined when $a_{j}=b_{j}=0$. Although the probability of sampling
such $a_{j}$ and $b_{j}$ is zero, the probability of sampling $a_{j}$ and
$b_{j}$ within machine precision is not. Thus in the case where $a_{j}$ and
$b_{j}$ are sufficiently small, we instead use
$\mathcal{T}_{j}^{2}|\mathcal{A},\mathcal{B},\mathcal{T}_{-j}^{2},\upeta^{-2}\sim\Gamma(1/2,\eta^{-2}/2)$
following (35).

#### Sampling on hyperparameter $\upeta^{-2}$

The conditional posterior of $\upeta^{-2}$ depends exclusively on
$\mathcal{T}^{2}$, hence
$f_{\upeta^{-2}|\mathcal{A},\mathcal{B},\mathcal{Y},\mathcal{T}^{2}}=f_{\upeta^{-2}|\mathcal{T}^{2}}$.
This density then has the form of (29), from which we are able to sample
directly.

### CVBL algorithm for sparse magnitude signals

The conditional distributions in (34a), (34b), and (36) are now combined to
form the CVBL Gibbs sampling method provided in Algorithm 2. Methods to
efficiently sample $\bm{a}$ and $\bm{b}$ in Algorithm 2 are discussed in
Section 5.1.

Algorithm 2 CVBL for complex-valued signal with sparse magnitude

0: Input Data vector $\bm{y}$, noise variance $\sigma^{2}$, hyperparameters
$r=1$ and $\delta=10^{-3}$, chain length $N_{M}$, and burn-in length $B$.

0: Output samples $\bm{z}^{(s-B+1)}=\bm{a}^{(s)}+i\bm{b}^{(s)}$ for
$s=B,\dots,N_{M}$.

0: Set $\bm{a}^{(0)}=\mathrm{Re}(A^{H}\bm{y})$,
$\bm{b}^{(0)}=\mathrm{Im}(A^{H}\bm{y})$, and
$\left(\tau_{j}^{2}\right)^{(0)}=\left(\eta^{-2}\right)^{(0)}=1$ for
$j=1,\dots,k$.

0: For $\ell=1,\dots,N_{M}$ do

  * i.

Sample $\bm{a}^{(\ell)}$ from
$\mathcal{A}|\mathcal{Y}=\bm{y},\mathcal{B}=\bm{b}^{(\ell-1)},\mathcal{T}^{2}=\left(\bm{\tau}^{2}\right)^{(\ell-1)}$
(34a).

  * ii.

Sample $\bm{b}^{(\ell)}$ from
$\mathcal{B}|\mathcal{Y}=\bm{y},\mathcal{A}=\bm{a}^{(\ell)},\mathcal{T}^{2}=\left(\bm{\tau}^{2}\right)^{(\ell-1)}$
(34b).

  * iii.

Sample $\left(\tau_{j}^{-2}\right)^{(\ell)}$ from
$\mathcal{T}_{j}^{-2}|\mathcal{A}=\bm{a}^{(\ell)},\mathcal{B}=\bm{b}^{(\ell)},\upeta^{-2}=\left(\eta^{-2}\right)^{(\ell-1)}$
(36) for $j=1,\dots,n$.

  * iv.

Sample $\left(\eta^{-2}\right)^{(\ell)}$ using
$\upeta^{-2}|\mathcal{T}^{2}=\left(\bm{\tau}^{2}\right)^{(\ell)}$ (23d).

###  4.2 Sparsity in a Transform Domain of $|\mathcal{Z}|$

We now turn our attention to the case where sparsity is expected in some
transform domain of $|\mathcal{Z}|$. For simplicity we assume the signal
magnitude is piecewise constant, so that there is sparsity in the
corresponding gradient domain. Hence we define $L$ in (12) to be the TV
operator (see Section 2.3 for more discussion).

Starting from (11) and assuming $f_{\mathcal{Y}}(\bm{y})>0$, Bayes’ Theorem
yields

(37) $\displaystyle
f_{\mathcal{G},\Phi,\upeta|\mathcal{Y}}(\bm{g},\bm{\phi},\eta|\bm{y})$
$\displaystyle\propto\exp\left(-\frac{1}{\sigma^{2}}\norm{\bm{y}-A\left(\bm{g}\odot
e^{i\bm{\phi}}\right)}^{2}_{2}-\frac{1}{\eta}\norm{L\bm{g}}_{1}\right)\bm{1}_{\mathbb{R}^{+}}(\bm{g})\bm{1}_{[-\pi,\pi)}(\bm{\phi}).$

Analogously rewriting the prior as a scale mixture of normals as in (19)
provides the posterior density function

(38) $\displaystyle
f_{\mathcal{G},\Phi,\mathcal{T}^{2},\upeta^{-2}|\mathcal{Y}}(\bm{g},\bm{\phi},$
$\displaystyle\bm{\tau}^{2},\eta^{-2}|\bm{y})\propto
f_{\mathcal{Y}|\mathcal{G},\Phi}(\bm{y}|\bm{g},\bm{\phi})f_{\Phi}(\bm{\phi})f_{\mathcal{G}|\mathcal{T}^{2}}(\bm{g}|\bm{\tau}^{2})f_{\upeta^{-2}}(\eta^{-2})\prod_{j=1}^{k}f_{\mathcal{T}_{j}^{2}|\upeta^{-2}}(\bm{\tau}_{j}^{2}|\eta^{-2}),$

where

$\displaystyle f_{\mathcal{Y}|\mathcal{G},\Phi}(\bm{y}|\bm{g},\bm{\phi})$
$\displaystyle\propto$
$\displaystyle\exp(-\frac{1}{\sigma^{2}}\norm{y-F(\bm{g}\odot
e^{i\bm{\phi}})}_{2}^{2})\bm{1}_{+}(\bm{g}),$ $\displaystyle
f_{\Phi}(\bm{\phi})$ $\displaystyle\propto$
$\displaystyle\bm{1}_{[-\pi,\pi)}(\bm{\phi}),$ $\displaystyle
f_{\mathcal{G}|\mathcal{T}^{2}}(\bm{g}|\bm{\tau}^{2})$ $\displaystyle\propto$
$\displaystyle\frac{1}{\sqrt{|D(\bm{\tau}^{2})|}}\exp\left(-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})\right),$
$\displaystyle f_{\upeta^{-2}}(\eta^{-2})$ $\displaystyle\propto$
$\displaystyle(\eta^{-2})^{r-1}\exp(-\delta\eta^{-2}),$ $\displaystyle
f_{\mathcal{T}_{j}^{2}|\upeta^{-2}}(\bm{\tau}_{j}^{2}|\eta^{-2})$
$\displaystyle\propto$
$\displaystyle\eta^{-2}\exp(-\frac{\tau_{j}^{2}\eta^{-2}}{2}).$

Below we introduce a four-step Gibbs sampling process to sample from (38) that
combines (1) the use of Metropolis-within-Gibbs to sample from
$\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2},\upeta^{-2}$; (2) Gibbs steps to
sample from $\mathcal{T}^{2}|\mathcal{Y},\mathcal{G},\Phi,\upeta^{-2}$ and
$\upeta^{-2}|\mathcal{Y},\mathcal{G},\Phi,\mathcal{T}^{2}$; and (3) rejection
sampling to sample from
$\Phi|\mathcal{Y},\mathcal{G},\mathcal{T}^{2},\upeta^{-2}$.

#### Sampling on magnitude $\mathcal{G}$

The conditional distribution on $\mathcal{G}$ is a nonnegatively-constrained
Gaussian density. For computational efficiency in sampling the magnitude, we
use as the posterior density the untruncated density
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\mathcal{T}^{2},\Phi}(\bm{g}|\bm{y},\bm{\tau}^{2},\bm{\phi})$,
where

$\displaystyle
f_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}(\bm{g}|\bm{y},\bm{\phi},\bm{\tau}^{2})=\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}(\bm{g}|\bm{y},\bm{\phi},\bm{\tau}^{2})\bm{1}_{+}(\bm{g}).$

Thus

(39)
$\displaystyle\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}(\bm{g}|\bm{y},\bm{\phi},\bm{\tau}^{2})\propto\exp\left(-\frac{\norm{\bm{y}-F_{1}\bm{g}}^{2}_{2}}{\sigma^{2}}-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})\right),$

where $F_{1}=FD(e^{i\bm{\phi}})$. For consistency with (38), we simply reject
any proposed samples of $\mathcal{G}$ that contain a negative entry. This
technique can be categorized as acceptance-rejection sampling, and to this
end, Theorem 4.1 shows that (39) is a conditional Gaussian distribution,
implying that we can directly draw its samples.

######  Theorem 4.1.

Let $F_{1}=FD(e^{i\phi})$, $\tilde{F}_{1}=[\mathrm{Re}(F_{1})^{T}\
\mathrm{Im}(F_{1})^{T}]^{T}$, and $\tilde{\bm{y}}=[\mathrm{Re}(\bm{y})^{T}\
\mathrm{Im}(\bm{y})^{T}]^{T}$. Assume that $L$ has rank $n$. The function
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}(\bm{g}|\bm{y},\bm{\phi},\bm{\tau}^{2})$
in (39) defines a Gaussian density over $\mathbb{R}^{n}$ with center and
precision given by

(40a)
$\bar{\bm{g}}=\Gamma^{-1}\left(\frac{2}{\sigma^{2}}\tilde{F}_{1}^{T}\tilde{\bm{y}}\right),$
(40b)
$\Gamma=L^{T}[D(\bm{\tau}^{2})]^{-1}L+\frac{2}{\sigma^{2}}\tilde{F}_{1}^{T}\tilde{F}_{1}.$

The proof of Theorem 4.1 is provided in Appendix A.

When $F$ is unitary, the resulting Gaussian density attains a simpler form
provided by Corollary 4.2 which enables faster sampling.

######  Corollary 4.2.

Let $F_{1}=FD(e^{i\phi})$ and suppose $F$ is unitary and $L$ has rank $n$. For
$F_{1}^{H}\bm{y}=\bm{q}\odot e^{i\bm{\varphi}}$ where
$\bm{q}\in\left(\mathbb{R}^{+}\right)^{n}$ and
$\bm{\varphi}\in[-\pi,\pi)^{n}$, the function
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\mathcal{T}^{2},\Phi}(\bm{g}|\bm{y},\bm{\tau}^{2},\bm{\phi})$
in (39) defines a Gaussian density over $\mathbb{R}^{n}$ with center and
precision given by

(41a)
$\bar{\bm{g}}=\Gamma^{-1}\left(\frac{2}{\sigma^{2}}\bm{q}\odot\cos(\bm{\varphi})\right),$
(41b)
$\Gamma=L^{T}[D(\bm{\tau}^{2})]^{-1}L+\frac{2}{\sigma^{2}}\mathbb{I}_{n}.$

The proof of Corollary 4.2 can be found in Appendix A. A couple of remarks are
in order.

######  Remark 4.3.

If the density (39) has little mass in the region where $\bm{g}\geq 0$, this
rejection method may become computationally infeasible. Other techniques exist
for sampling from truncated multivariate normal distributions, such as the
exact Hamiltonian Monte Carlo method in [27] or a Gibbs sampling technique
that updates each $g_{i}$ separately for $i=1,\dots,n$. These methods are
often more computationally expensive and may become cost-prohibitive for high
dimensional problems, however.

######  Remark 4.4.

When $\bar{\bm{g}}$ in (40a) and (41a) has mostly nonnegative entries with a
few negative elements, sampling
$f_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}$ using the rejection
technique involving $\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}$
may become computationally inefficient. In this case the rejection sampling
from the mode (RSM) [24] provides a possible alternative. In short, RSM
generates an exact sample of
$f_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}$ by sampling a shifted
truncated Gaussian density followed by an acceptance-rejection step. In our
numerical experiments RSM is implemented when the rejection technique
involving (39) fails to generate a sample after $N_{s}$ attempts (in our
experiments $N_{s}=10$).

#### Updating the scale mixture parameter $\mathcal{T}^{2}$

From (38) we have the conditional posterior

(42) $\displaystyle
f_{\mathcal{T}_{j}^{2}|\mathcal{Y},\mathcal{G},\Phi,\mathcal{T}_{-j}^{2},\upeta^{-2}}(\tau_{j}^{2}|\bm{y},\bm{g},\bm{\phi},\tau_{-j}^{2},\eta^{-2})$
$\displaystyle\propto\left(\tau_{j}^{2}\right)^{-\frac{1}{2}}\exp(-\frac{[L\bm{g}]_{j}^{2}}{2\tau_{j}^{2}}-\frac{\tau_{j}^{2}\eta^{-2}}{2}),\quad{j=1,\dots,k},$

which is equivalent to the density function in the real-valued signal case
(26) with $\mathcal{X}=\mathcal{G}$. Analogously to what followed there, we
recognize (42) as an inverse Gaussian distribution with mean parameter
$\sqrt{\eta^{-2}/[L\bm{g}]_{j}^{2}}$ and shape parameter $\eta^{-2}$. The
conditional posterior on $\upeta^{-2}$ is given by (29).

#### Updating the phase $\Phi$

Lastly from (38) we have the phase posterior distribution

(43) $\displaystyle
f_{\Phi|\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\bm{\phi}|\bm{y},\bm{g},\bm{\tau}^{2})\propto\exp\left(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F_{2}e^{i\bm{\phi}}}^{2}_{2}\right)\bm{1}_{[-\pi,\pi)}(\bm{\phi}),$

where $F_{2}=FD(\bm{g})$. The nature of (43) can be better understood by
defining random variable $\Theta=e^{i\Phi}$ and corresponding realization
$\bm{\theta}=e^{i\bm{\phi}}$. The posterior is then expressed as

(44) $\displaystyle
f_{\Theta|\mathcal{Y},\mathcal{G}}(\bm{\theta}|\bm{y},\bm{g})\propto\exp\left(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F_{2}\bm{\theta}}^{2}_{2}\right)\bm{1}_{CS_{1}}(\bm{\theta}),$

where $CS_{1}$ is the unit circle in the complex plane. Clearly (44), and by
extension (43), are probability density functions of complex Gaussian
distributions restricted to the unit circle. For a general forward operator
$F$, each
$f_{\Phi_{i}|\Phi_{-i},\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\phi_{i}|\phi_{-i},\bm{y},\bm{g},\bm{\tau}^{2})$
is conditionally von Mises (45) for $i=1,\dots,n$, which we now state in
Theorem 4.5.

######  Theorem 4.5.

Let $F_{2}=FD(\bm{g})$ as in (43) and (44) and define $A=F_{2}^{H}F_{2}$ with
elements $[A]_{j,k}=a_{j,k}e^{i\alpha_{j,k}}$, where $a_{j,k}\in\mathbb{R}$
and $\alpha_{j,k}\in[-\pi,\pi)$. Further denote $F_{2}^{H}\bm{y}=\bm{q}\odot
e^{i\bm{\varphi}}$, where $\bm{q}\in\mathbb{R}_{+}^{n}$ and
$\bm{\varphi}\in[-\pi,\pi)^{n}$, with

$\displaystyle
u_{i}=\frac{2q_{i}}{\sigma^{2}}\cos\varphi_{i}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{n}\frac{2a_{i,k}}{\sigma^{2}}\cos(\tilde{\alpha}_{i,k}+\phi_{k}),\quad\quad
v_{i}=\frac{2q_{i}}{\sigma^{2}}\sin\varphi_{i}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{n}\frac{2a_{i,k}}{\sigma^{2}}\sin(\tilde{\alpha}_{i,k}+\phi_{k}),$

and $\tilde{\alpha}_{i,k}=\operatorname{sgn}(k-i)\alpha_{i,k}$ for
$k=1,\dots,n$. Then

$\displaystyle
f_{\Phi_{i}|\Phi_{-i},\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\phi_{i}|\phi_{-i},\bm{y},\bm{g},\bm{\tau}^{2})\propto
f_{vM}\left(\phi_{i}\Big{|}\mu_{i},\kappa_{i}\right),\quad i=1,\dots,n,$

where $f_{vM}(x|\mu,\kappa)$ is the von Mises probability density function,
given as

(45) $\displaystyle f_{vM}(x|\mu,\kappa)=\frac{\exp(\kappa\cos(x-\mu))}{2\pi
I_{0}(\kappa)},$

with location $\mu$ and concentration $\kappa$. Here

(46)
$\displaystyle\kappa_{i}=\sqrt{u_{i}^{2}+v_{i}^{2}},\quad\quad\mu_{i}=\begin{cases}\arctan\left(-\frac{v_{i}}{u_{i}}\right)&\text{if
}u_{i}>0\\\ \pi/2&\text{if }u_{i}=0\\\
\arctan\left(-\frac{v_{i}}{u_{i}}\right)+\pi&\text{if }u_{i}<0.\end{cases}$

Note that in (45), $I_{0}$ is the zeroth order modified Bessel function of the
first kind.

The proof to Theorem 4.5 is given in Appendix A.

Since the calculation of $\\{\kappa_{i}\\}_{i=1}^{n}$ and
$\\{\mu_{i}\\}_{i=1}^{n}$ in Eq. 46 require only scalar operations, our method
does not suffer from the curse of dimensionality as $n$ grows large.
Furthermore, when $F$ is unitary,
$f_{\Phi|\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\bm{\phi}|\bm{y},\bm{g},\bm{\tau}^{2})$
is precisely a product of von Mises density functions, as is told in Theorem
4.6.

######  Theorem 4.6.

Suppose $F_{2}=FD(\bm{g})$, where $F$ is unitary. Let
$F_{2}^{H}\bm{y}=\bm{q}\odot e^{i\bm{\varphi}}$ where
$\bm{q}\in\mathbb{R}_{+}^{n}$ and $\bm{\varphi}\in[-\pi,\pi)^{n}$. Then

$\displaystyle
f_{\Phi|\mathcal{Y},\mathcal{G}}(\bm{\phi}|\bm{y},\bm{g})=\prod_{i=1}^{n}f_{vM}\left(\phi_{i}\Big{|}\varphi_{i},\frac{2}{\sigma^{2}}q_{i}\right).$

The proof to Theorem 4.6 is provided in Appendix A.

When $F$ is unitary, the conditional independence shown in Theorem 4.6 allows
us to update each $\Phi_{j}$ independently of $\Phi_{-j}$, increasing the
opportunity for parallelization in our technique. Even when $F$ is not
unitary, Theorem 4.5 allows us to update each $\Phi_{i}$ sequentially using a
Gibbs sampling scheme, although we do not benefit from the same
parallelization opportunities as when $F$ is unitary.

To sample the von Mises distribution, we use the wrapped Cauchy distribution,
which has probability density function

$\displaystyle
f_{WC}(\theta;\mu,\gamma)=\sum_{k=-\infty}^{\infty}\frac{\gamma}{\pi(\gamma^{2}+(\theta-\mu+2\pi
k)^{2})},\quad-\pi<\theta<\pi.$

Here $\gamma$ is the scale factor and $\mu$ is the mode of the unwrapped
distribution. The acceptance-rejection method introduced in [4] utilizes a
wrapped Cauchy density as an envelope for sampling from the von Mises
distribution (45).

Algorithm 3 Sampling the von Mises density using the wrapped Cauchy

0: Input Location $\mu$ and concentration $\kappa$ of von Mises distribution.

0: Output Sample $\theta$.

0: Set $\tau=1+(1+4\kappa^{2})^{\frac{1}{2}}$,
$\rho=(\tau-(2\tau)^{\frac{1}{2}})/(2\kappa)$, and $r=(1+\rho^{2})/(2\rho)$.

0: Generate $u_{1}\sim U(0,1)$, then set $z=\cos(\pi u_{1})$,
$f=(1+rz)/(r+z)$, $c=\kappa(r-f)$.

0: Generate $u_{2}\sim U(0,1)$, then if $c(2-c)-u_{2}>0$, go to step $5$.

0: If $\ln(c/u_{2})+1-c<0$, return to step $2$.

0: Generate $u_{3}\sim U(0,1)$, then set
$\theta=[\text{sgn}(u_{3}-0.5)]\cos^{-1}(f)$.

Algorithm 4 CVBL for a signal with sparsity in the transformed magnitude

0: Input data $\bm{y}$, noise variance $\sigma^{2}$, hyperparameters $r=1$ and
$\delta=10^{-3}$, sparse transform operator $L$, chain length $N_{M}$, and
burn-in length $B$.

0: Output samples $\bm{z}^{(s-B+1)}=\bm{g}^{(s)}\odot\exp(i\bm{\phi}^{(s)})$
for $s=B,\dots,N_{M}$.

0: Set $\bm{g}^{(0)}=|A^{H}\bm{y}|$, $(\bm{\tau}^{2})^{(0)}=1$, and
$\bm{\phi}^{(0)}=\arg(A^{H}\bm{y})$

0: For $\ell=1,\dots,N_{M}$ do

  * i.

Draw a sample $\bm{g}^{*}$ from
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}\left(\bm{g}\big{|}\bar{\bm{y}},\bm{\phi}^{(\ell-1)},\left(\bm{\tau}^{2}\right)^{(\ell-1)}\right)$
using Theorem 4.1.

  * ii.

If $\bm{g}^{*}$ contains a negative element, return to (i). Otherwise, set
$\bm{g}^{(\ell)}=\bm{g}^{*}$.

  * iii.

Sample each $\left(\tau_{j}^{-2}\right)^{(\ell)}$ from
$\mathcal{T}_{j}^{-2}|\mathcal{G}=\bm{g}^{(\ell)},\upeta^{-2}=\left(\eta^{-2}\right)^{(\ell-1)}$
(42) for $j=1,\dots,k$.

  * iv.

Sample $\left(\eta^{-2}\right)^{(\ell)}$ from
$\mathcal{T}_{j}^{-2}=\left(\tau^{2}\right)^{(\ell)}$ (23d).

  * v.

For each $i=1,\dots,n$, draw a sample $\phi_{i}^{*}$ from
$\Phi_{i}|\mathcal{Y}=\bm{y},\mathcal{G}=\bm{g}^{(\ell)}$ (43) using Algorithm
3 and set $\phi_{i}^{(\ell)}=\phi_{i}^{*}$.

We are now ready to sample from the joint distribution
$f_{\mathcal{G},\Phi|\mathcal{Y},\bm{\eta}}(\bm{g},\bm{\phi}|\bm{y},\eta)$ in
(38). After initializing our chain, the three-stage Gibbs sampler is
implemented, where $\bm{g}$ is updated first, followed by the $\bm{\tau}^{2}$
update, and concluded by the $\bm{\phi}$ update. This is done for some
predetermined number of iterations $N_{M}$, after which the output chain is
formed using all the samples generated after the burn-in period $B$. This
method is summarized in Algorithm 4.777 Both Algorithm 2 and Algorithm 4 may
be easily modified to use hyperprior (14) instead of (15) by fixing the value
of $\left(\eta^{-2}\right)^{(\ell)}=\hat{\eta}^{-2}$ for all
$\ell=1,\dots,N_{M}$.

##  5 Numerical Results

We now demonstrate the efficacy of the CVBL algorithm given in Algorithm 2 for
the sparse magnitude case and Algorithm 4 for the sparse transform of the
magnitude case by performing experiments in 1D and 2D using three different
forward operators for $F$ in (1):

  * $\blacksquare$

$F_{F}\in\mathbb{C}^{n\crossproduct n}$: The discrete Fourier transform
operator with entries (in 1D) given by

(47) $\displaystyle[F_{F}]_{j,k}=\frac{1}{\sqrt{n}}\exp(-2\pi
i\frac{jk}{n}),\quad\quad j,k=0,\dots,n-1.$

Observe that $F_{F}$ is unitary, so Theorem 4.6 applies.

  * $\blacksquare$

$F_{B}\in\mathbb{R}^{n\crossproduct n}$: A blurring operator with entries (in
1D) given by

(48)
$\displaystyle[F_{B}]_{j,k}=\begin{cases}\frac{1}{\sqrt{26}}(2^{2-|j-k|})&\text{
if }|j-k|\leq 2\\\ 0&\text{ else}\end{cases},\quad\quad j,k=0,\dots,n-1.$

Observe that $F_{B}$ is a banded Toeplitz matrix, and although non-singular,
it is notoriously ill-conditioned.

  * $\blacksquare$

$F_{U}\in\mathbb{C}^{m\crossproduct n}$: A random under-sampled Fourier
transform matrix which has entries (in 1D) given by (47) but with randomly
zeroed-out rows so that $m=\lceil\nu n\rceil$, where $0<\nu\leq 1$. For
$\mathcal{M}_{\nu}\subseteq\\{2,3,\dots,n\\}$ such that
$|\mathcal{M}_{\nu}|=m$, $F_{U}$ is defined as

(49) $\displaystyle[F_{U}]_{j,k}=\begin{cases}[F_{F}]_{j,k}&\text{ if
}j\in\mathcal{M}_{\nu}\\\ 0&\text{ else}\end{cases},\quad\quad j,k=1,\dots,n.$

Note that the zeroth frequency term is never zeroed-out in the construction of
$\mathcal{M}_{\nu}$.

######  Remark 5.1.

The extension to 2D for $F_{F}$ and $F_{U}$ is straightforward. For $F_{B}$,
the 2D operator is defined such that $F_{B}\bm{z}$ convolves the image
$\bm{z}$ with the kernel $K$ given by

$\displaystyle K=\frac{1}{2\sqrt{70}}\begin{bmatrix}0&0&1&0&0\\\ 0&1&2&1&0\\\
1&2&16&2&1\\\ 0&1&2&1&0\\\ 0&0&1&0&0\end{bmatrix},$

causing the blurring effect to occur in both dimensions.888The 1D and 2D
forward operators used in our experiments are explicitly stated for
reproducibility purposes.

In the numerical examples that follow, $L=\mathbb{I}_{n}$ in the sparse signal
case or the first order differencing operator in the case of the magnitude
having a sparse gradient, where we enforce zero boundary conditions to ensure
that it is of rank $n$. For 1D signals this amounts to

(50) $\displaystyle[L]_{i,j}=\begin{cases}1&\text{if }i=j,\\\ -1&\text{if
}i=j+1,\\\ 0&\text{else}.\end{cases}$

Finally we assume that the variance $\sigma^{2}$ in our observation model (1)
is known, with corresponding signal-to-noise (SNR) given by

(51)
$\displaystyle\text{SNR}=10\log_{10}\left(\frac{\norm{F\bm{z}^{exact}}_{2}^{2}}{m\sigma^{2}}\right),$

where $\bm{z}^{exact}\in\mathbb{C}^{m}$ is the exact solution to (1). The SNR
values chosen in our experiments highlight the effectiveness of the CVBL
method in noisy environments.

In all experiments, a total of $5000$ samples are drawn with a burn-in period
of $B=200$. To encourage the sampler to move towards high-mass areas of the
probability density, for $\ell<B$ we set $\bm{g}^{(\ell)}=|\bm{g}^{*}|$,
regardless of whether or not $\bm{g}^{*}$ contains negative elements. Our 1D
analysis includes figures showing the means and $90\%$ credibility intervals
(CI) for the marginal magnitude of the posterior distribution. We also
approximate the marginal density function of the phase of random individual
pixels using a kernel density estimation technique [5]. In 2D we provide the
magnitude means and the size of the $90\%$ CIs. Finally, for each choice of
$F$ we also compare the MAP estimate of our CVBL posterior to the
corresponding classical LASSO solution $\bm{z}^{LASSO}$:

  * $\blacksquare$

sparse signal case: We compute the solution to the objective function

(52)
$\displaystyle\bm{z}^{LASSO}=\operatorname*{arg\,min}_{\bm{z}}\norm{F\bm{z}-\bm{y}}_{2}^{2}+\lambda\norm{\bm{z}}_{1}$

by considering the real and imaginary parts of $\bm{z}$, namely $\bm{a}$ and
$\bm{b}$, respectively, giving the objective function as

(53) $\displaystyle\begin{bmatrix}\bm{a}^{LASSO}\\\
\bm{b}^{LASSO}\end{bmatrix}=\operatorname*{arg\,min}_{\bm{a},\bm{b}}\left(\norm{\begin{bmatrix}F_{R}&-F_{I}\\\
F_{I}&F_{R}\end{bmatrix}\begin{bmatrix}\bm{a}\\\
\bm{b}\end{bmatrix}-\begin{bmatrix}\bm{y}_{R}\\\
\bm{y}_{I}\end{bmatrix}}_{2}^{2}+\lambda\sum_{j=1}^{n}\sqrt{a_{j}+b_{j}}\right),$

which we solve using the alternating direction method of multipliers (ADMM)
[6].

  * $\blacksquare$

sparse transform case: We again employ ADMM to compute the solution in the
sparse transform case to the generalized LASSO problem. Here, however, since
$|\bm{z}|$ is not differentiable, in order to solve

(54)
$\displaystyle\bm{z}^{LASSO}=\operatorname*{arg\,min}_{\bm{z}}\norm{F\bm{z}-\bm{y}}_{2}^{2}+\lambda\norm{L|\bm{z}|}_{1},$

we follow what was done in [11, 31] and instead use the diagonal matrix
$\Theta^{(k)}$ with non-zero entries
$\Theta^{(k)}_{jj}=z^{(k)}_{j}/|z^{(k)}_{j}|$ at each iteration $k$ of the
algorithm. The objective function then becomes

$\displaystyle\bm{z}^{LASSO}=\operatorname*{arg\,min}_{\bm{z}}\norm{F\bm{z}-\bm{y}}_{2}^{2}+\lambda\norm{L\left(\Theta^{(k)}\right)^{H}\bm{z}}_{1}.$

In both cases we test $\lambda=\alpha\sigma^{2}/\bar{\eta}$, where
$\bar{\eta}$ is the mean of the samples of $\upeta$ generated by the CVBL
method. In some sense $\alpha=1$ represents the “best case” scenario for
selecting suitable parameters, while other values of $\alpha$ allow us to test
for robustness. We use $\alpha=1$ in each experiment unless otherwise
specified.

###  5.1 Numerical Efficiency

Sampling Gaussians using matrix (e.g. Cholesky) factorization can be
prohibitively expensive in high dimensions. This is because new matrices must
be factorized for each re-sampling of the hyper-parameters $\tau^{2}$,
yielding a general cost $\mathcal{O}(n^{3})$ flops per sample. The approach
detailed in Algorithm 5 provides an efficient way to generate samples of
multivariate Gaussian distributions [28, 37], which we use to sample
$\mathcal{A}|\mathcal{Y},\mathcal{B},\mathcal{T}^{2}$ and
$\mathcal{B}|\mathcal{Y},\mathcal{A},\mathcal{T}^{2}$ in Algorithm 2 and
$\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}$ in Algorithm 4. Moreover, we
can efficiently solve the system in Step $3$ of Algorithm 5 using the
conjugate gradient method [17].

Algorithm 5 Gaussian sampling by perturbation-optimization

0: Input Mean $\bm{\mu}_{p}\in\mathbb{R}^{n}$, measurement
$\bm{y}\in\mathbb{R}^{m}$, variances $\Sigma_{p}$, $\Sigma_{\ell}$,
sparsifying transform $L\in\mathbb{R}^{k\crossproduct n}$, and forward model
$F$ corresponding to Gaussian prior
$f_{\mathcal{X}}=\mathcal{N}(L\bm{\mu}_{p},\Sigma_{p})$ and likelihood
$f_{\mathcal{Y}|\mathcal{X}}=\mathcal{N}(F\bm{x},\Sigma_{\ell})$.

0: Output Sample $\bm{x}_{s}$ of the posterior $f_{\mathcal{X}|\mathcal{Y}}$.

0: Perturb the prior mean
$\tilde{\bm{\mu}}_{p}\sim\mathcal{N}(\bm{\mu}_{p},\Sigma_{p})$.

0: Perturb the data $\tilde{\bm{y}}\sim\mathcal{N}(\bm{y},\Sigma_{\ell})$.

0: Solve
$(L^{T}\Sigma_{p}^{-1}L+F^{T}\Sigma_{\ell}^{-1}F)\bm{x}_{s}=(L^{T}\Sigma_{p}^{-1}\tilde{\bm{\mu}}_{p}+F^{T}\Sigma_{\ell}^{-1}\tilde{\bm{y}})$
with respect to $\bm{x}_{s}$.

Finally we point out that when the forward operator is given by $F_{B}$ (48)
and the magnitude sparsity is in the gradient domain, we use a block sampling
approach to generate samples of the phase (43), allowing us to takes advantage
of the inherent sparsity in $F_{B}$ and to update portions of the phase in
parallel.

###  5.2 Sparsity in Magnitude

For the sparse magnitude case, we perform experiments with Algorithm 2 using
$5000$ samples and the forward operators in (47), (48), and (49). Figure 11
shows the magnitude means and $90\%$ CIs for the 1D noisy experiments for
$\text{SNR}=20\operatorname{dB}$, (51).

.28

Figure 5: $F_{F}$ in (47)

.28

Figure 6: $F_{B}$ in (48)

.28

Figure 7: $F_{U}$ in (49); $\nu=0.8$

[b][2.3cm][t].13

.28

Figure 8: $F_{F}$ in (47)

.28

Figure 9: $F_{B}$ in (48)

.28

Figure 10: $F_{U}$ in (49); $\nu=0.8$

[b][2.1cm][t].13

Figure 11: Recovered magnitude values for complex-valued signals with sparse
magnitude. (a)-(c) The magnitude of the LASSO solution and the mean of the
CVBL output; (d)-(f) $90\%$ CI for the magnitude reconstruction. Here SNR
$=20\operatorname{dB}$.

Although the recovered magnitude means are not sparse, we see in Fig. 18 that
the real and imaginary components of the signal are close to zero outside of
the signal support. Combined with Fig. 11 our results show that, as expected,
the CVBL provides mean information similar to that of the classical LASSO
point-estimate technique (52).

.28

Figure 12: $F_{F}$ in (47)

.28

Figure 13: $F_{B}$ in (48)

.28

Figure 14: $F_{U}$ in (49); $\nu=0.8$

[b][2.3cm][t].11

.28

Figure 15: $F_{F}$ in (47)

.28

Figure 16: $F_{B}$ in (48)

.28

Figure 17: $F_{U}$ in (49); $\nu=0.8$

[b][2.3cm][t].13

Figure 18: Means of (top) real and (bottom) imaginary parts of the signal in
the sparse magnitude case. Shaded regions indicate $90\%$ CIs for the CVBL
recovery. SNR $=20\operatorname{dB}$.

Figure 22 displays the phase results for each forward operator at a randomly
selected “on” pixel (where ${\bm{g}}_{j}>0$).999The probability density plots
in Fig. 22, Fig. 33, and Fig. 75 are formed using kernel density estimation
techniques [5]. In all three experiments, the CVBL method recovers the support
of the signal as well as uncertainty information for both the magnitude and
the phase.

.21  .21  .21  [b][2.4cm][t].13

Figure 19: $F_{F}$ in (47) Figure 20: $F_{B}$ in (48) Figure 21: $F_{U}$ in
(49); $\nu=0.8$ Figure 22: Phase estimate in the sparse magnitude case at a
randomly selected “on” pixel (${\bm{g}}_{j}=1$). True phase, LASSO phase point
estimate, and CVBL marginal posterior for the phase. SNR =
$20\operatorname{dB}$.

###  5.3 Sparsity in Transform of Magnitude

We now consider the case where sparsity is expected in some linear transform
domain of the magnitude. Here the magnitude $\bm{g}=\\{g_{j}\\}_{j=1}^{200}$
is given by $\\{f(t_{j})\\}_{j=1}^{200}$ for the function
$f\mathrel{\mathop{\ordinarycolon}}[-\pi,\pi)\to\mathbb{R}_{+}$ defined as

(55) $\displaystyle f(t)=\begin{cases}2,&-2.8\leq t\leq-2.1\\\ 1.5&-1.6\leq
t\leq-1.3\\\ 1+\frac{3}{2}\exp(\left(\frac{t-\pi/2}{2/3}\right)^{2})&t>0\\\
1&\text{else}\end{cases}$

and $t_{j}=-\pi+\frac{j\pi}{100}$. The corresponding phase
$\Phi=\\{\phi_{j}\\}_{j=1}^{200}$ is randomly chosen uniformly from
$[-\pi,\pi)$. Observe that while $|f|(t)$ is sparse in the gradient domain
when $t\leq 0$, this is not the case when $t>0$. Since our prior is based on
the assumption that the gradient domain is sparse, using Eq. 55 allows us to
test the effectiveness of the CVBL method in regions where the gradient is
nonzero. As already noted in Section 2.3, other choices of sparsifying
transform operators such as HOTV may be more suitable and our method is not
inherently limited to using the identity or TV operators. For simplicity, as
well as to emphasize robustness of our approach, we use the differencing
operator in (50) and leave other operators for future investigations.

Figure 29 compares results for recovering (55) from Algorithm 4 where $5000$
samples were generated to the LASSO solution in (54) with
$\lambda=\alpha\sigma^{2}/\bar{\eta}$, where $\bar{\eta}$ is again the mean of
the smaples of $\upeta$ generated by the CVBL method and $\alpha=.1,1,10$.
Each of the three transforms, $F_{F}$ (47), $F_{B}$ (48), and $F_{U}$ (49)
with $\nu=0.8$ were considered for $SNR=20\operatorname{dB}$. The LASSO method
is clearly sensitive to the choice of regularization parameter, and is most
accurate for $\alpha=1$.

.28

Figure 23: $F_{F}$ in (47)

.28

Figure 24: $F_{B}$ in (48)

.28

Figure 25: $F_{U}$ in (49); $\nu=0.8$

[b][2.3cm][t].13

.28

Figure 26: $F_{F}$ in (47)

.28

Figure 27: $F_{B}$ in (48)

.28

Figure 28: $F_{U}$ in (49); $\nu=0.8$

[b][2.1cm][t].13

Figure 29: Magnitude recovery for ${\bm{g}}$ given by (55). (top) The CVBL
mean along with the magnitude of LASSO solutions using different
regularization parameters; (bottom) $90\%$ CI for the magnitude
reconstruction. Here SNR $=20\operatorname{dB}$.

.21  .21  .21  [b][2cm][t].13

Figure 30: $F_{F}$ in (47) Figure 31: $F_{B}$ in (48) Figure 32: $F_{U}$ in
(49); $\nu=0.8$ Figure 33: Phase estimate in the sparse transform case of a
random pixel. True phase (black dash), LASSO phase point estimate (red dash-
dot), and CVBL marginal posterior for the phase. SNR $=20\operatorname{dB}$.

Figure 33 demonstrates that the CVBL method recovers marginal posterior
density functions for the phase that are consistent with the estimates
calculated by (54). The true phase is also located in the regions of large
mass generated by the corresponding kernel density estimation. Figure 37 shows
the mean of the $\tau^{2}$ samples (42) for each experiment, which as expected
is largest in support regions of the sparse domain.

.21  .21  .21  [b][2.2cm][t].13

Figure 34: $F_{F}$ in (47) Figure 35: $F_{B}$ in (48) Figure 36: $F_{U}$ in
(49); $\nu=0.8$ Figure 37: Mean of $\tau^{2}$ (right axis) compared with the
true magnitude (left axis). SNR $=20\operatorname{dB}$.

###  5.4 Noise Study

We now analyze the effect of noise on Algorithm 4. To this end for both the
sparse magnitude case and the signal corresponding to (55) we consider
$10\operatorname{dB}\leq$ SNR $\leq 30\operatorname{dB}$.

.4  .4  [b][2.4cm][t].08

Figure 38: Sparse signal Figure 39: Sparse transform of magnitude, (55) Figure
40: Average error in the sample mean generated using Algorithm 4 for
increasing SNR. The vertical axes are in logarithmic scale.

Figure 40 demonstrates that the average error is similar for all three forward
operators across the range of SNR values, with a greater overall range of
errors with respect to the forward operator in the sparse transform case.
These results are consistent with what is observed in Fig. 29 as well as what
is apparent in Fig. 29 (b). In particular the approximation in the region
$[-1.7,-1]$ is not well resolved. More insight is provided in Fig. 37 (b),
where it is evident that the support in the sparse domain in that range is not
clearly identified.

.21  .21  .21  .21  [b][2.4cm][t].12

Figure 41: $F_{F}$ in (47) Figure 42: $F_{B}$ in (48) Figure 43: $F_{U}$ in
(49); $\nu=0.8$ Figure 44: $F_{U}$; $0.5\leq\nu\leq 1$ Figure 45: Comparison
of the average phase error in mean of samples generated using CVBL with the
LASSO solution (54) for the test function (55) at (a)-(c) increasing SNR and
(d) increasing sampling rate $0.5\leq\nu\leq 1$ with SNR
$=20\operatorname{dB}$. The vertical axes are in logarithmic scale.

Figure 45 (a)-(c) compares the average phase error over a range of SNR values
for the CVBL and generalized LASSO point estimates. While the LASSO method
outperforms Algorithm 4 for each choice of $F$ with increasing SNR, the error
difference is neglible for low SNR values for $F_{B}$ and $F_{U}$. Hence we
see that the optimal point estimate recovery algorithm essentially depends on
the SNR and sampling rate of the observable data. Uncertainty information,
however, is only acquired when using CVBL, as the generalized LASSO technique
does not infer the phase information.

Finally, Fig. 45 (d) compares the average phase error in the sample mean for
the function with magnitude given by (55) using CVBL with the LASSO solution
for the undersampled observable data case. Specifically, the data are obtained
using $F_{U}$ in (49) for a range of sampling rates $\nu$ while the SNR is
held constant at $20\operatorname{dB}$. While consistent with Fig. 45 (c),
this result also suggests that for smaller values of $\nu$ (more
undersampling) the CVBL method provides on average a lower average phase error
than the generalized LASSO.

###  5.5 2D Experiments

Our 2D experiments consider sparsity in the signal magnitude with SNR
$=0\operatorname{dB}$ and sparsity in the signal magnitude gradient with SNR
$=15\operatorname{dB}$ and SNR $=25\operatorname{dB}$.

.24  .24  .24

Figure 46: True magnitude Figure 47: MLE Figure 48: Mean posterior Figure 49:
Maximum likelihood and mean posterior point estimates of a sparse magnitude
complex-valued signal for forward operator $F_{F}$ in (47). SNR
$=0\operatorname{dB}$.

Figure 49 displays point estimates for the magnitude using $F_{F}$ while
Figure 54 shows the corresponding results for $F_{B}$ and $F_{U}$. In each
case the same $64\crossproduct 64$ sparse images is used as input and we
simulate $5000$ samples of the posterior. We observe that much of the
background noise is suppressed while the fidelity of the support of the signal
is maintained when compared to the maximum likelihood point estimate.

.24  .24  .24  .24

Figure 50: MLE; $F_{B}$ Figure 51: Mean posterior; $F_{B}$ Figure 52: MLE;
$F_{U}$ Figure 53: Mean posterior; $F_{U}$ Figure 54: Maximum likelihood and
mean posterior point estimates recoveries of a sparse magnitude complex-valued
signal. Forward operator (a-b) $F_{B}$ in (48) (c-d) $F_{U}$ $\nu=0.8$ in
(49). SNR = $0\operatorname{dB}$.

We now consider the case of promoting sparsity in the magnitude gradient for
the Shepp-Logan phantom image [34] depicted in Fig. 63 (a) on a
$256\crossproduct 256$ grid. A random phase is then added to each pixel, so
that the resulting image $\bm{z}$ is modeled by (1). Figure 63 and Fig. 70
compare the posterior means of the CVBL to the generalized LASSO estimates for
$F_{F}$ in Eq. 47 and $F_{B}$ in (48).101010The dense nature of $F_{U}$ in
(49) makes exact sampling of the phase $\bm{\phi}$ impractical for large
signals. Approximate or stochastic techniques may expedite sampling, but is
beyond the scope of the current investigation. To avoid the issue discussed in
Remark 2, in the Shepp Logan experiments we fix $\eta$ to be
$\hat{\eta}=10^{-2}$, which was chosen heuristically, with
$\lambda=\sigma^{2}/\hat{\eta}$ used for the generalized LASSO regularization
parameter. The hyperprior in turn is given by (14) rather than (15).

.23  .23  .23  .23  .23  .23  .23  .23

Figure 55: Ground truth Figure 56: Generalized LASSO Figure 57: Posterior mean
Figure 58: $90\%$ CI Figure 59: Horizontal differences Figure 60: Generalized
LASSO Figure 61: Posterior mean Figure 62: $90\%$ CI Figure 63: Magnitude
recovery of the complex-valued Shepp Logan phantom using $F_{F}$ in (47). (a)
Shepp-Logan phantom; (e) horizontal differences of the ground truth. (b)-(d)
and (f)-(h) correspond to SNR values of $15\operatorname{dB}$ and
$25\operatorname{dB}$, respectively. $\hat{\eta}=10^{-2}$.

Fig. 63 and Fig. 70 indicate that the generalized LASSO estimate (54) is
smoother when compared to the CVBL posterior mean. In particular we observe
that the three lower circles of the generalized LASSO estimate of the phantom
are “blurred” together but remain distinct for the CVBL. These results seem to
indicate that when using the CVBL method, a smaller value of $\hat{\eta}$,
that is, increasing the reliance on the prior, is needed to match the amount
of regularization apparent the generalized LASSO approach.

.23  .23  .23
.23  .23  .23

Figure 64: Generalized LASSO Figure 65: Posterior mean Figure 66: $90\%$ CI
Figure 67: Generalized LASSO Figure 68: Posterior mean Figure 69: $90\%$ CI
Figure 70: Magnitude recovery of the complex-valued Shepp Logan phantom using
$F_{B}$ in (48). (a)-(c) and (d)-(f) correspond to SNR values of
$15\operatorname{dB}$ and $25\operatorname{dB}$, respectively.
$\hat{\eta}=10^{-2}$.

.21  .21  .21  .21  [b][2.4cm][t].12

Figure 71: $F_{F}$, $15\operatorname{dB}$ Figure 72: $F_{F}$,
$25\operatorname{dB}$ Figure 73: $F_{B}$, $15\operatorname{dB}$ Figure 74:
$F_{B}$, $25\operatorname{dB}$ Figure 75: Recovered phase at a randomly chosen
pixel of the complex-valued Shepp Logan image using $F_{F}$ in (47) and
$F_{B}$ in (48) with SNR of $15\operatorname{dB}$ and $25\operatorname{dB}$.
$\hat{\eta}=10^{-2}$.

Finally, Fig. 75 demonstrates that the true phase of a randomly selected pixel
is within an area of high density of the kernel density estimation of the
marginal posterior density of the phase, and that the kernel density
estimation mode is consistent with those obtained using the generalized LASSO
technique.

##  6 Concluding remarks and future work

This investigation extends the real-valued Bayesian LASSO (RVBL), which was
originally designed to promote sparsity in a sparse signal, in two ways. We
first show that it can be modified to promote sparsity in a chosen transform
(here the gradient) domain. We then demonstrate that the RVBL can be further
expanded to include complex-valued signals and images. We call our method the
complex-valued Bayesian LASSO (CVBL). Our numerical experiments show that the
CVBL can efficiently recover samples from the entire complex-valued posterior
density function, enabling uncertainty quantification of both the magnitude
and phase of the true signal.

The CVBL is practical for coherent imaging problems with unitary or sparse
forward operators, since it is easily parallelizable. Developing surrogates
for dense forward operators will be necessary to efficiently sample large
problems, and will be the focus of future investigations. We will also
consider different sparse transform operators (such at HOTV) along with
adaptive empirical hyperparameters to avoid the pitfulls of over-
regularization. Recent work in [42] may be useful in this regard.

##  Appendix A Proofs

######  Theorem 4.1.

Let $F_{1}=FD(e^{i\phi})$, $\tilde{F}_{1}=[\mathrm{Re}(F_{1})^{T}\
\mathrm{Im}(F_{1})^{T}]^{T}$, and $\tilde{\bm{y}}=[\mathrm{Re}(\bm{y})^{T}\
\mathrm{Im}(\bm{y})^{T}]^{T}$. Assume that $L$ is rank $n$. The function
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\Phi,\mathcal{T}^{2}}(\bm{g}|\bm{y},\bm{\phi},\bm{\tau}^{2})$
in (39) defines a Gaussian density over $\mathbb{R}^{n}$ with center and
precision given by

(56a)
$\bar{\bm{g}}=\Gamma^{-1}\left(\frac{2}{\sigma^{2}}\tilde{F}_{1}^{T}\tilde{\bm{y}}\right),$
(56b)
$\Gamma=L^{T}[D(\bm{\tau}^{2})]^{-1}L+\frac{2}{\sigma^{2}}\tilde{F}_{1}^{T}\tilde{F}_{1}.$

######  Proof A.1.

Let
$\Gamma=L^{T}[D(\bm{\tau}^{2})]^{-1}L+\frac{2}{\sigma^{2}}\tilde{F}_{1}^{T}\tilde{F}_{1}$,
and suppose $\bm{x}\in\operatorname{Ker}(\Gamma)$ so that
$\bm{x}^{T}\Gamma\bm{x}=0$. By (40b) we then have

$\displaystyle\bm{x}^{T}\Gamma\bm{x}=\bm{x}^{T}L^{T}[D(\bm{\tau}^{2})]^{-1}L\bm{x}+\frac{2}{\sigma^{2}}\bm{x}^{T}\tilde{F}_{1}^{T}\tilde{F}_{1}\bm{x}=0,$

that is $\bm{x}\in\operatorname{Ker}(L)\cap\operatorname{Ker}(\tilde{F}_{1})$.
Since $\operatorname{Ker}(L)=\\{0\\}$, we have immediately that
$\operatorname{Ker}(\Gamma)=\\{0\\}$ so that $\Gamma$ is invertible. We
furthermore observe that $\Gamma$ is also symmetric positive definite. To
explicitly determine the center and precision as ${\bar{\bm{g}}}$ and $\Gamma$
in (40a) and (40b), consider the following from (39):

$\displaystyle-\frac{\norm{\bm{y}-F_{1}\bm{g}}^{2}_{2}}{\sigma^{2}}-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})$
$\displaystyle=-\frac{\norm{\tilde{\bm{y}}-\tilde{F}_{1}\bm{g}}^{2}_{2}}{\sigma^{2}}-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})$
$\displaystyle=-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})-\frac{1}{\sigma^{2}}(\tilde{\bm{y}}-\tilde{F}_{1}\bm{g})^{T}(\tilde{\bm{y}}-\tilde{F}_{1}\bm{g})$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}L^{T}[D(\bm{\tau}^{2})]^{-1}L\bm{g}-\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{\bm{y}}+\frac{1}{\sigma^{2}}\bm{g}^{T}\tilde{F}_{1}^{T}\tilde{\bm{y}}$
$\displaystyle\quad+\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{F}_{1}\bm{g}-\frac{1}{\sigma^{2}}\bm{g}^{T}\tilde{F}_{1}^{T}\tilde{F}_{1}\bm{g}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\frac{1}{\sigma^{2}}\bm{g}^{T}\tilde{F}_{1}^{T}\tilde{\bm{y}}+\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{F}_{1}\bm{g}-\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{\bm{y}}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\frac{2}{\sigma^{2}}\bm{g}^{T}\tilde{F}_{1}^{T}\tilde{\bm{y}}-\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{\bm{y}}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\bm{g}^{T}\Gamma\bar{\bm{g}}-\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}+\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}-\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{\bm{y}}$
$\displaystyle=-\frac{1}{2}(\bm{g}-\bar{\bm{g}})^{T}\Gamma(\bm{g}-\bar{\bm{g}})+c(\tilde{\bm{y}}),$

where
$c(\tilde{\bm{y}})=\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}-\frac{1}{\sigma^{2}}\tilde{\bm{y}}^{T}\tilde{\bm{y}}$
is mutually independent of $\bm{g}$. Thus we have mean $\bar{\bm{g}}$ and
precision $\Gamma$, yielding the desired result.

######  Corollary 4.2.

Let $F_{1}=FD(e^{i\phi})$ and suppose $F$ is unitary and $L$ is rank $n$. For
$F_{1}^{H}\bm{y}=\bm{q}\odot e^{i\bm{\varphi}}$ where
$\bm{q}\in\left(\mathbb{R}^{+}\right)^{n}$ and
$\bm{\varphi}\in[-\pi,\pi)^{n}$, the function
$\tilde{f}_{\mathcal{G}|\mathcal{Y},\mathcal{T}^{2},\Phi}(\bm{g}|\bm{y},\bm{\tau}^{2},\bm{\phi})$
in (39) defines a Gaussian density over $\mathbb{R}^{n}$ with center and
precision given by

(57a)
$\bar{\bm{g}}=\Gamma^{-1}\left(\frac{2}{\sigma^{2}}\bm{q}\odot\cos(\bm{\varphi})\right),$
(57b)
$\Gamma=L^{T}[D(\bm{\tau}^{2})]^{-1}L+\frac{2}{\sigma^{2}}\mathbb{I}_{n}.$

######  Proof A.2.

Since $F$ is unitary, we have in (39)

$\displaystyle-\frac{\norm{\bm{y}-F_{1}\bm{g}}^{2}_{2}}{\sigma^{2}}-\frac{1}{2}(L\bm{g})^{T}[D(\bm{\tau}^{2})]^{-1}(L\bm{g})$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\frac{1}{\sigma^{2}}\bm{g}^{T}\left(\bm{q}\odot
e^{i\bm{\varphi}}\right)$
$\displaystyle\quad+\frac{1}{\sigma^{2}}\left(\bm{q}\odot
e^{-i\bm{\varphi}}\right)^{T}\bm{g}-\frac{1}{\sigma^{2}}\bm{y}^{H}\bm{y}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\frac{1}{\sigma^{2}}\sum_{j=1}^{n}g_{j}q_{j}\left(e^{i\varphi_{j}}+e^{-i\varphi_{j}}\right)-\frac{1}{\sigma^{2}}\bm{y}^{H}\bm{y}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\frac{2}{\sigma^{2}}\sum_{j=1}^{n}g_{j}q_{j}\cos(\varphi_{j})-\frac{1}{\sigma^{2}}\bm{y}^{H}\bm{y}$
$\displaystyle=-\frac{1}{2}\bm{g}^{T}\Gamma\bm{g}+\bm{g}^{T}\Gamma\bar{\bm{g}}-\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}+\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}-\frac{1}{\sigma^{2}}\bm{y}^{H}\bm{y}$
$\displaystyle=-\frac{1}{2}(\bm{g}-\bar{\bm{g}})^{T}\Gamma(\bm{g}-\bar{\bm{g}})+c(\bm{y})$

where
$c(\bm{y})=\frac{1}{2}\bar{\bm{g}}^{T}\Gamma\bar{\bm{g}}-\frac{1}{\sigma^{2}}\bm{y}^{H}\bm{y}$
is mutually independent of $\bm{g}$. Thus we have mean $\bar{\bm{g}}$ and
precision $\Gamma$, yielding the desired result.

######  Theorem 4.3.

Let $F_{2}=FD(\bm{g})$ as in (43) and (44) and define $A=F_{2}^{H}F_{2}$ with
elements $[A]_{j,k}=a_{j,k}e^{i\alpha_{j,k}}$, where $a_{j,k}\in\mathbb{R}$
and $\alpha_{j,k}\in[-\pi,\pi)$. Further denote $F_{2}^{H}\bm{y}=\bm{q}\odot
e^{i\bm{\varphi}}$, where $\bm{q}\in\mathbb{R}_{+}^{n}$ and
$\bm{\varphi}\in[-\pi,\pi)^{n}$, with

$\displaystyle
u_{i}=\frac{2q_{i}}{\sigma^{2}}\cos\varphi_{i}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{n}\frac{2a_{i,k}}{\sigma^{2}}\cos(\tilde{\alpha}_{i,k}+\phi_{k}),\quad\quad
v_{i}=\frac{2q_{i}}{\sigma^{2}}\sin\varphi_{i}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{n}\frac{2a_{i,k}}{\sigma^{2}}\sin(\tilde{\alpha}_{i,k}+\phi_{k}),$

where $\tilde{\alpha}_{i,k}=\operatorname{sgn}(k-i)\alpha_{i,k}$ for
$k=1,\dots,n$. Then

$\displaystyle
f_{\Phi_{i}|\Phi_{-i},\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\phi_{i}|\phi_{-i},\bm{y},\bm{g},\bm{\tau}^{2})\propto\pi_{vM}\left(\phi_{i}\Big{|}\mu_{i},\kappa_{i}\right),\quad
i=1,\dots,n,$

where

$\displaystyle\kappa_{i}=\sqrt{u_{i}^{2}+v_{i}^{2}},\quad\quad\mu_{i}=\begin{cases}\arctan\left(-\frac{v_{i}}{u_{i}}\right)&\text{if
}u_{i}>0\\\ \pi/2&\text{if }u_{i}=0\\\
\arctan\left(-\frac{v_{i}}{u_{i}}\right)+\pi&\text{if }u_{i}<0\end{cases}$

and $\pi_{vM}(x|\mu,\kappa)$ is the von Mises probability density function,
given as

$\displaystyle\pi_{vM}(x|\mu,\kappa)=\frac{\exp(\kappa\cos(x-\mu))}{2\pi
I_{0}(\kappa)},$

with location $\mu$ and concentration $\kappa$, and where $I_{0}$ is the
modified Bessel function of the first kind of order zero.

######  Proof A.3.

By direct calculation using (43), we have

$\displaystyle f_{\Phi_{i}|\Phi_{-i},\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}$
$\displaystyle(\phi_{i}|\phi_{-i},\bm{y},\bm{g},\bm{\tau}^{2})\propto\exp(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F_{2}e^{i\bm{\phi}}}^{2}_{2})$
$\displaystyle\propto\exp(-\frac{1}{\sigma^{2}}\left(e^{-i\bm{\phi}^{T}}Ae^{i\bm{\phi}}-\bm{y}^{H}F_{2}e^{i\bm{\phi}}-e^{-i\bm{\phi}^{T}}F_{2}^{H}\bm{y}\right))$
$\displaystyle=\exp\Bigg(-\frac{1}{\sigma^{2}}\bigg{(}\sum\limits_{\begin{subarray}{c}j,k=1\\\
k>j\end{subarray}}^{n}\left(a_{j,k}e^{i(\alpha_{j,k}-\phi_{j}+\phi_{k})}+a_{j,k}e^{-i(\alpha_{j,k}-\phi_{j}+\phi_{k})}\right)-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})\bigg{)}\Bigg{missing})$
$\displaystyle=\exp\Bigg(-\frac{1}{\sigma^{2}}\bigg{(}\sum\limits_{\begin{subarray}{c}j,k=1\\\
k>j\end{subarray}}^{n}2a_{j,k}\cos(\alpha_{j,k}-\phi_{j}+\phi_{k})-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})\bigg{)}\Bigg{missing})$
$\displaystyle\propto\exp\Bigg(-\frac{1}{\sigma^{2}}\bigg{(}\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
i\end{subarray}}^{n}2a_{i,k}\left(\cos\phi_{i}\cos(\phi_{k}+\tilde{\alpha}_{i,k})+\sin\phi_{i}\sin(\phi_{k}+\tilde{\alpha}_{i,k})\right)-2q_{i}\cos(\phi_{i}-\varphi_{i})\bigg{)}\Bigg{missing})$
$\displaystyle=$ $\displaystyle=\exp(u_{i}\cos\phi_{i}+v_{i}\sin\phi_{i})$
$\displaystyle=\exp(\kappa_{i}\cos(\phi_{i}-\mu_{i})).$

######  Theorem 4.4.

Suppose $F_{2}=FD(\bm{g})$, where $F$ is unitary. Let
$F_{2}^{H}\bm{y}=\bm{q}\odot e^{i\bm{\varphi}}$ where
$\bm{q}\in\mathbb{R}_{+}^{n}$ and $\bm{\varphi}\in[-\pi,\pi)^{n}$. Then

$\displaystyle
f_{\Phi|\mathcal{Y},\mathcal{G}}(\bm{\phi}|\bm{y},\bm{g})=\prod_{i=1}^{n}\pi_{vM}\left(\phi_{i}\Big{|}\varphi_{i},\frac{2}{\sigma^{2}}q_{i}\right).$

######  Proof A.4.

Consider the following:

$\displaystyle f_{\Phi|\mathcal{Y},\mathcal{G}}(\bm{\phi}|\bm{y},\bm{g})$
$\displaystyle=C\exp(-\frac{1}{\sigma^{2}}\norm{\bm{y}-F_{2}e^{i\bm{\phi}}}^{2}_{2})\bm{1}_{[-\pi,\pi)}(\bm{\phi})$
$\displaystyle=C\exp(-\frac{1}{\sigma^{2}}\left(e^{-i\bm{\phi}}F_{2}^{H}F_{2}e^{i\bm{\phi}}-\bm{y}^{H}F_{2}e^{i\bm{\phi}}-e^{-i\bm{\phi}}F_{2}^{H}\bm{y}+\bm{y}^{H}\bm{y}\right))\bm{1}_{[-\pi,\pi)}(\bm{\phi})$
$\displaystyle=C\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})+\bm{y}^{H}\bm{y}\right))\bm{1}_{[-\pi,\pi)}(\bm{\phi}),$

where $C$ is a normalization constant. To find $C$, we integrate the above
expression over $\mathbb{R}^{n}$ as follows:

$\displaystyle\frac{1}{C}$
$\displaystyle=\int_{\mathbb{R}^{n}}\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})+\bm{y}^{H}\bm{y}\right))\bm{1}_{[-\pi,\pi)}(\bm{\phi})\mathrm{d}\bm{\phi}$
$\displaystyle=\int_{[-\pi,\pi)^{n}}\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})+\bm{y}^{H}\bm{y}\right))\mathrm{d}\bm{\phi}$
$\displaystyle=\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}+\bm{y}^{H}\bm{y}\right))\int_{[-\pi,\pi)^{n}}\exp(\frac{2}{\sigma^{2}}\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi}))\mathrm{d}\bm{\phi}$
$\displaystyle=\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}+\bm{y}^{H}\bm{y}\right))\prod_{i=1}^{n}2\pi
I_{0}\left(\frac{2}{\sigma^{2}}q_{i}\right).$

Thus, we have

$\displaystyle
f_{\Phi|\mathcal{Y},\mathcal{G},\mathcal{T}^{2}}(\bm{\phi}|\bm{y},\bm{g},\bm{\tau}^{2})$
$\displaystyle=\frac{\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}-2\bm{q}^{T}\cos(\bm{\phi}-\bm{\varphi})+\bm{y}^{H}\bm{y}\right))\bm{1}_{[-\pi,\pi)}(\bm{\phi})}{\exp(-\frac{1}{\sigma^{2}}\left(\bm{g}^{T}\bm{g}+\bm{y}^{H}\bm{y}\right))\prod_{i=1}^{n}2\pi
I_{0}\left(\frac{2}{\sigma^{2}}q_{i}\right)}$
$\displaystyle=\prod_{j=1}^{n}\frac{\exp(\frac{2}{\sigma^{2}}q_{j}\cos(\phi_{j}-\varphi_{j}))}{2\pi
I_{0}\left(\frac{2}{\sigma^{2}}q_{j}\right)}\bm{1}_{[-\pi,\pi)}(\phi_{j})$
$\displaystyle=\prod_{j=1}^{n}\pi_{vM}\left(\phi_{j}\Big{|}\varphi_{j},\frac{2}{\sigma^{2}}q_{j}\right),$

which completes the proof.

## References

  * [1] E. Aguilera, M. Nannini, and A. Reigber, Wavelet-based compressed sensing for SAR tomography of forested areas, IEEE transactions on geoscience and remote sensing, 51 (2013), pp. 5283–5295.
  * [2] R. Archibald, A. Gelb, and R. B. Platte, Image reconstruction from undersampled Fourier data using the polynomial annihilation transform, Journal of Scientific Computing, 67 (2016), pp. 432–452.
  * [3] J. Barber, A generalized likelihood ratio test for coherent change detection in polarimetric SAR, IEEE Geoscience and Remote Sensing Letters, 12 (2015), pp. 1873–1877.
  * [4] D. Best and N. I. Fisher, Efficient simulation of the von Mises distribution, Journal of the Royal Statistical Society: Series C (Applied Statistics), 28 (1979), pp. 152–157.
  * [5] Z. I. Botev, J. F. Grotowski, and D. P. Kroese, Kernel density estimation via diffusion, (2010).
  * [6] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine learning, 3 (2011), pp. 1–122.
  * [7] E. J. Candes, M. B. Wakin, and S. P. Boyd, Enhancing sparsity by reweighted $\ell_{1}$ minimization, Journal of Fourier Analysis and Applications, 14 (2008), pp. 877–905.
  * [8] M. Cetin, W. C. Karl, and A. S. Willsky, Feature-preserving regularization method for complex-valued inverse problems with application to coherent imaging, Optical Engineering, 45 (2006), pp. 017003–017003.
  * [9] M. Cheney, A mathematical tutorial on synthetic aperture radar, SIAM review, 43 (2001), pp. 301–312.
  * [10] V. Churchill and A. Gelb, Sampling-based spotlight SAR image reconstruction from phase history data for speckle reduction and uncertainty quantification, SIAM/ASA Journal on Uncertainty Quantification, 10 (2022), pp. 1225–1249.
  * [11] V. Churchill and A. Gelb, Sub-aperture SAR imaging with uncertainty quantification, Inverse Problems, 39 (2023), p. 054004.
  * [12] D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory, 52 (2006), pp. 1289–1306.
  * [13] H. Duan, L. Zhang, J. Fang, L. Huang, and H. Li, Pattern-coupled sparse bayesian learning for inverse synthetic aperture radar imaging, IEEE Signal Processing Letters, 22 (2015), pp. 1995–1999.
  * [14] R. Gens and J. L. Van Genderen, Review article SAR interferometry—issues, techniques, applications, International journal of remote sensing, 17 (1996), pp. 1803–1835.
  * [15] J. Glaubitz, A. Gelb, and G. Song, Generalized sparse bayesian learning and application to image reconstruction, SIAM/ASA Journal on Uncertainty Quantification, 11 (2023), pp. 262–284.
  * [16] D. Green, J. Jamora, and A. Gelb, Leveraging joint sparsity in 3D synthetic aperture radar imaging, Applied Mathematics for Modern Challenges, 1 (2023), pp. 61–86.
  * [17] M. R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, Journal of research of the National Bureau of Standards, 49 (1952), pp. 409–436.
  * [18] J. Immerkaer, Fast noise variance estimation, Computer vision and image understanding, 64 (1996), pp. 300–302.
  * [19] C. V. Jakowatz, D. E. Wahl, P. H. Eichel, D. C. Ghiglia, and P. A. Thompson, Spotlight-mode synthetic aperture radar: a signal processing approach, Springer Science & Business Media, 2012.
  * [20] M. Järvenpää and R. Piché, Bayesian hierarchical model of total variation regularisation for image deblurring, arXiv preprint arXiv:1412.4384, (2014).
  * [21] J. P. Kaipio, V. Kolehmainen, E. Somersalo, and M. Vauhkonen, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse Problems, 16 (2000), p. 1487.
  * [22] I. Loris, G. Nolet, I. Daubechies, and F. Dahlen, Tomographic inversion using $\ell$1-norm regularization of wavelet coefficients, Geophysical Journal International, 170 (2007), pp. 359–370.
  * [23] M. Lustig, D. Donoho, and J. M. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 58 (2007), pp. 1182–1195.
  * [24] H. Maatouk and X. Bay, A new rejection sampling method for truncated multivariate Gaussian random variables restricted to convex sets, in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014, Springer, 2016, pp. 521–530.
  * [25] J. R. Michael, W. R. Schucany, and R. W. Haas, Generating random variates using transformations with multiple roots, The American Statistician, 30 (1976), pp. 88–90.
  * [26] L. M. Novak, Coherent change detection for multi-polarization SAR, in Conference Record of the Thirty-Ninth Asilomar Conference on Signals, Systems and Computers, 2005., IEEE, 2005, pp. 568–573.
  * [27] A. Pakman and L. Paninski, Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians, Journal of Computational and Graphical Statistics, 23 (2014), pp. 518–542.
  * [28] G. Papandreou and A. L. Yuille, Gaussian sampling by local perturbations, Advances in Neural Information Processing Systems, 23 (2010).
  * [29] T. Park and G. Casella, The Bayesian Lasso, Journal of the American Statistical Association, 103 (2008), pp. 681–686.
  * [30] A. Rajagopal, J. Hilton, D. Boutte, A. P. Brown, and J. R. Jamora, Enhanced compressed sensing 3D SAR imaging via cross-modality EO-SAR joint-sparsity priors, in Algorithms for Synthetic Aperture Radar Imagery XXX, vol. 12520, SPIE, 2023, pp. 13–23.
  * [31] T. Sanders, A. Gelb, and R. B. Platte, Composite SAR imaging using sequential joint sparsity, Journal of Computational Physics, 338 (2017), pp. 357–370.
  * [32] T. Sanders, R. B. Platte, and R. D. Skeel, Effective new methods for automated parameter selection in regularized inverse problems, Applied Numerical Mathematics, 152 (2020), pp. 29–48.
  * [33] T. Scarnati and A. Gelb, Joint image formation and two-dimensional autofocusing for synthetic aperture radar data, Journal of Computational Physics, 374 (2018), pp. 803–821.
  * [34] L. A. Shepp and B. F. Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science, 21 (1974), pp. 21–43.
  * [35] R. Tibshirani, Regression shrinkage and selection via the Lasso, Journal of the Royal Statistical Society Series B: Statistical Methodology, 58 (1996), pp. 267–288.
  * [36] R. Touzi, A. Lopes, J. Bruniquel, and P. W. Vachon, Coherence estimation for SAR imagery, IEEE transactions on geoscience and remote sensing, 37 (1999), pp. 135–149.
  * [37] M. Vono, N. Dobigeon, and P. Chainais, High-dimensional Gaussian sampling: a review and a unifying approach based on a stochastic proximal point algorithm, SIAM Review, 64 (2022), pp. 3–56.
  * [38] D. E. Wahl, D. A. Yocky, C. V. Jakowatz, and K. M. Simonson, A new maximum-likelihood change estimator for two-pass SAR coherent change detection, IEEE Transactions on Geoscience and Remote Sensing, 54 (2016), pp. 2460–2469.
  * [39] L. P. Yaroslavskii and N. S. Merzlyakov, Methods of digital holography, Springer, 1980.
  * [40] J. T. Ylitalo and H. Ermert, Ultrasound synthetic aperture imaging: monostatic approach, IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 41 (1994), pp. 333–339.
  * [41] S. Yu, A. S. Khwaja, and J. Ma, Compressed sensing of complex-valued data, Signal Processing, 92 (2012), pp. 357–362.
  * [42] J. Zhang, A. Gelb, and T. Scarnati, Empirical Bayesian inference using a support informed prior, SIAM/ASA Journal on Uncertainty Quantification, 10 (2022), pp. 745–774.