remarkRemark hypothesisHypothesis claimClaim exampleExample

# Empirical Bayesian Inference using Joint Sparsity

Jiahui Zhang Department of Mathematics, University of Minnesota, Minneapolis,
MN<EMAIL_ADDRESS>Anne Gelb Department of Mathematics, Dartmouth
College, Hanover, NH<EMAIL_ADDRESS>Theresa Scarnati Qualis
Corporation, Huntsville, AL<EMAIL_ADDRESS>
###### Abstract

This paper develops a new empirical Bayesian inference algorithm for solving a
linear inverse problem given multiple measurement vectors (MMV) of under-
sampled and noisy observable data. Specifically, by exploiting the joint
sparsity across the multiple measurements in the sparse domain of the
underlying signal or image, we construct a new support informed sparsity
promoting prior. Several applications can be modeled using this framework, and
as a prototypical example we consider reconstructing an image from synthetic
aperture radar (SAR) observations using nearby azimuth angles. Our numerical
experiments demonstrate that using this new prior not only improves accuracy
of the recovery, but also reduces the uncertainty in the posterior when
compared to standard sparsity producing priors.

###### keywords:

Empirical Bayesian inference, uncertainty quantification, sparsity promoting
priors, multiple measurement vectors, joint sparsity

62C12,65C40,68U10

##  1 Introduction

We consider the linear inverse problem modeled as
$\mathbf{Y}=A\mathbf{X}+\mathbf{E}$ where $\mathbf{X},\mathbf{Y},\mathbf{E}$
are random variables defined over the common probability space
$(\Omega,\mathcal{F},\mathbb{P})$, and $\mathbf{X}$ and $\mathbf{E}$ are
assumed to be independent. Here $\mathbf{X}$ represents the underlying unknown
we seek to recover, $\mathbf{Y}$ is the observable data,
$A\in\mathbb{C}^{n\times n}$ is a deterministic and known forward operator,
and $\mathbf{E}\sim\mathcal{CN}(0,\sigma^{2}\mathbb{I}_{{n}})$ models the
additive complex Gaussian noise. The realizations of the observable data
$y_{j}=Ax+e_{j}$, $j=1,\dots,J$, are obtained as a collection of $J>1$
measurements, and are often referred to as multiple measurement vectors
(MMVs), [13]. Many real world phenomena may be modeled in this way, and our
study is motivated by applications in synthetic aperture radar (SAR) imaging.
In this case, the observable data are multiple sets of phase history data, $A$
is a (non-uniform) Fourier transform matrix, and we seek to recover a complex
SAR image ${\mathbf{X}}$. To simplify our presentation and analysis, our
investigation is limited to one-dimensional real-valued signals. It is
straightforward to extend our approach to multi-dimensions as well as complex-
valued signals and images, however. Additionally, our approach is conducive to
parallelizable algorithms.

A common assumption made in signal111We use the terms signal, image, and
function interchangeably throughout our manuscript. recovery, which we also
use in this investigation, is that the underlying signal is sparse in some
domain (e.g. the gradient or edge domain). In such cases compressive sensing
(CS) algorithms, [9, 18], are frequently employed, as they are specifically
designed to exploit these sparsity assumptions. Constructed as a convex
optimization problem, a CS algorithm seeks to minimize the $\ell_{2}$ norm of
the data fidelity regularized by the $\ell_{1}$ norm of the solution in its
sparse domain. For example, given a particular observation $y_{j}=Ax+e_{j}$,
the solution is obtained as

(1)
$x^{\ast}={\underset{x\in\mathbb{R}^{n}}{\operatorname{argmin}}}\Big{\\{}\lambda||\mathcal{L}x||_{1}+\frac{1}{2}||y_{j}-Ax||_{2}^{2}\Big{\\}},$

where $\mathcal{L}$ is a pre-determined sparsifying transform and the
regularization parameter $\lambda>0$ is chosen to balance each term’s
contribution. For presumably less error in the fidelity term,
$||y_{j}-Ax||_{2}^{2}$, a smaller $\lambda$ is chosen and vice versa. In this
context the $\ell_{1}$ norm is typically viewed as a surrogate for the
$\ell_{0}$ pseudo-norm since using the latter yields an intractable problem,
[11]. The additional information from the multiple observations collected in
the MMVs can be further exploited by employing the joint sparsity assumption.
In a nutshell, joint sparsity leverages the assumption that the sparse domain
of the underlying signal should be similar across all collected measurements
for that same signal.

Many techniques have been introduced to improve the robustness and efficacy of
CS algorithms in both the single and multiple measurement cases, [10, 16, 1,
21, 33]. However, since the recovery is limited to point estimate solutions,
it is not possible to quantify the uncertainty. This disadvantage is non-
trivial because in practice it is often crucial to know how reliable the
recovery is, especially in the case of noisy or limited data availability
within in each measurement vector. Moreover, even within the same application,
CS algorithms inevitably require hand tuning of the regularization parameter,
with small changes often resulting in drastically different recovery outcomes.

Due to these inherent limitations, the Bayesian approach, [25], is often used
to recover the posterior distribution of the unknown random variable ${\bf X}$
with respect to ${\mathbf{Y}}$ as related by the model
${\mathbf{Y}}=A{\mathbf{X}}+{\mathbf{E}}$. We note that from the Bayesian
perspective, the CS methodology may be viewed as obtaining a maximum a
posteriori (MAP) estimate [27], which represents a mode of the posterior
distribution. However, the MAP estimate of the posterior distribution may not
always provide the most descriptive approximation of the true underlying
signal and furthermore does not yield the required information about the
posterior distribution to quantify the uncertainty of the recovery. Therefore
this investigation aims to recover the posterior distribution representing the
unknown $\mathbf{X}$ using sampling methods.

In the Bayes formulation, the posterior density of the unknown $\mathbf{X}$ is
proportional to the product of the likelihood and prior densities, with the
likelihood corresponding to the fidelity term in a convex optimization problem
and the prior to the regularization term. Determining an appropriate prior to
quantify assumptions about the unknown signal is well known to be difficult.
To this end, the sparsity assumption suggests that the $\ell_{1}$ (Laplace)
prior may be suitable for this purpose, and indeed the corresponding MAP
estimate results in the same convex CS-based optimization problem in (1).
Although effective in many applications, the standard Laplace prior does not
really include all of the prior beliefs about the underlying image.
Specifically, no consideration is given to the locations of the support within
the sparse domain. We note that this is also the case in standard CS
algorithms, as the regularization parameter is constant throughout the
sparsity domain.

This issue has been addressed in CS algorithms using various forms of weighted
$\ell_{1}$ regularization, [16, 10, 1, 21]. In all cases, the spatially
varying weighting parameters are designed to promote solutions that are sparse
only in regions that are likely to be sparse, as opposed to promoting global
sparsity. Unfortunately, many of the proposed re-weighted $\ell_{1}$ (and
$\ell_{2}$) regularization algorithms are iterative in nature and sensitive to
noise. A review of the advantages and shortcomings of some of these algorithms
may be found in [14].

In [1, 21], where the MMV case was specifically considered, the spatially
varying weighting parameter was designed by exploiting the joint sparsity
across the measurements in the sparsity domain. In particular, the regions of
support in the sparse domain were determined to be those in which there were
large variances in the measurements. By contrast, small variances in the
measurements indicated sparse regions. The spatially varying weights for the
weighted $\ell_{1}$ (or $\ell_{2}$) regularization were then constructed
accordingly. The technique, coined the variance based joint sparsity (VBJS)
method, proved to be efficient and robust in recovering images from noisy and
under-sampled data. The method can also be used for a single set of
measurements by processing it multiple times to obtain MMVs, [33].

Inspired by these results, in the empirical Bayesian approach used in this
investigation, we introduce a new prior that is constructed by exploiting the
joint sparsity assumption. As is typically done, our prior assumes that the
signal is sparse in the sparse domain. But in our new support informed sparse
prior, we also incorporate the likely regions of support in the sparse domain,
thereby more accurately describing our prior beliefs about the underlying
signal. As in the VBJS method, our technique uses the variance of the
measurements in the sparsity domain to determine regions of support. We then
use this information to construct a binary edge mask for the prior. That is,
the weighting parameter in the prior is binary with a zero value assigned
whenever an edge is detected and one otherwise.222Binary edge masks were also
used in the edge-adaptive $\ell_{2}$ regularization image reconstruction
method developed in [14]. In this way we obtain a truly hierarchical prior,
which will be described in more detail throughout Section 2. Since our new
method uses observable data to inform the prior it may be categorized as an
empirical Bayesian inference technique.

For context, we can compare our method to another empirical Bayesian inference
technique, sparse Bayesian learning (SBL), [34, 35]. In SBL, one seeks to
learn the distributions of the hyperparameters of the prior. For a conjugate
prior (given in [34, 35] by a Gaussian prior with Gamma distributed
hyperparameters), these hyperparameters can be approximated using expectation
maximization (EM) [17]. SBL with EM has the same effect as defining a sparse
prior that is additionally informed by the supports in the sparse domain. In
this regard it is important to point out that (i) SBL is designed specifically
to recover sparse signals. That is, the sparse prior is in the domain of
recovery. By contrast our new method is more general, since it also considers
priors that are sparse in a transformed space; (ii) EM is not computationally
efficient for higher dimensional data. While other more efficient methods may
be used to learn the hyperparameters, the resulting accuracy is considerably
reduced, [35]; and (iii) due to the use of conjugate priors, the density
yields a closed form representation of the posterior distribution, and thus
sampling can be done using the Gibbs sampling method, an algorithm in the
family of Markov Chain Monte-Carlo (MCMC) techniques. Our proposed technique,
on the other hand, is not restricted to using only conjugate priors.
Consequently, the posterior distribtion does not have a closed form solution,
which means that we cannot use the Gibbs sampling method. Hence instead we
employ the celebrated Metropolis Hastings (MH) MCMC algorithm, [20, 30]. For
one dimensional problems, the efficiency of the MH algorithm is comparable to
Gibbs sampling. In higher dimensions, Gibbs sampling is inevitably a faster
sampling technique. The evaluation of the posteriors derived in this paper is
easily parallelized across each dimension, however, so that for many
applications enforcing conjugate priors to enable Gibbs sampling does not seem
necessary.

Our new method may also be compared those of [5, 7, 8], which in sequel we
will refer to as the hierarchical Gaussian prior (HGP) approach. In contrast
to SBL, HGP techniques are not limited to the recovery of sparse signals, as
long as the sparsifying transform in the prior is symmetric positive definite.
This will ensure a closed form solution for the posterior which in turn allows
Gibbs sampling to be used.333Based on its derivation, it seems that EM could
also be used in the case of SPD sparsifying transforms. This would make SBL
applicable to a broader range of signals (not only sparse signals). While such
applications should be considered in future investigations, it does not negate
SBL’s issue with implementation cost in multi-dimensions. As in SBL, the prior
in HGP is typically Gaussian with Gamma distributed hyperparameters. Unlike
SBL, these hyperparameters are not approximated by EM. In order to maintain a
closed form solution for the posterior, the hyperparameters are scalar, that
is, the prior is not informed by the signal’s support in the sparsity domain.
In this regard, a point estimate of the HGP is more similar to the original CS
approach which does not use spatially varying weighted $\ell_{1}$
regularization. Of course HGP has the benefit of being able to quantify the
uncertainty of the solution. We note again that in our new approach we use a
hierarchical support informed prior, where the support in the sparse domain is
determined by exploiting the joint sparsity of the multiple measurements in
the sparse domain.

The rest of this paper is organized as follows. In Section 2 we review the
Bayesian approach for solving inverse problems and introduce our new support
informed sparse prior. We discuss how the resulting posterior is evaluated via
the MH MCMC algorithm in Section 3. Numerical experiments used to evaluate our
method are performed in Section 4. Finally, Section 5 provides some concluding
remarks.

##  2 The Bayesian approach for solving inverse problems

We consider the linear forward model of the form

(2) $\mathbf{Y}=A\mathbf{X}+\mathbf{E},$

where $\mathbf{X},\mathbf{Y},\mathbf{E}$ are random variables defined over a
common probability space $(\Omega,\mathcal{F},\mathbb{P})$, and $\mathbf{X}$
and $\mathbf{E}$ are assumed to be mutually independent. Here $\mathbf{X}$
represents the underlying unknown we seek to recover, $\mathbf{Y}$ is the
observable data, $A\in\mathbb{C}^{n\times n}$ is a known forward operator, and
$\mathbf{E}\sim\mathcal{CN}(0,\sigma^{2}\mathbb{I}_{{n}})$ models additive
complex Gaussian noise. A particular observation $y_{j}$ of (2), $j=1,\dots
J$, where $J$ is the given number of MMVs, is given by

(3) $y=Ax+e,$

where we have dropped the subscript $j$ to avoid cumbersome notation. As noted
in the introduction, sparsity based (CS) reconstruction algorithms are often
used to obtain the optimal solution to (3) as a point estimate solution of
(2). Also as discussed in the introduction, exploring the entire posterior
probability of $\bm{X}$ is valuable in assessing the uncertainty of the our
recovery. This can be achieved within a Bayesian framework, which we briefly
review below. A complete overview of Bayesian inverse problems may be found in
[25].

###  2.1 Bayes’ Theorem

Let $x,y,e\in{\mathbb{R}^{n}}$ and $\lambda\in[0,1]$ be realizations of the
random variables $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{E}$ and $\bm{\lambda}$
respectively. Bayes’ theorem is given by

(4)
$f_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)=\frac{f_{\mathbf{Y|X}}(y|x){f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda)}{f_{\mathbf{Y}}(y)},$

where ${f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)$ is the posterior
density function we seek to recover, $f_{\mathbf{Y}|\mathbf{X}}(y|x)$ is the
likelihood function, and
$\tilde{f}_{\mathbf{X},\bm{\lambda}}(x,\lambda)={f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda)$
is the prior probability density function. The prior is broken into two parts,
the conditional prior ${f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)$ and the
hyper-prior $f_{\bm{\lambda}}(\lambda)$. Further, $f_{\mathbf{Y}}(y)$, often
called the evidence, is given by the total law of probability as

(5)
$f_{\mathbf{Y}}(y)=\int_{\mathbb{R}^{n}}\int_{[0,1]}f_{\mathbf{X},\mathbf{Y},\bm{\lambda}}(x,y,\lambda)d\lambda
dx=\int_{\mathbb{R}^{n}}\int_{[0,1]}f_{\mathbf{Y|X}}(y|x){f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda)d\lambda
dx.$

It is reasonable to assume that (5) is nonzero because otherwise the
observation $y\in\mathbb{R}^{n}$ has a probability of zero and is therefore
irrelevant in the applications we are considering. Moreover, since calculating
(5) is often computationally intractable, particularly in high dimensional
domains, a closed form posterior density is often unavailable. Thus sampling
methods such as those described in Section 3.2 are typically used to
approximate the un-normalized version of the posterior (4) given by

(6)
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)=f_{\mathbf{Y|X}}(y|x){f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda).$

In our numerical experiments we use the Metropolis-Hastings (MH) algorithm to
draw samples of (6).

In what follows, Section 2.2 defines the likelihood $f_{\mathbf{Y|X}}(y|x)$
with respect to the proposed linear model (2), while Section 2.3 describes how
we estimate the hyper-prior $\lambda$. Section 2.4 focuses on the definition
of the conditional prior ${f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)$. There we
will first discuss how sparsity enforcing priors are often employed in
practice and then introduce our novel support informed sparse prior that is
empirically constructed by exploiting the joint sparsity across the multiple
measurements in the sparse domain of the underlying signal. Finally, in
Section 2.5 we define and explore the posterior
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)$.

###  2.2 Estimating the likelihood

Returning to the model in (2), observe that

$f_{\mathbf{Y}|\mathbf{X},\mathbf{E}}(y|x,e)=\delta(y-Ax-e),$

where $\delta(\cdot)$ is the Dirac delta function. Marginalizing with respect
to the noise $e$ yields

$\displaystyle f_{\mathbf{Y}|\mathbf{X}}(y|x)$
$\displaystyle=\int_{\mathbb{R}^{n}}f_{\mathbf{Y}|\mathbf{X},\mathbf{E}}(y|x,e)f_{\mathbf{E}|\mathbf{X}}(e|x)de$
$\displaystyle=\int_{\mathbb{R}^{n}}\delta(y-Ax-e)f_{\mathbf{E}|\mathbf{X}}(e|x)de$
$\displaystyle=f_{\mathbf{E}|\mathbf{X}}(y-Ax|x)$ (7)
$\displaystyle=f_{\mathbf{E}}(y-Ax),$

with the last equality resulting from the independence between $\mathbf{X}$
and $\mathbf{E}$, [25]. Hence we see that the likelihood function,
$f_{\mathbf{Y}|\mathbf{X}}(y|x)$, depends on the noise variable $\mathbf{E}$.
Since $\mathbf{E}$ is Gaussian with mean $0$ and covariance matrix
$\Sigma=\sigma^{2}\mathbb{I}_{{n}}$, the noise density is given by

(8)
$f_{\mathbf{E}}(e)\propto\exp\\{-\frac{1}{2\sigma^{2}}(e-0)^{T}(e-0)\\}=\exp\\{-\frac{||e||^{2}_{2}}{2\sigma^{2}}\\}.$

Thus from (7) and (8) we obtain the likelihood density

(9)
$f_{\mathbf{Y}|\mathbf{X}}(y|x)=f_{\mathbf{E}}(y-Ax)\propto\exp\\{-\frac{1}{2\sigma^{2}}||y-Ax||^{2}_{2}\\}.$

###  2.3 Estimating the hyper-prior

We use $K$-fold cross validation to estimate the hyper-prior, [27]. To
summarize this procedure, we begin by approximating

(10) $\displaystyle f_{\bm{\lambda}}(\lambda)\approx
f_{\bm{\lambda}|\mathbf{Y}}(\lambda|y)=\frac{f_{\mathbf{Y}|\bm{\lambda}}(y|\lambda)\tilde{f}_{\bm{\lambda}}(\lambda)}{f(y)}\propto
f_{\mathbf{Y}|\bm{\lambda}}(y|\lambda)\tilde{f}_{\bm{\lambda}}(\lambda),$

where $\tilde{f}_{\bm{\lambda}}(\lambda)$ is the prior on the hyper-prior. By
assuming $\lambda\sim U[0,1]$ so that $\tilde{f}_{\bm{\lambda}}(\lambda)=1$,
we have

(11) $\displaystyle f_{\bm{\lambda}}(\lambda)\propto
f_{\mathbf{Y}|\bm{\lambda}}(y|\lambda).$

The mode of the distribution (11) can be approximated using the MAP estimate
as

(12)
$\hat{\lambda}={\underset{\lambda}{\operatorname{argmax}}}\\{f_{\mathbf{Y}|\bm{\lambda}}(y|\lambda)\\}={\underset{\lambda}{\operatorname{argmax}}}\int_{\mathbb{R}^{n}}f_{\mathbf{Y}|\mathbf{X}}(y|x)f_{\mathbf{X}|\bm{\lambda}}(x|\lambda)dx.$

Integration with respect to $x$ over the space $\mathbb{R}^{n}$ is not
computationally feasible. We therefore instead employ $K$-fold cross-
validation to approximate $\hat{\lambda}$ given the data $y_{j}\in\Omega$ for
$j=1,\dots,J$. The general process for $K$-fold cross-validation is to first
calculate $M$ MAP estimates of the unknown from a subset of $M<J$ observable
measurements using sample candidates of the hyper-prior $\hat{\lambda}$. These
are then used to generate $M$ training vectors, which are in turn compared to
the remaining $J-M$ observable measurements, or testing vectors. The whole
process is then repeated for $K$ independent trials, and $\hat{\lambda}$ is
chosen to minimize the mean square error of (2) between all training and
testing vectors. The procedure is summarized in Algorithm 1.

Algorithm 1 $K$-Fold Cross-Validation to Estimate $\hat{\lambda}$.

1: Input: Set of multiple measurements (MMVs) ${\bf
y}=[y_{1},y_{2},y_{3},\ldots,y_{J}]$.

2: Output: Point Estimate of the hyper-prior $\hat{\lambda}$.

3: for $k=1$ to $K$ do

4: Randomly partition MMVs into $M$ training vectors and $J-M$ testing
vectors.

5: for $i=1$ to $M$ do

6: Sample candidate hyper-prior from $\tilde{\lambda}_{i,k}\sim U[0,1]$.

7: Calculate the MAP estimate of the unknown as
$\hat{x}_{i,k}={\underset{x}{\operatorname{argmax}}}\\{f_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\tilde{\lambda}_{i,k}|y_{i})\\}.$

8: Compute training data according to $\hat{y}_{i,k}=A\hat{x}_{i,k}.$

9: Evaluate the training data by calculating
$E_{i,k}=\frac{1}{J-M}\sum_{l=M+1}^{J}\text{MSE}(\hat{y}_{i,k},y_{l}),$ where
$\\{y_{l}\\}_{l=M+1}^{J}$ are the partitioned $J-M$ testing vectors.

10: end for

11: end for

12: Choose $(i^{*},k^{*})={\underset{i,k}{\operatorname{argmin}}}\hskip
2.84526ptE_{i,k}$ and set $\hat{\lambda}=\tilde{\lambda}_{i^{*},k^{*}}$.

The hyper-prior $\hat{\lambda}$ recovered from Algorithm 1 approximates the
solution to (12) such that

(13) $f_{\bm{\lambda}}(\lambda)\approx\delta_{\hat{\lambda}}(\lambda).$

Intuitively, (13) implies that we are simply concentrating the density (11) on
the value $\hat{\lambda}$, which is our estimate of the mode of the hyper-
prior distribution (12). Note that the $K$ fold validation process assumes
that the distribution for $\lambda$ is unimodal. Previous efforts, see e.g.
[27], have also demonstrated that (13) is a good approximation of (12).

######  Remark 2.1.

The hyper-prior $\hat{\lambda}$ recovered in Algorithm 1 is consistent with
how one ideally chooses the regularization parameter in the CS approximation
given by (1), or as we will again see in (22), where the fidelity and
regularization parameters are combined to form a single parameter.
Specifically this parameter should be chosen to offset the variance of the
noise, and the $K$-fold cross validation method provides an empirical approach
to calculating this variance.

###  2.4 A support informed sparsity prior

The main focus of this investigation is to construct a support informed
sparsity prior
$\tilde{f}_{\mathbf{X},\bm{\lambda}}(x,\lambda)={f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda)$
in (4). As will be demonstrated in what follows, this will reduce the
uncertainty of the posterior for models given by (2).

Using the approximation in (13) we can make the approximation

$\tilde{f}_{\mathbf{X},\bm{\lambda}}(x,\lambda)={f}_{\mathbf{X}|\bm{\lambda}}(x|\lambda)f_{\bm{\lambda}}(\lambda)\approx{f}_{\mathbf{X}|\bm{\lambda}}(x|\hat{\lambda}),$

and it remains to define the prior through
${f}_{\mathbf{X}|\bm{\lambda}}(x|\hat{\lambda})$. To do so, we make two
assumptions: (i) the prior enforces the sparsity of some transformation of
$\mathbf{X}$ and (ii) there are multiple observations $y$ of $\mathbf{Y}$ for
each realization $x$ of $\mathbf{X}$.444Similarly to what was done in [33, 1],
it is possible to relax this second assumption and instead have a single
measurement vector processed multiple times. For purposes of brevitity we
exclude such considerations in this investigation.

There are many applications for which it is reasonable to assume that the
desired solution is sparse in some domain. In our investigation, the
underlying signal is assumed to be piecewise smooth, implying that the edge
domain is sparse. It is important to note that as a departure from the
traditional sparse priors commonly used in formulating the posterior, our new
support informed sparse prior considers both the existence of the sparse
domain as well as the locations of the non-zero values within it, i.e., its
support.

####  2.4.1 Sparsity promoting priors

Within the Bayesian formulation, sparsity may be enforced through the Laplace
prior given by

(14)
$f_{\mathbf{X}|\bm{\lambda}}(x|\hat{\lambda})\propto\exp\\{-\hat{\lambda}||\mathcal{L}x||_{1}\\},$

where the operator $\mathcal{L}$ is a transform operator used to approximate
the sparse domain. In this investigation we choose $\mathcal{L}$ to be the
polynomial annihilation (PA) transform, [3, 21], which is designed to
approximate the corresponding edge function of a piecewise smooth function
when acting on the unknown (given by $x$ in (14)) at specified grid points in
the spatial domain.555A short description of the PA transform is provided in
Appendix A. We note that our new method is not defined by the choice of
sparsifying transform (here ${\mathcal{L}}$), nor by the sparsity domain (here
the edge domain), and that other choices may be more appropriate for different
applications.

Regardless of how both the hyper-parameter $\hat{\lambda}$ and sparsifying
transform matrix $\mathcal{L}$ are determined, the Laplace prior given in (14)
may not provide the most information about what is known about the underlying
signal. As discussed already, a better prior would make use of the support
locations in the sparse domain. Therefore, following the ideas in [1, 21], we
develop a new support informed sparse prior. In so doing, we better capture
the piecewise smoothness of the underlying signal. A couple of remarks are in
order:

######  Remark 2.2.

Our numerical experiments confirm that using the PA transform with $m>1$ is
appropriate when constructing the support informed sparse prior for a
piecewise smooth function. In particular using $m=1$, which is equivalent to
total variation (TV), does not perform as well. In general when the underlying
function is not constant in smooth domains, a higher degree polynomial
annihilation is needed to recover an accurate representation of the sparse
domain. This issue is extensively discussed for the compressive sensing
framework in [3], where higher order PA is used to avoid the “staircasing”
effect caused by TV regularization. Similarly, the high order PA describes a
better prior in our approach used here.

######  Remark 2.3.

To be clear, there are a variety of sparsity promoting priors to choose from,
including but not limited to the spike and slab prior and the horseshoe prior
[6, 28]. Hierarchical Gaussian priors, which provides the basic framework for
the methods described in [5, 34], may also be used. Regardless of the type of
prior chosen, some transform operator ${\mathcal{L}}$ to the sparse domain is
still needed. This investigation does not intend to provide a complete survey
of the advantages and drawbacks of using various sparsity promoting priors.
Instead we demonstrate that our proposed prior is a departure from the
aforementioned priors because it explicitly includes the support locations in
the sparse domain. Indeed, such knowledge may be incorporated into some of the
aforementioned choices. Even further, our approach may be used to form mixed
priors. That is, because our spatially adaptive prior encourages the
separation of scales in the sparse domain, different priors may be used in
different regions of the domain. The Laplace prior in (14) is mentioned as a
prototype, and is a common choice when interested in constructing the MAP
estimate of the posterior (CS approach). As will be demonstrated by our
numerical examples, our support informed sparse prior yields better accuracy
and reduced uncertainty than other priors typically used. More information
about sparse priors may be found in [28].

####  2.4.2 Joint sparsity

As noted in the introduction, many applications collect multiple measurement
vectors (MMVs) of observable data. By exploiting the joint sparsity in the
sparse domain of the underlying signal that can be obtained from the MMVs, a
technique was developed in [1, 21] to improve the robustness of the MAP
estimate.666We note that the method was derived from the CS perspective. For
reasons that will become apparent later in this section, the technique was
coined the Variance Based Joint Sparsity (VBJS) method. We now adapt VBJS to
construct a new support informed sparse prior. We then demonstrate that this
new prior yields reduced uncertainty and improved efficiency in the sampling
of the posterior in (6) when compared to other typical prior definitions, such
as the Laplace prior in (14).

As previously mentioned, the domain of the problem is $\mathbb{R}^{n}$ so that
the unknown $x$ may be represented by a $n$ length vector. A $d$-dimensional
problem will then yield an $n^{d}$ length vector representation of the
unknown, and as already noted our method is easily parallelizable since it
samples the posterior for each pixel of interest. Suppose that we have
multiple vector realizations of $y_{j}\sim\mathbf{Y}\in\mathbb{R}^{n}$ to form
$Y=[y_{1},y_{2},y_{3},\ldots,y_{J}]$. In practice, the data vectors $y_{j}$
are often noisy indirect measurements of the same underlying unknown signal
represented by a random variable $\mathbf{X}$. It is therefore realistic to
suppose that they should share similar support in a chosen sparse domain, for
instance in the edge domain. This is the joint sparsity assumption. As noted
previously, other sparse domains may also be used for this purpose.

In order to make use of the joint sparsity concept, we first provide some
definitions and notation, which may be found in [1]. Let $p\in\mathbb{R}^{n}$.
As is conventional, we write ${\left\|\cdot\right\|}_{0}$ for the
$\ell^{0}$-‘norm’, i.e.

${\left\|p\right\|}_{0}=\left|\mathrm{supp}(p)\right|,$

where $\mathrm{supp}(p)$ is the support of $p=(p_{i})^{n}_{i=1}$ defined by

$\mathrm{supp}(p)=\\{i:p_{i}\neq 0\\}.$

Given a vector $\bm{w}=(w_{i})^{n}_{i=1}$ of positive weights, we define the
weighted $\ell^{1}_{\bm{w}}$-norm as

(15) ${\left\|p\right\|}_{1,\bm{w}}=\sum^{n}_{i=1}w_{i}|p_{i}|.$

It is also possible to define weighted $\ell^{q}_{\bm{w}}$-norms for $q\neq
1$, but this is not needed for our investigation.

If $P=(p_{i,j})^{n,J}_{i,j=1}\in\mathbb{R}^{n\times J}$ is a matrix, we define
the $\ell^{q_{1},q_{2}}$-norms by

${\left\|P\right\|}_{q_{1},q_{2}}=\left(\sum^{n}_{i=1}\left(\sum^{J}_{j=1}|p_{i,j}|^{q_{1}}\right)^{q_{2}/q_{1}}\right)^{1/q_{2}}.$

Again by convention the $\ell^{2,0}$-‘norm’ is defined as

${\left\|P\right\|}_{2,0}=\left|\left\\{i:\sum^{J}_{j=1}|p_{i,j}|^{2}\neq
0\right\\}\right|.$

We are now able to provide the following definition:

######  Definition 2.4.

A vector $p\in\mathbb{R}^{n}$ is $s$-sparse for some $1\leq s\leq N$ if

${\left\|p\right\|}_{0}=|\mathrm{supp}(p)|\leq s.$

A collections of vectors ${p}_{1},\ldots,p_{J}\in\mathbb{R}^{n}$ is $s$-joint
sparse if

${\left\|P\right\|}_{\alpha,0}=\left|\bigcup^{J}_{j=1}\mathrm{supp}\left(p_{j}\right)\right|\leq
s,\quad\quad{{\alpha>0}}$

where $P=\left[p_{1}|\cdots|p_{J}\right]$. In typical applications $\alpha$ is
chosen to be $2$.

Based on our discussion following (14), the sparsifying operator $\mathcal{L}$
is the PA transform where $p_{j}=\mathcal{L}x_{j}$ approximates the solution
$x$ in the edge domain for each $j=1,\cdots,J$. It follows that
$P=\left[\mathcal{L}x_{1}|\cdots|\mathcal{L}x_{J}\right]$ in Definition 2.4.

If the edges were explicitly known, the prior in (14) could be tailored to
promote a posterior that not only reflects the existence of the sparse domain
but also the regions of its support. Specifically, a more accurate prior would
be dependent on $\mathrm{supp}(\mathcal{L}x)$. Although this information is
not directly available in most applications, it is possible to extract it
(approximately) from the observable data, as was demonstrated in [1, 21, 33]
using the VBJS method. In the context of this investigation, the VBJS approach
is an example of empirical Bayesian inference, and results in a support
informed sparse prior.

####  2.4.3 Variance based joint sparsity (VBJS)

To use the VBJS approach, we first require the calculation of the sample
variance $v\in\mathbb{R}^{n}$ across the rows of $P$, where
$P=[\mathcal{L}x_{1}|\dots|\mathcal{L}x_{J}]$ is the sparse matrix defined by
Definition 2.4. The (spatial) entries of $v$ are given by

(16)
${v}_{i}=\frac{1}{J}\displaystyle\sum_{j=1}^{J}{P}_{i,j}^{2}-\left(\frac{1}{J}\displaystyle\sum_{j=1}^{J}{P}_{i,j}\right)^{2},\quad
i=1,\dots,n.$

We note that because we do not have the true solution vector, $x$, we first
must approximate each $\mathcal{L}x_{j}$ from the observable MMV data $y_{j}$,
$j=1,\dots,J$. This may be done in a variety of ways, and in some cases does
not require an initial approximation of the solution itself, [33].

The VBJS technique constructs a weighting matrix $W=diag(w_{1},\dots,w_{n})$
that weighs the importance of the prior for different regions of the domain
based on the sparse representations of the data. For the compressive sensing
approach used in [1, 21], these spatially varying weights were scaled to
essentially be inversely proportional to the reciprocal of (16). That is,

$w_{i}=\frac{1}{v_{i}+\epsilon},$

where $\epsilon>0$ is some small parameter that offsets zero variance. As
discussed in [21], the weights may be further scaled to account for the
different magnitudes of the non-zero values in the sparse domain. While such
separation of scales is important for determining the spatially varying
weighting parameter in the CS environment, it is not needed in the Bayesian
framework for reasons that will become clear below. What is important is to
note that the VBJS technique’s reliance on the variance of the measurements in
the sparse domain is motivated by the provable convergence rates of edge
detection methods near and away from jumps, [4, 22, 23, 33]. In particular, if
the variance in the sparse domain is high in a region of $\mathbb{R}^{n}$,
then the corresponding entries of $W$ are small. This follows the assumption
that high variance values indicate regions where the measurements share
support in the sparse domain. Conversely, low variance values suggest truly
sparse regions, and in this case the corresponding entries of $W$ should be
large to strongly enforce the sparse prior.

For the purpose of empirical Bayesian inference, and specifically to define
the support informed sparse prior, we simply threshold the variance so that
$W$ becomes a diagonal binary mask
$M=diag(m_{1},\dots,m_{n})\in\mathbb{R}^{n\times n}$ that enforces a prior in
regions that are truly sparse, but assumes no prior in regions of support.
More specifically, for threshold $\tau=\frac{1}{n}$ we define the diagonal
entries of $M$ as

(17) $m_{i}=\begin{cases}{1},&\quad w_{i}\geq\tau\\\ {0},&\quad
w_{i}<\tau,\end{cases}$

where $w_{i}$, $i=1,\dots,n$, are provided above. Figure 1 demonstrates the
process for obtaining (17). Observe that because the weights $w_{i}$ exhibit a
large separation of scale, our masking technique is not sensitive to the
thresholding parameter $\tau$, which was chosen simply to be consistent with
the resolution of the data.

Figure 1: Demonstration of the calculation of the mask used in the proposed
support informed sparsity prior. (top-left) True function. (top-middle)
Jointly sparse vectors generated from indirect, noisy measurements $y_{j}$,
$j=1,\dots,10$, of $x$. (top-right) Spatial variance $v$ (16) across joint
sparsity vectors. (bottom-left) Weights calculated using the VBJS technique,
[1, 21]. (bottom-right) Mask (17) generated for use in the support informed
sparsity prior. The red crosses indicate locations in the spatial domain $s$
where we calculate evaluation metrics described in Section 4

Given the binary mask defined in (17), Algorithm 2 describes the VBJS process
for constructing a support informed sparsity prior as it is used for this
investigation.

Algorithm 2 Support Informed Sparsity Prior

1: Input: Set of measurement vectors, ${{Y}}=[y_{1},y_{2},\dots,y_{J}]$.

2: Output: Support informed sparsity prior.

3: Estimate the hyper-parameter $\hat{\lambda}$ using Algorithm 1.

4: Approximate the joint sparsity matrix $P$ defined in Definition 2.4 from
the $J$ measurements.

5: Compute the spatial variance $v$ in (16) from the joint sparsity matrix
$P$.

6: Use the spatial variance $v$ to determine the weighting vector $w$.

7: Use $w$ to determine the binary mask $M$ according to (17).

8: Define the support informed sparsity prior density as  (18)
$\tilde{f}_{\mathbf{X}|\bm{\hat{\lambda}},M}(x|\hat{\lambda},M)=C\exp\\{-{\hat{\lambda}}||{M}\mathcal{L}x||_{1}\\},{\quad
C>0}.$

A couple of remarks are in order:

######  Remark 2.5.

In the original construction of the weighting matrix $W$ in [21], the hyper-
parameter $\hat{\lambda}$ was intrinsically fixed. Specifically, in the CS
solution (1), the regularization parameter $\lambda$ was set to $1$.
Therefore, to separate the scales of the non-zero components in the sparse
domain, it was important that $W$ not be binary, that is, it had to vary in
magnitude to accommodate the various scales of the non-zero components in the
sparse domain. Furthermore, although not derived this way, one may understand
this scaling to be informed by the error estimate of the fidelity term, much
like what is obtained through Algorithm 1. By contrast, in the current
Bayesian framework we consider the hyper-parameter in the prior to be
independent from the support in the sparse domain. This may be advantageous in
circumstances where it is difficult to determine regions of support. In such
cases the resulting mask is close to the identity matrix, and using a hyper-
parameter determined through methods such as cross validation may still reduce
uncertainty when compared to standard priors that have fixed hyper-parameters
(which is analogous to choosing the best regularization parameter in the CS
approach). Hence we set the weighting matrix to be a binary mask matrix
indicating the support regions and learn the hyper-parameter $\hat{\lambda}$
separately, rather than combining these parameters into one as was done in
[21].

######  Remark 2.6.

The Laplace prior is commonly used to promote sparsity in the sparse domain of
the underlying signal, especially when constructing the MAP estimate of the
posterior. Analogously, $\ell_{1}$ regularization is chosen in the compressive
sensing formulation since using $\ell_{2}$ does not effectively promote
sparsity, [11]. However, when the mask matrix $M$ is applied, the sparse
regions of $\mathbf{X}$ that are to be regularized are identically zero. Hence
there is no fundamental reason to choose the $\ell_{1}$ norm. Indeed, the
$\ell_{2}$ norm was employed in the compressive sensing algorithms designed in
[14] and [21], and in both cases, the resulting algorithms more efficient.
Adopting the same rationale here, for comparative purposes, we also consider a
support informed $\ell_{2}$ prior.

###  2.5 Defining the Posterior

Armed with both a likelihood (9) and the support informed sparsity prior (18),
we are able to define the full posterior (6) for the model in (2).

Based on our previous discussion, we will consider two posterior estimates
which both use (9) and hyper-parameter $\hat{\lambda}$ calculated in Algorithm
1:

(19a)
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\hat{\lambda}|y)=C_{L}\exp\\{-{\hat{\lambda}}||\mathcal{L}x||_{1}-\frac{1}{2\sigma^{2}}||y-Ax||_{2}^{2}\\},\quad\quad
C_{L}>0,$ (19b)
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\hat{\lambda}|y)=C_{M}\exp\\{-{\hat{\lambda}}||M\mathcal{L}x||_{1}-\frac{1}{2\sigma^{2}}||y-Ax||_{2}^{2}\\},\quad\quad
C_{M}>0.$

The first uses the Laplace prior (14), while the second uses the support
informed sparse $\ell_{1}$ prior (18) constructed in Algorithm 2. Although
parameters $C_{L}$ and $C_{M}$ are not explicitly known, they are not needed
in calculation. Note that within this MMV framework, $\sigma{{}^{2}}$ can
simply be estimated as the average sample variance across all MMVs in the data
space, obtained as

${\sigma}^{2}\approx\frac{1}{n}\left(\frac{1}{J}\displaystyle\sum_{j=1}^{J}{Y}_{i,j}^{2}-\left(\frac{1}{J}\displaystyle\sum_{j=1}^{J}{Y}_{i,j}\right)^{2}\right),$

where ${Y}=\left[y_{1}|\cdots|y_{J}\right]$. For comparative purposes, we will
also consider the following posteriors,

(20a)
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\hat{\lambda}|y)=C_{L}\exp\\{-\frac{\hat{\lambda}}{2}||\mathcal{L}x||_{2}^{2}-\frac{1}{2\sigma^{2}}||y-Ax||_{2}^{2}\\},\quad\quad
C_{L}>0,$ (20b)
$\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\hat{\lambda}|y)=C_{M}\exp\\{-\frac{\hat{\lambda}}{2}||M\mathcal{L}x||_{2}^{2}-\frac{1}{2\sigma^{2}}||y-Ax||_{2}^{2}\\},\quad\quad
C_{M}>0,$

for which we use the same approximation of $\sigma^{2}$ as above. As discussed
in the introduction, the method in [5] uses a hierarchical Gaussian prior
(HGP) with the scalar hyper-priors chosen to be Gamma distributed.
Furthermore, the sparse transform operator in [5] is symmetric positive
definite. Because of the resulting conjugacy relationships, the posteriors are
known in closed form and the Gibbs sampling method can be efficiently applied,
which is a primary objective in [5]. It is important to note that employing
the mask matrix $M$ to ${\mathcal{L}}$ (even if ${\mathcal{L}}$ is constructed
to be SPD), would not yield an SPD transform. Without considering the
convergence properties of the corresponding numerical implementations, the
posterior provided in (20a) is comparable to the one in [5]. Indeed it is
suggested there that additional information may be used to determine the
hyper-prior. The $K$-fold cross validation algorithm in Algorithm 1 provides
such information. On the other hand, when the underlying signal is itself
sparse so that ${\mathcal{L}}=\mathcal{I}$ the identity matrix, (20b), which
uses the support informed $\ell_{2}$ prior, can be viewed as alternative
approach to SBL (using EM), [34, 35], since both our new algorithm and SBL
incorporate the support locations in the sparse domain to inform the prior.
Additional contrasts and comparisons are provided in the introduction. Finally
we note that while using Gibbs sampling is generally more efficient, it is not
consistent with employing the support informed prior which reduces the
uncertainty of the solution. Moreover, our method will allow downstream
processing to be more efficient, as we can more efficiently probe regions of
interest.

Finally, we point out that we can incorporate the MMV framework directly into
the likelihood estimate that appears in each of the posterior estimates, (19)
and (20). However, as we are only considering one data source given by (2),
this information is somewhat redundant. Furthermore, the added MMV information
is already accounted for in determining the support of the underlying signal.
In our numerical experiments we use $y=\frac{1}{J}\sum_{1}^{J}y_{j}$, i.e. the
mean of the observations, in each posterior estimate.

##  3 Evaluating the Posterior

Obtaining a maximum a posteriori (MAP) estimate of (19) or (20) is relatively
straightforward. However, because it provides only a mode of the desired
distribution, the MAP estimate is unable to quantify the uncertainty of the
results. Moreover, a mode of the posterior distribution may be an atypical
solution in the sense that even if it is unique, it may differ substantially
from its mean or median. By contrast, the sampling techniques described in
Section 3.2 allow for both uncertainty quantification in the signal recovery
as well as the ability to discern modes of the posterior distribution.

###  3.1 Point Estimation

While the main goal of this investigation is to draw samples of the full
posterior, it is instructive to first demonstrate how the MAP point estimate
is used to recover a mode of the unknown distribution, [27].

For any realization of (2), $y_{j}\sim\mathbf{Y}$, $j=1,\dots,J$, the MAP
estimate for either posterior in (19) or (20) is computed as

(21)
$x_{MAP}^{j}\equiv{\underset{x\in\mathbb{R}^{n}}{\operatorname{argmax}}}\left\\{\log\left[\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\hat{\lambda}|y_{j})\right]\right\\},\quad
j=1,\dots,J.$

A corresponding convex optimization problem is then easily derived for the
posteriors in (19) or (20) as

$\displaystyle x_{MAP}^{j}$
$\displaystyle={\underset{x\in\mathbb{R}^{n}}{\operatorname{argmax}}}\Big{\\{}\log(C)-\Big{(}\frac{\hat{\lambda}}{\nu}||\mathcal{T}x||_{\nu}^{\nu}+\frac{1}{2\sigma^{2}}||y_{j}-Ax||_{2}^{2}\Big{)}\Big{\\}}$
(22)
$\displaystyle={\underset{x\in\mathbb{R}^{n}}{\operatorname{argmin}}}\Big{\\{}\frac{\hat{\lambda}}{\nu}||\mathcal{T}x||_{\nu}^{\nu}+\frac{1}{2\sigma^{2}}||y_{j}-Ax||_{2}^{2}\Big{\\}},\quad
j=1,\dots,J,$

where $\nu=1$ or $2$, $\mathcal{T}=\mathcal{L}$ or $\mathcal{ML}$, and
$C=C_{L}$ or $C_{M}$. Since the MAP estimation can be reduced to an
optimization problem, (22) is often understood from a non-Bayesian
perspective. In particular, note that the MAP estimate used in Step 7 of
Algorithm 1 is similarly constructed, with $\hat{\lambda}$ given by
$\tilde{\lambda}$ in Step 6, $\mathcal{T}=\mathcal{L}$, and $\nu=1$. We will
also use the average MAP estimate given by (22) over the $J$ measurement
vectors to set the initial state $x^{0}$ for the Metropolis-Hastings Algorithm
given in Algorithm 3. Finally, although not discussed in this current
investigation, we point out that there may be situations for which there are
regions in a signal domain where efficiency is more important than quantifying
uncertainty. In such cases using (22) to construct an approximation within
these regions may be more appropriate.

###  3.2 Drawing samples of the posterior distribution

We seek to sample the posterior probability distribution of $\mathbf{X}$
described by (19) or (20) for the model given by (2). This allows us both to
quantify the uncertainty in our signal recovery as well as to discern modes of
the posterior distribution that are unavailable in a MAP estimate. By Bayes’
rule, (4), it is clear that sampling $\mathbf{X}$ in (2) from observations in
(3) requires an approximation to the probability density $f_{\mathbf{Y}}(y)$,
which typically does not have a closed form solution. Moreover, numerical
quadrature rules are impractical whenever the dimension of the sample space
$\mathbb{R}^{n}$ is high [25]. We therefore choose to use the Markov Chain
Monte Carlo (MCMC) method which is designed to construct an ergodic Markov
chain that converges to the stationary probability density
${f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)$ in (19), which is
proportional to the posterior we seek. Overviews of MCMC may be found in [20,
30, 25]. Because our approach does not assume conjugate priors, the posterior
density (4) we derive does not represent a known posterior distribution. Hence
we employ the particular MCMC method known as the Metropolis-Hastings (MH)
algorithm [24, 26], which we briefly describe below.

####  3.2.1 The Metropolis Hastings (MH) Algorithm

In the MH algorithm, a Markov chain in ${\Omega}$ is constructed so that it
converges to the conditional posterior density
${f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)$ in (4). The MH
algorithm first defines a _proposal distribution_
$q:\Omega\times\Omega\rightarrow[0,1]$ such that when the Markov chain is at
state $x^{k-1}$, it proposes a _candidate_ state $x^{cand}\sim q$. Its
acceptance of the next state, $x^{k}$, is determined by an _acceptance
probability_ $\alpha:\Omega\times\Omega\rightarrow[0,1]$ given by

$\displaystyle\alpha(x^{cand}|x^{k-1})$
$\displaystyle=\min\Bigg{\\{}1,\frac{q(x^{k-1}|x^{cand}){f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x^{cand},\lambda|y)}{q(x^{cand}|x^{k-1}){f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x^{k-1},\lambda|y)}\Bigg{\\}}$
(23)
$\displaystyle{\approx}\min\Bigg{\\{}1,\frac{q(x^{k-1}|x^{cand})\hat{f}_{\mathbf{X},\bm{{{\lambda}}}|\mathbf{Y}}(x^{cand},{\hat{\lambda}}|y)}{q(x^{cand}|x^{k-1})\hat{f}_{\mathbf{X},\bm{{{\lambda}}}|\mathbf{Y}}(x^{k-1},{\hat{\lambda}}|y)}\Bigg{\\}},$

where $\hat{f}_{\mathbf{X},{\bm{\lambda}}|\mathbf{Y}}$ is defined in (19). In
this way the MH algorithm is able to evaluate the full posterior (4) without
needing to calculate the evidence $f_{\mathbf{Y}}(y)$ in (5). The general MH
algorithm is summarized in Algorithm 3.

Algorithm 3 The Metropolis-Hastings Algorithm

1: Input: The initial state $x^{0}$, the choice of proposal distribution $q$,
and its corresponding variance $\sigma_{q}$.

2: Output: Markov chain $\\{x^{1},x^{2},x^{3},\ldots,x^{M}\\}$.

3: Initialize $x^{0}$ and $k=1$.

4: while $k\leq M$ do

5: Propose $x^{cand}\sim q(x^{k}|x^{k-1})$.

6: Compute $\alpha(x^{cand}|x^{k-1})$ using (23) or (24).

7: Sample $u\sim\mathcal{U}[0,1]$.

8: if $u<\alpha$ then

9: Accept the candidate: set $x^{k}=x^{cand}$.

10: else

11: Reject the candidate: set $x^{k}=x^{k-1}$.

12: end if

13: end while

Algorithm 3 requires user inputs of the initial state $x^{0}$ and the proposal
distribution $q$. For our experiments we set the initial state $x^{0}$ as the
average MAP estimate in (22) over all measurement vectors $y_{j}$,
$j=1,\dots,J$, and $q$ to be Gaussian, such that
$x^{cand}\sim\mathcal{N}(x^{k-1},\sigma_{q})$. As noted previously, we also
use $y=\frac{1}{J}\sum_{j=1}^{J}y_{j}$, the mean of the observable data. Since
the Gaussian distribution is symmetric, i.e.
$q(x^{k-1}|x^{cand})=q(x^{cand}|x^{k-1})$, the acceptance probability (23) is
simplified to

(24)
$\alpha(x^{cand}|x^{k-1})=\min\Bigg{\\{}1,\frac{\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x^{cand},\hat{\lambda}|y)}{\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(x^{k-1},\hat{\lambda}|y)}\Bigg{\\}}.$

While the ideal acceptance ratio is an open question for the Metropolis-
Hastings algorithm, a typical range for an acceptance ratio is $[.2,.8]$. This
enables the method to “reasonably” explore the space, preventing the chain
from getting stuck in one place or moving too quickly from one component to
another, [30]. Because the acceptance ratio explicitly depends on
$\sigma_{q}$, for the purposes of this investigation, we simply chose
$\sigma_{q}$ so that it fell within this interval. Specifically, for all
posteriors given in (19) and (20), we ran a few iterations of Algorithm 3 and
observed the behavior of the acceptance ratio. If it fell outside the
interval, we adjusted the input parameter $\sigma_{q}$ accordingly.

It is important to note that $x^{0}$ may not provide a good initial guess.
This is a standard problem in all MCMC algorithms and is addressed by
implementing what is known a burn-in period. That is, if $M$ is the total
number of iterations computed by Algorithm 3, the first $B<M$ (called the
burn-in rate) iterations are discarded, so that the total length of the chain
is $M-B$. We chose $M=50,000$ and $B=25,000$ in all of our experiments, the
same values used in studies in which convergence is consistently observed,
[19]. In this regard we point out that our strategy for choosing $x^{0}$, and
more importantly our use of the support informed sparse prior, allowed us to
pick smaller burn-in rate and chain length while still reducing the
uncertainty in all of our experiments. More studies are needed to confirm this
is the case in general.

The iterations of the Markov chain obtained using Algorithm 3 may be
represented by a matrix $\mathcal{M}\in\mathbb{R}^{n\times M}$ given by

(25)
$\mathcal{M}=\begin{bmatrix}x_{1}^{0}&x_{1}^{1}&x_{1}^{2}&\dots&x_{1}^{M}\\\
x_{2}^{0}&x_{2}^{1}&\dots&x_{2}^{3}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
x_{n}^{0}&x_{n}^{1}&x_{n}^{2}&\dots&x_{n}^{M}\end{bmatrix},$

where each column of the matrix $\mathcal{M}$ represents one state in the
Markov chain. In particular, an entry $x^{t}_{i}$ is the $i^{th}$ component of
the $t^{th}$ state in the Markov chain. As in [5], we then calculate a final
estimate for the unknown $x$ at each component as the mean of the samples
after discarding the first $B$ states in the Markov chain:

(26) $x_{i}=\frac{1}{M-B}\sum_{t=B+1}^{M}x_{i}^{t}.$

A summary of our proposed empirical Bayesian inference approach to solving the
statistical inverse problem defined in (2) is provided in Algorithm 4. The
mean in (26) is obtained after discarding the first $B<M$ burn-in iterations.
The Algorithm is written for the support informed sparse $\ell_{1}$ prior, but
the $\ell_{2}$ prior case given in (20b) follows analogously.

Algorithm 4 Support Informed Sparse Prior Empirical Bayesian Inference
Algorithm

1: Input: Multiple measurement vectors, $y_{j}$, for $j=1,...,J$.

2: Output: Markov chain $\mathcal{M}$ in (25) describing samples from
$f_{\textbf{X},\bm{\lambda}|\mathbf{Y}}(x,\lambda|y)$.

3: Estimate the hyper-prior $\hat{\lambda}$ according to Algorithm 1.

4: Calculate the support informed sparsity prior using Algorithm 2.

5: Sample the posterior (19b) using Algorithm 3.

##  4 Numerical Experiments

We are now ready to demonstrate the utility of our support informed sparsity
prior (18) in sampling the full posterior (19b) and the posterior derived
using the support informed $\ell_{2}$ prior used in (20b). In analyzing our
results, we will calculate the relative error

(27) $\frac{||x-x^{true}||}{||x^{true}||}.$

Here $x$ is determined by (26) and $x^{true}$ is the underlying signal. We
will also consider measurements that are acquired with increasing levels of
noise. The signal to noise ratio (SNR) computed in decibals (dB) according to

(28)
$\text{SNR}=10\log_{10}\left(\frac{\mathbb{E}[\mathbf{Y}]}{\sigma^{2}}\right),$

where $\mathbb{E[\mathbf{Y}]}$ is the expected value of the collected data and
$\sigma^{2}$ is the variance of the complex additive Gaussian white noise
(AGWN) in (2).

To quantify the uncertainty of our recovery, we approximate the credibility
interval (also called the Bayesian interval) of the resulting distribution,
[27, 25, 12], which is defined for $x\in\mathbb{R}$ as:

######  Definition 4.1.

Given $\eta\in(0,1)$, the credibility set $C_{\eta}=(\ell,u)$ is an
equiprobability hypersurface that contains $100(1-\eta)$% of the mass of the
posterior distribution and is defined according to

$\int_{C_{\eta}}\hat{f}_{\mathbf{X},\bm{\lambda}}(x,\hat{\lambda}|y)dx=1-\eta,\quad\hat{f}_{\mathbf{X},\bm{\lambda}}(x,\hat{\lambda}|y)|_{x\in\partial
C_{\eta}}=\text{constant}.$

Following Definition 4.1, we calculate $C_{\eta}$ as

(29) $C_{\eta}=\\{(\ell,u)\hskip 2.84526pt|\hskip
2.84526pt\hat{f}_{\mathbf{X},\bm{\lambda}|\mathbf{Y}}(\ell\leq x\leq
u,\hat{\lambda}|y)=1-\eta\\}.$

For example, $C_{0.05}$ is a $95\%$ credibility interval of a probability
density and is interpreted as a sub-interval of the domain that contains
$95\%$ of the probability density. We note that $C_{\eta}$ is generally not
unique because there is typically more than one interval that satisfies (29).
This fact emphasizes the importance of initiating the MCMC chain with an
appropriate point estimate, so that we may constrain our recovered credibility
interval to contain the mode, [27]. Hence, as noted previously, we choose the
beginning state of the chain $x^{0}$ in Algorithm 3 to be the average MAP
estimate (22) over all measurements $y_{j}$ for $j=1,\dots,J$, which is
approximately the peak of the mode of the posterior probability density. In
this way we ensure that the Markov chain explores a meaningful region of the
posterior density, [30].

We also analyze the convergence of our MCMC chains to steady state through the
use of autocorrelations, trace plots and acceptance ratios. Trace plots show
the state of the chain over each iteration. Theoretically, if the Markov chain
has converged to the posterior density, the states of the chain should not be
correlated. Therefore, the chain is said to be well-mixed if the trace plot
exhibits behavior indicating that the states are independently and identically
distributed (iid). On the other hand, any portion of the Markov chain having
high auto-correlation between states would imply that the Markov chain has not
(yet) converged to the stationary distribution. As another indicator of
mixing, we use the Autocorrelation Function (ACF) package in Matlab, [29], to
compute the correlogram of each Markov Chain recovered through the proposed
MCMC sampling techniques.

The acceptance ratio determined from line 9 in Algorithm 3 is the number of
times $x^{cand}$ is accepted as the next sample relative to the total number
of iterations. Therefore, a large acceptance ratio (typically $\alpha>.8$) may
indicate the Markov chain is stuck in a certain region of the posterior
probability density and has not explored the entire support of the density. In
contrast, a small acceptance ratio (e.g. $\alpha<.2$) may indicate that the
chain may be moving too quickly within the domain of the proposal density.
Further, with a low acceptance ratio, it remains possible that the chain does
not explore certain isolated modes of the posterior probability density [30].
While the ideal acceptance ratio remains an open question for the Metropolis-
Hastings algorithm, as a general rule of thumb, it should be close to $.5$ for
a one-dimensional problem [30]. This metric was primarily derived for the
Gaussian proposal distribution, which is also used in this paper.

###  4.1 Sparse Signal Recovery in One Dimension

We first consider the recovery of a sparse signal. In this case the
sparsifying transform matrix operator ${\mathcal{L}}={\mathcal{I}}$ in (19)
and (20). This example eliminates any complications in evaluating the
posterior caused by different choices of sparsifying transform operators.
Moreover, sparse signal recovery is perhaps the simplest case to analyze,
since conjugate priors are easily obtained, and indeed there is a closed form
posterior for (20a).777As noted previously, the focus of our investigation is
not to compare the efficiency of sampling methods for closed form posteriors.
Hence we will still employ the MH algorithm even though Gibbs sampling is more
appropriate in this particular case. As noted previously, in the case of
sparse signals, SBL with EM yields a comparable posterior, in the sense that
it uses an empirically based support informed prior, to that in (20b). It is
not computationally efficient in higher dimensions, however, [15, 34, 35]. To
further simplify the analysis we choose the non-zero values of the underlying
signal to be exactly $1$.

Figure 2: The mean and $C_{0.05}$ credibility intervals when the input data
are noisy Fourier samples with $\text{SNR}=0$dB. (left) comparison of (19a)
and (19b); (right) comparison of (20a) and (20b).

Figure 2 displays the results and the $C_{0.05}$ credibility intervals where
the observations in (2) are noisy Fourier samples with $\text{SNR}=0$dB, so
that $A$ is modeled by the discrete Fourier transform, and the complex
additive white Gaussian noise is given by
$\mathbf{E}\sim\mathcal{CN}(0,\sigma^{2}\mathbb{I}_{{n}})$ with $\sigma\approx
4.5$. Other credibility intervals (e.g. $C_{0.01}$ and $C_{0.1}$) for
$\text{SNR}=0$dB produce comparable results, with reduced uncertainty
especially apparent in the sparse regions of the signal when our new support
informed prior is used. As shown respectively in Figure 2(left) and (right),
this is true for both (19b) and (20b), although the results are more striking
for the sparse prior ($\ell_{1}$) case.

Figure 3: Relative error of the mean (27) for $\text{SNR}\in[-10,10]$ at
(left) a location where the sparse signal equals zero and (right) at a
location where the sparse signal is non-zero.

Figure 3 compares the relative error of the mean (27) using (19a), (19b),
(20a) and (20b) for decreasing SNR. These comparisons are made at locations in
the domain where the signal value is zero (left) and where it is $1$ (right).
From Figure 3(right) we observe that by using the support informed prior, that
is by “masking out” the sparse prior at the non-zero function locations, the
target density is equivalent to likelihood in regions of support, i.e. we are
not enforcing an incorrect sparse prior in those locations. This is true for
all levels of noise. The relative error in sparse regions, depicted in Figure
3(left), demonstrates comparable results for $\text{SNR}\geq 0$. Once the
noise dominates the signal, the mean of the posterior does not yield
meaningful results using either prior.888In this case it is probably better to
employ (19a) using a Gibbs sampler, [30]. For low SNR, the results using the
standard (not support informed) priors are essentially the same regardless of
whether or not the signal is zero or non-zero. When combining these results
with those in Figure 2, we see that even though the support informed sparse
prior may not be as accurate in the mean when the signal value is zero, the
uncertainty is still reduced in those regions. Thus we conclude that using the
support informed sparse prior yields improved accuracy in the mean at the non-
zero signal locations while also reducing the uncertainty everywhere.

###  4.2 Signal Recovery in One Dimension

We now consider reconstructing the signal defined as

######  Example 4.2.

$x(s)=\begin{cases}40,&\quad 0.1\leq s\leq 0.25\\\ 10,&\quad 0.35\leq s\leq
0.325\\\ 20\sqrt{2\pi}e^{-\left(\frac{s-0.75}{0.05}\right)^{2}},&\quad
s>0.5\\\ 0,&\quad\mathrm{otherwise,}\end{cases}$

where $s\in[0,1]$ is the domain of interest. Similar to the example used in
[5], the domain is uniformly discretized so that $s_{i}=\frac{i}{n}$,
$i=1,\dots,n$ with $n=80$. For all of our experiments in one dimension, we
simulate $J=20$ measurements. We choose $20$ since it was shown to be a
sufficient number of measurements for reconstructing a signal using the VBJS
method in [21]. We also choose $K=20$ and $M=10$ for Algorithm 1.

Because we are ultimately interested in imaging applications such as MRI, SAR
and ultrasound, in our experiments we assume that the observable data $y$
consists of multiple, noisy Fourier measurements of $x$. That is, each
measurement is collected according to (2), with $A$ defined as the discrete
Fourier transform matrix and $\sigma>0$. We point out that the discrete
Fourier transform matrix is better conditioned than the blurring model used in
[5]. Here we allow for the measurements to be collected with significantly
higher noise levels, i.e., $\sigma>>0$ or SNR $<<1$, however. In what follows
we show how results for $SNR=5$ dB and $SNR\approx-.2$ dB.

Figure 4: Approximate reconstructions (26) of Example 4.2 displayed along with
credibility intervals $C_{0.05}$ when using (left) the Laplace prior (19a)
(middle) the proposed support informed sparsity prior (19b) and (right) the
proposed support informed prior described in (20b). (top) $\sigma\approx 5.5$
($SNR=5$ dB), (bottom) $\sigma=10$ ($SNR\approx-.2$ dB).

Figure 4 displays the mean of the Markov chain (26) recovered on each grid
point in the domain when sampling the posteriors described in (19) with (left)
the Laplace prior (14) and (right) the proposed support informed sparsity
prior (18) for (top) $\sigma\approx 5.5$ ($SNR=5$ dB) and (bottom) $\sigma=10$
($SNR\approx-.2$ dB). Along with the mean, in Figure 4 we also display
$C_{0.05}$ for each reconstruction. Notice that the proposed support informed
sparsity prior results in tighter credibility intervals and more accurate
reconstructions, especially in the $\sigma=10$ case.

Figure 5: Estimated posterior distributions at the points in the domain
indicated by the red crosses in Figure 1. Table 1 reports each distribution
and its corresponding parameters. (top) $\sigma\approx 5.5$ ($SNR=5$ dB);
(bottom) $\sigma=10$ ($SNR\approx-.2$ dB).

Figure 5 shows the resulting distributions at the points in the domain
indicated by the red crosses in Figure 1. Each distribution is estimated given
all samples of the Markov chain at the desired point in the domain after
discarding the first $B$ states. Distributions are optimized based on the
Pearson Chi-Squared goodness of fit metric and the mean-squared error when
compared to the histogram. The resulting distributions and their corresponding
parameters are displayed in Table 1. Consistent with the remainder of our
results, we see that the distributions of the support informed $\ell_{1}$
prior at edge locations is more closely centered around the true function
values. In smooth regions, we see compelling evidence to use the support
informed $\ell_{2}$ prior. This motivates future work into exploring the use
of different priors in different regions of the domain.

Posterior | Noise | Location 1 | Location 2 | Location 3 | Location 4
---|---|---|---|---|---
(19a) | 5 dB | Beta | $\chi^{2}$ | Normal | Log-Normal
(12.45,11.94) | (247.60) | (0.46, 1.50) | (0.01)
(19b) | 5 dB | $\chi^{2}$ | Pearson type III | Normal | Log-Normal
(95.34) | (-0.07) | (-1,1.29) | (0.06)
(20b) | 5 dB | Beta | Beta | Normal | Normal
(22.15,18.08) | ($5.13\times 10^{6}$, 196.99) | (0.15, 1.61) | (48.48, 1.53)
(19a) | -.2 dB | Beta | Weibull (min) | Beta | Pearson type III
(23.53,34.94) | (3.61) | (399.66,$1.15\times 10^{7}$) | (-0.11)
(19b) | -.2 dB | Log-Normal | Weibull (max) | Log-Normal | Weibull (min)
(0.03) | (3.32) | (0.01) | (3.2)
(20b) | -.2 dB | Beta | Beta | Weibull (max) | Beta
(325.45,$3.32\times 10^{3}$) | (16.56,19.12) | (3.9) | (48.37,44.73)

Table 1: The distributions and corresponding parameters resulting from a best-
fit analysis of the last $M-B$ Markov chain samples (25) at the locations in
the domain indicated by the red crosses seen in Figure 1. The distributions
along with the histogram of samples are visualized in Figure 5.

######  Remark 4.3.

We do not report on results using the posterior given by (20a) for Example 4.2
because the MH algorithm as implemented in Algorithm 3 did not converge in
this case. We speculate that the (unmasked) $\ell_{2}$ norm in the prior is
prohibitively large for numerical computation. It may be possible to calculate
(23) so as to avoid this issue, for example by using logarithms.
Alternatively, one could use a Gibbs sampler with an appropriate SPD operator
to replace ${\mathcal{L}}$ in (20a), see for example the discussion in [5].
Since our previous observations in the sparse signal case suggest that the
resulting confidence interval would still be larger than those displayed in
Figure 4, and as comparing different sampling algorithms is not the focus of
this investigation, we give no further consideration to (20a).

Figure 6: Acceptance ratios resulting from sampling the posterior defined
using the Laplace prior (19a), support informed sparsity prior (19b) and
support informed $\ell_{2}$ prior (20b). (left) $\sigma\approx 5.5$ ($SNR=5$
dB); (right) $\sigma=10$ ($SNR\approx-.2$ dB).

We now evaluate the mixing and convergence of the MCMC chains. Figure 6
displays the acceptance ratio across all MH iterations. With our strategy to
adjust the variance of the proposal distribution, we see that each chain
exhibits an acceptable acceptance ratio, with the acceptance ratio for the
support informed sparsity prior tending closer to the desirable $.5$ value. We
note that we did not try to optimize $\sigma_{q}$ in Algorithm 3, and that its
initial tuning as described in the paragraph below Algorithm 3 was done in the
same way for each choice of posterior. Figure 7 displays the correlograms
associated with sampling the posterior distributions at each grid point marked
with a red cross in Figure 1. The results affirm that prior information about
the support locations in the sparse domain, which is achievable $SNR=5$ dB,
does indeed reduce the uncertainty. When comparing the correlograms resulting
from the posteriors in (19a), (19b), and (20b) for lower SNR, it is evident
that the proposed support informed priors yield less or similarly correlated
samples near discontinuities than the more standard Laplace prior does.
Consistent with the results in the sparse signal case, the Laplace prior
demonstrates earlier mixing in smooth regions than the support informed sparse
prior, but the mixing is comparable to that achieved using the support
informed $\ell_{2}$ prior. These results, confirmed in Figure 10, suggest that
as the SNR decreases a better result overall may be obtained when the choice
of posterior is itself domain dependent; that is, when the support informed
$\ell_{2}$ prior and the support informed sparse prior are adopted in
different regions once the support in the sparse domain has been approximated.
This will be discussed in future investigations, as such a “splitting” may
also yield more efficient computation.

Figure 7: The autocorrelation function (ACF) calculated for various lag
lengths at the points in the domain indicated by the red crosses in Figure 1.
Samples are taken from the posterior distributions defined using the Laplace
prior (19a), support informed sparsity prior (19b) and support informed
$\ell_{2}$ prior (20b). (top) $\sigma\approx 5.5$ ($SNR=5$ dB), (bottom)
$\sigma=10$ ($SNR\approx-.2$ dB).

Figure 8 plots the value of each sample drawn from the posterior distribution
determined via Algorithm 3 described by (top) (19a), (middle) (19b), and
(bottom) (20b) at the values in the domain indicated by the red crosses in
Figure 1 for $SNR=5$ dB ($\sigma\approx 5.5$). While each trace plot yields
similar behavior, it is evident that at all points considered, using the
posterior defined with the support informed sparsity prior (19b) results in
samples that are more uniformly distributed about the true function value,
indicated by the horizontal red line. We note that comparable results were
obtained for $\sigma=10$ ($SNR\approx-.2$ dB) and are not displayed here.

Figure 8: Trace plots calculated at the points described by the red crosses in
Figure 1. Samples are taken from the posterior distributions defined using the
(top) Laplace prior (19a); (second row) support informed sparsity prior (19b);
and (bottom) support informed $\ell_{2}$ prior (20b). Here $\sigma=10$
($SNR\approx-.2$ dB). The trace plots display similar behavior when
$\sigma\approx 5.5$ ($SNR=5$ dB).

As was done for the sparse signal recovery case, we also evaluate the accuracy
of the reconstructions (26) for Example 4.2 given discrete Fourier
measurements acquired with increasing levels of complex additive Gaussian
white noise in (2). Figure 9 displays the relative error, (27), and Figure 10
displays the absolute error, $|x_{i}-x^{true}_{i}|$, at the points in the
spatial domain described by the red crosses in Figure 1. Each point $x_{i}$
across the domain is computed by averaging the results over 20 independent
MCMC trials. Once again we see that the support informed priors consistently
yield either better or comparable results for all levels of SNR. This is true
in smooth regions as well as near discontinuities, in which case it is clear
that using the Laplace prior in regions where the underlying signal is not
sparse in the presumably sparse domain yields dramatically worse results.

Figure 9: The relative reconstruction error (27) calculated across the entire
spatial domain for decreasing SNR (28). Figure 10: Absolute error comparison
for decreasing SNR (28) at each of the red crosses in Figure 1

##  5 Conclusions

This paper developed a new empirical Bayesian inference method using joint
sparsity. In particular, the joint sparsity of multiple measurement vectors
(MMVs) of observable data in the sparse domain is used to determine the likely
support in the sparse domain. From this information a new support informed
sparse prior is developed, which is then used to generate a posterior and
subsequently sampled using the MCMC algorithm. The support in the sparse
domain is obtained using provable convergence rates of edge detection methods
from the observable data (in this case Fourier data) near and away from that
support, [2, 4, 22, 23, 33].

Our numerical experiments for sparse signals and piecewise smooth functions
demonstrate that the new support informed sparse prior yields reduced
uncertainty when compared to the Laplace prior, and improved accuracy in the
mean wherever the underlying signal was non-zero in the sparse domain. These
findings strongly suggest that employing a any global (unmasked) sparse prior,
regardless of the type, leads to less accurate results. We also observed
faster convergence in the ACF plots in smooth regions and near
discontinuities. This was more apparent for signals with higher SNR, but still
true when $SNR\approx 0$ dB.

The support informed sparse prior does not yield a closed form solution for
the posterior, so we therefore use the MH algorithm to sample it. This is of
course not as efficient as Gibbs sampling. However, because our method samples
at each individual point in the domain, meaning that separate chains are
independently constructed for each pixel, extensions to multi-dimensional
images is straightforward. Our approach may also be expanded to include
complex valued signals and images. In this case, it was shown in [3, 14, 31,
32] how the PA transform may be efficiently constructed in multi-dimensions to
promote sparsity in the real part of the respective signal or image. Finally,
a main advantage of our technique is that once the support in the sparse
domain is determined, different priors can be used in each smooth region, and
in some cases a more efficient MAP estimate will suffice. As a main motivation
for this investigation is to reduce the uncertainty in applications such as
SAR imaging, we will explore these ideas in future work.

Acknowledgments. Anne Gelb’s work is supported in part by the NSF grants DMS
#1502640 and DMS #1912685, AFOSR grant #FA9550-18-1-0316, and ONR MURI grant
#N00014-20-1-2595. Theresa Scarnati’s work is supported in part by AFOSR LRIR
#18RYCOR011.

##  Appendix A Polynomial Annihilation (PA) Matrix Construction

While a variety of sparsifying operators can be used for (19), for this
investigation we employ the polynomial annihilation (PA) operator,
${\mathcal{L}}={\mathcal{L}}^{m}$, which is based on the PA edge detection
method, [4]. The primary advantage in using the PA operator is that it is
constructed to yield higher order approximations whenever $m>1$. (When $m=1$,
it is equivalent to TV.) This means that the approximation of the
corresponding edge function of a piecewise smooth function is more likely to
be sparse, and in particular yields values $\approx 0$ whenever the
corresponding grid point falls in a smooth region of the function. This is
especially important when the smooth region has a lot of variation (as opposed
to being piecewise constant) or if the measurements are sparsely sampled.
Below we describe the construction of ${\mathcal{L}}^{m}$. More information
can be found in [3].

For ease of presentation and for consistency with our manuscript we consider
$f(x)$ to be a piecewise smooth function on $[-\pi,\pi]$ (the finite domain is
arbitrary) given on uniform grid point values $x_{j}=-\pi+\frac{2\pi j}{n}$,
$j=1,\cdots,n$. We note that this is not a requirement for the PA method, and
furthermore that the method is designed to approximate the edge domain for
$d>1$ dimensional piecewise smooth functions for any given set of scattered
points of data, [4]. However, since images are usually defined on a Cartesian
grid, for the purposes of defining a sparsifying transform operator to be used
in the prior given by (14), as well as all of the priors defined in sequel
throughout the manuscript, it was demonstrated in [3] that applying the PA
transform operator dimension by dimension was accurate and efficient.

With this in mind we now define the $m$th order PA edge detection
approximation as

(30) $L^{m}f(x)=\frac{1}{q^{m}(x)}\sum_{x_{j}\in S_{x}}c_{j}(x)f(x_{j}),$

where $S_{x}$ is the local set of $m+1$ grid points from the set of given grid
points about $x$, $c_{j}(x)$ are the polynomial annihilation edge detection
coefficients, (31), and $q^{m}(x)$ is the normalization factor, (32). Each
parameter of the method can be further described as:

  * •

$S_{x}$: For any particular cell $I_{j}=[x_{j},x_{j+1})$, there are $m$
possible stencils, $S_{x}$ of size $m+1$, that contain the interval $I_{j}$.
For simplicity, we assume that the stencils are centered around the interval
of interest, $I_{j}$, and are given by

$S_{I_{j}}=\\{x_{j-\frac{m}{2}},\cdots,x_{j+\frac{m}{2}}\\},\hskip
14.45377ptS_{I_{j}}=\\{x_{j-\frac{m+1}{2}},\cdots,x_{j+\frac{m-1}{2}}\\}$

for $m$ even and odd respectively. For non-periodic or non-zero-padded
solutions the stencils are adapted to be more one sided as the boundaries of
the interval are approached, [4]. To avoid cumbersome notation, we write
$S_{x}$ as the generic stencil unless further clarification is needed.

  * •

$c_{j}(x)$: The polynomial annihilation edge detection coefficients,
$c_{j}(x),\;j=1,\dots,m+1$, are constructed to annihilate polynomials up to
degree $m$. They are obtained by solving the system

(31) $\sum_{x_{j}\in
S_{x}}c_{j}(x)p_{\ell}(x_{j})=p_{\ell}^{(m)}(x),\;j=1,\dots,m+1,$

where $p_{\ell}$, $\ell=0,\dots,m,$ is a basis of for the space of polynomials
of degree $\leq m$.

  * •

$q^{m}(x)$: The normalization factor, $q^{m}(x)$, normalizes the approximation
to assure the proper convergence of $L^{m}f$ to the jump value at each
discontinuity. It is computed as

(32) $q^{m}(x)=\sum_{x_{j}\in S_{x}^{+}}c_{j}(x),$

where $S_{x}^{+}$ is the set of points $x_{j}\in S_{x}$ such that $x_{j}\geq
x$.

Since the data are given on uniform points, there is an explicit formula for
the PA edge detection coefficients, independent of location $x$, computed as
([4])

(33) $c_{j}=\frac{m!}{\prod_{k=1,k\neq j}^{m+1}(j-k)\Delta x},\quad
j=1,\dots,m+1,\quad\Delta x=\frac{2\pi}{n}.$

We can now define the PA transform matrix ${\mathcal{L}}^{m}$ used in the
sparse prior in terms of its elements as

(34) ${\mathcal{L}}^{m}_{j,l}=\frac{c(j,l)}{q^{m}(x_{l})},\hskip
14.45377pt0\leq l\leq n,\hskip 14.45377pt0\leq j<n,$

where

$c(j,l)=\left\\{\begin{array}[]{lr}c_{j-l-\lfloor\frac{m}{2}\rfloor},&0<j-l-\lfloor\frac{m}{2}\rfloor+s(j,l)\leq
m+1\\\ 0&\textrm{otherwise}\end{array}\right.$

and

$s(j,l)=\left\\{\begin{array}[]{lr}l-\lfloor\frac{m}{2}\rfloor,&\textrm{$l\leq\lfloor\frac{m}{2}\rfloor$}\\\
l+m-\lfloor\frac{m}{2}\rfloor-n,&\textrm{$l+m-\lfloor\frac{m}{2}\rfloor>n$}\\\
0&\textrm{otherwise.}\end{array}\right.$

If the underlying image is known to be periodic, or is zero padded at the
boundaries, a centered stencil can be used throughout the domain producing a
circulant matrix ${\mathcal{L}}^{m}$. A reduction of accuracy is expected near
the boundaries for the non-periodic case due to the one-sided stencils.

For example, assuming periodicity, the banded matrix ${\mathcal{L}}^{1}$ and
${\mathcal{L}}^{3}$ are given by

(35) $\mathcal{L}^{1}=\begin{bmatrix}1&-1&&&\\\ &1&-1&&\\\ &&\ddots&\ddots&\\\
&&&1&-1\\\
-1&&&&1\end{bmatrix},\quad\mathcal{L}^{3}=\frac{1}{2}\begin{bmatrix}3&-3&1&&&-1\\\
-1&3&-3&1&&\\\ &\ddots&\ddots&\ddots&\ddots&\\\ &&1&3&-3&1\\\ 1&&&1&3&-3\\\
-3&1&&&-1&3\end{bmatrix}.$

## References

  * [1] B. Adcock, A. Gelb, G. Song, and Y. Sui, Joint sparse recovery based on variances, SIAM Journal on Scientific Computing, 41 (2019), pp. A246–A268, https://doi.org/10.1137/17M1155983, https://doi.org/10.1137/17M1155983, https://arxiv.org/abs/https://doi.org/10.1137/17M1155983.
  * [2] R. Archibald and A. Gelb, Reducing the effects of noise in image reconstruction, Journal of Scientific Computing, 17 (2002), pp. 167–180.
  * [3] R. Archibald, A. Gelb, and R. Platte, Image reconstruction from undersampled fourier data using the polynomial annihilation transform, Journal of Scientific Computing, (2015), https://doi.org/10.1007/s10915-015-0088-2.
  * [4] R. Archibald, A. Gelb, and J. Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, 43 (2005), pp. 259–279.
  * [5] J. M. Bardsley, Mcmc-based image reconstruction with uncertainty quantification, SIAM Journal on Scientific Computing, 34 (2012), pp. A1316–A1332, https://doi.org/10.1137/11085760X, https://doi.org/10.1137/11085760X, https://arxiv.org/abs/https://doi.org/10.1137/11085760X.
  * [6] A. Bhadra, J. Datta, N. G. Polson, B. Willard, et al., Lasso meets horseshoe: A survey, Statistical Science, 34 (2019), pp. 405–427.
  * [7] D. Calvetti, M. Pragliola, and E. Somersalo, Hybrid solver for hierarchical Bayesian inverse problems, arXiv preprint arXiv:2003.06532, (2020).
  * [8] D. Calvetti, M. Pragliola, E. Somersalo, and A. Strang, Sparse reconstructions from few noisy data: analysis of hierarchical Bayesian models with generalized gamma hyperpriors, Inverse Problems, 36 (2020), p. 025010.
  * [9] E. J. Candes, J. Romberg, and T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory, 52 (2006), pp. 489–509.
  * [10] 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.
  * [11] E. J. Candès, J. K. Romberg, and T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Communications on Pure and Applied Mathematics, 59 (2006), pp. 1207–1223, https://doi.org/10.1002/cpa.20124, https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.20124, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.20124.
  * [12] G. Casella and R. Berger, Statistical Inference, Duxbury advanced series, Duxbury Thomson Learning, 2002, https://books.google.com/books?id=ZpkPPwAACAAJ.
  * [13] J. Chen and X. Huo, Theoretical results on sparse representations of multiple-measurement vectors, IEEE Transactions on Signal Processing, 54 (2006), pp. 4634–4643.
  * [14] V. Churchill, R. Archibald, and A. Gelb, Edge-adaptive $\ell_{2}$ regularization image reconstruction from non-uniform Fourier data, arXiv preprint arXiv:1811.08487, (2018).  https://arxiv.org/abs/1811.08487.
  * [15] V. Churchill and A. Gelb, Detecting edges from non-uniform Fourier data via sparse Bayesian learning, Journal of Scientific Computing, (2019), pp. 1–22.
  * [16] I. Daubechies, R. DeVore, M. Fornasier, and C. S. Güntürk, Iteratively reweighted least squares minimization for sparse recovery, Communications on Pure and Applied Mathematics, 63 (2010), pp. 1–38.
  * [17] A. P. Dempster, N. M. Laird, and D. B. Rubin, Maximum likelihood from incomplete data via the em algorithm, Journal of the royal statistical society. Series B (methodological), (1977), pp. 1–38.
  * [18] D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory, 52 (2006), pp. 1289–1306.
  * [19] J. Ellis, A practical guide to mcmc part 1: Mcmc basics, 2018, https://jellis18.github.io/post/2018-01-02-mcmc-part1/ (accessed 2020-04-29).
  * [20] D. Gamerman and H. F. Lopes, Morkov chain monte carlo, 68 (2006).  MCMC.
  * [21] A. Gelb and T. Scarnati, Reducing effects of bad data using variance based joint sparsity recovery, Journal of Scientific Computing, 78 (2019), pp. 94–120.
  * [22] A. Gelb and E. Tadmor, Detection of edges in spectral data, Applied and Computational Harmonic Analysis, 7 (1999), pp. 101–135.
  * [23] A. Gelb and E. Tadmor, Adaptive edge detectors for piecewise smooth data based on the minmod limiter, J. Sci. Comput., 28 (2006), pp. 279–306.
  * [24] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57 (1970), pp. 97–109, https://doi.org/10.1093/biomet/57.1.97, https://doi.org/10.1093/biomet/57.1.97, https://arxiv.org/abs/https://academic.oup.com/biomet/article-pdf/57/1/97/23940249/57-1-97.pdf.
  * [25] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Applied Mathematical Sciences, Springer New York, 2006, https://books.google.com/books?id=h0i-Gi4rCZIC.
  * [26] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, The Journal of Chemical Physics, 21 (1953), pp. 1087–1092, https://doi.org/10.1063/1.1699114, https://doi.org/10.1063/1.1699114, https://arxiv.org/abs/https://doi.org/10.1063/1.1699114.
  * [27] K. Murphy, Machine Learning: A Probabilistic Perspective, Adaptive Computation and Machine Learning series, MIT Press, 2012, https://books.google.com/books?id=NZP6AQAAQBAJ.
  * [28] J. Piironen, A. Vehtari, et al., Sparsity information and regularization in the horseshoe and other shrinkage priors, Electronic Journal of Statistics, 11 (2017), pp. 5018–5051.
  * [29] C. Price, Autocorrelation function (ACF).  https://www.mathworks.com/matlabcentral/fileexchange/30540-autocorrelation-function-acf, 2020\.  accessed: October 2020.
  * [30] C. Robert and G. Casella, Introducing Monte Carlo Methods with R, Use R!, Springer, 2010, https://books.google.com/books?id=WIjMyiEiHCsC.
  * [31] T. Sanders, A. Gelb, and R. Platte, Composite sar imaging using sequential joint sparsity, Journal of Computational Physics, 338 (2017), https://doi.org/10.1016/j.jcp.2017.02.071.
  * [32] 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.
  * [33] T. Scarnati and A. Gelb, Accurate and efficient image reconstruction from multiple measurements of fourier samples, Journal of Computational Mathematics, 88 (2020), pp. 798–828.
  * [34] M. E. Tipping, Sparse bayesian learning and the relevance vector machine, Journal of machine learning research, 1 (2001), pp. 211–244.
  * [35] D. P. Wipf and B. D. Rao, Sparse bayesian learning for basis selection, IEEE Transactions on Signal processing, 52 (2004), pp. 2153–2164.