###### Abstract

Reconstructing images from downsampled and noisy measurements, such as MRI and
low dose Computed Tomography (CT), is a mathematically ill-posed inverse
problem. We propose an easy-to-use reconstruction method based on Expectation
Propagation (EP) techniques. We incorporate the Monte Carlo (MC) method,
Markov Chain Monte Carlo (MCMC), and Alternating Direction Method of
Multiplier (ADMM) algorithm into EP method to address the intractability issue
encountered in EP. We demonstrate the approach on complex Bayesian models for
image reconstruction. Our technique is applied to images from Gamma-camera
scans. We compare EPMC, EP-MCMC, EP-ADMM methods with MCMC only. The metrics
are the better image reconstruction, speed, and parameters estimation.
Experiments with Gamma-camera imaging in real and simulated data show that our
proposed method is convincingly less computationally expensive than MCMC and
produces relatively a better image reconstruction.

Keyword: Alternating Directional Method of Multiplier, Gamma-camera imaging,
Markov Chain Monte Carlo, Expectation Propagation.

Image Reconstruction by Splitting Expectation Propagation Techniques from
Iterative Inversion

Robert G. Aykroyd

Department of Statistics, School of Mathematics,

University of Leeds, UK

<EMAIL_ADDRESS>
Kehinde Olobatuyi

Department of Statistics and Mathematical Finance,

University of Milan-Bicocca, Italy

<EMAIL_ADDRESS>
##  1 Introduction

The Bayesian approach to modeling data provides a flexible mathematical
frameworks for incorporating uncertainty of unknown quantities within complex
statistical models. The Bayesian posterior distribution encodes the
probabilistic uncertainty in the model parameters and can be used for
predictions of new unobserved data. In general, the posterior distribution is
intractable and it is therefore necessary to approximate it. Deterministic
approximations, such as Laplace approximation (see Bishpop 2006), variational
Bayes (Blei, Kucukelbir, and McAuliffe 2017), and expectation-propagation
(Minka 2001) aim to approximate the posterior with a simpler tractable
distribution (e.g., a normal distribution). The deterministic approaches are
often fit using fast optimization techniques and trade-off exact posterior
inference for computational efficiency.

Markov chain Monte Carlo (MCMC) algorithms (Brooks et al. 2011) approximate
the posterior distributions with a discrete set of samples generated from a
Markov chain whose invariant distribution is the posterior distribution.
Simple MCMC algorithms, such as random-walk Metropolis (Metropolis et al.
1953), are easy to apply and only require that the unnormalized density of the
posterior can be evaluated point-wise. Under mild conditions, the samples
generated from the Markov chain converge to the posterior distribution
(Roberts and Rosenthal 2004) and for many popular MCMC algorithms, rates of
convergence based on geometric ergodicity have been established (Meyn and
Tweedie 1994; Roberts and Rosenthal 1997)

While MCMC algorithms have the advantage of providing asymptotic exact
posterior samples, this comes at the cost of beign computationally slow to
apply in practice. This issue is further exacerbated by the demand to store in
the memory and analyze large-scale datasets and to fit increasingly
sophisticated and complex models to the high-dimensional data. For example,
brain imaging (Anderson et al. 2018), natural language processing (Yogatama et
al. 2014), popular genetics (Raj, Stephens, and Pritchard 2014), commonly use
a Bayesian approach to data analysis, but the continual growth in the size of
the dataset in these fields prevents the use of traditional MCMC methods.
Computational challenges such as these have led to recent research interest in
scalable Monte Carlo algorithms. Broadly speaking, the new developments in
Monte Carlo algorithm achieve computational efficiency by either parallelizing
the MCMC scheme, or by subsampling the data.

If the data can be distributed across multiple computer cores then the
computational challenge of inference can be shared, with MCMC algorithm run on
each core to draw samples from a partial posterior that is conditional on the
subset of the full data. The major challenge is to combine these posterior
samples from each computer to generate an approximation to the full posterior
distribution. It is then possible to merge the sample that are exact if the
partial posteriors are Gaussian (Scott et al. 2016) with the update rules that
just depend on the mean and variance for each partial posterior. However, it
is difficult to quantify the level of approximation of the rules due to non-
Gaussianity of the partial posterior. Alternative merging procedures that are
more robust to non-Gaussianity, have been proposed (Neiswanger, Wang, and Xing
2014; Rabionovich, Angelino, and Jordan 2015; Minsker et al. 2017; Nemeth and
Sherlock 2018; Srivastava, Li, and Dunson 2018) but it is difficult to measure
the level of approximation accuracy. Bespoke methods are also needed when
interest in the joint posterior of the parameters relates to subsets of the
data, or individual data points (Zuanetti et al. 2019).

We briefly review the standard EP algorithm as it forms the basis for the
novel proposed methods presented in this paper.

##  2 Expectation Propagation

Suppose a dataset $\mathcal{D}$ is observed which consists of $n$ independent
and identically distributed components such that
$\mathcal{D}=\\{\mathbf{x}_{i}\\}_{i=1}^{n}$. Further, each component is from
the parametric probabilistic model $p(\mathbf{x}|\bm{\theta})$ parameterized
by an unknown $p$-dimensional vector $\bm{\theta}$ whose prior distribution is
$p_{0}(\bm{\theta})$. Given the data, the posterior distribution is formulated
as follows;

$p(\bm{\theta}|\mathcal{D})=\frac{p_{0}(\bm{\theta})}{p(\mathcal{D})}\;\,\prod_{i=1}^{n}p(\mathbf{x}_{i}|\bm{\theta}).$
(1)

Exact Bayesian inference then involves the computation of typically
intractable posterior distribution over the parameters. Note that in what
follows, to simplify notation, the posterior distribution is denoted
$p(\bm{\theta})$ rather than the more usual $p(\bm{\theta}|\mathcal{D})$.

The goal of EP is to capture the contribution of each data point, through the
likelihood term, to the posterior by refining approximate terms, that is;

$p(\bm{\theta})\propto
p_{0}(\bm{\theta})\prod_{i=1}^{n}t_{i}(\bm{\theta})\,\approx\,q(\bm{\theta})\propto
p_{0}(\bm{\theta})\prod_{i=1}^{n}\tilde{t}_{i}(\bm{\theta}),$ (2)

where $t_{i}(\bm{\theta})\propto p(\mathbf{x}_{i}|\bm{\theta})$ is the exact
term, $\tilde{t}_{i}(\bm{\theta})$ is the corresponding approximating term for
$t_{i}(\bm{\theta})$, and $q(\bm{\theta})$ is the approximate posterior
distribution that will be refined by EP. One approach is to search for the
approximate posterior distribution which minimizes the Kullback-Leibler (KL)
divergence between the exact posterior distribution and the distribution
formed by replacing, in turn, one of the likelihood components by the
approximate term $\tilde{t}_{i}(\bm{\theta})$, that is
$\text{KL}[p(\bm{\theta}|\mathcal{D})~{}||~{}p(\bm{\theta}|\mathcal{D})\tilde{t}_{i}(\bm{\theta})/p(\mathbf{x}_{i}|\bm{\theta})]$;
Even with this approach, however, the update remains intractable in high
dimensional space. Instead, EP approximates this procedure by removing the
observations one-by-one from the posterior distribution, that is
$p_{-i}(\bm{\theta})\propto p(\bm{\theta}|\mathcal{D})/t_{i}(\bm{\theta})$,
which is called the exact leave-one-out posterior, and the corresponding
approximate leave-one-out posterior, that is $q_{-i}(\bm{\theta})\propto
q(\bm{\theta})/\tilde{t}_{i}(\bm{\theta})$, which is called the cavity
distribution. Since this couples the updates for the approximating factors,
the updates must now be iterated until convergence.

In a more general perspective, EP first selects a component to update and then
removes it from the approximate posterior distribution to produce the cavity
distribution, and the corresponding exact likelihood term is combined with the
cavity distribution to produce the tilted posterior distribution,
$\tilde{p}_{i}(\bm{\theta})\propto q_{-i}(\bm{\theta})t_{i}(\bm{\theta})$.
Then, EP updates the approximate posterior distribution by minimizing the
tilted distribution, that is;

$\tilde{t}_{i}(\bm{\theta})=\operatorname*{argmin}_{q(\bm{\theta})}\text{KL}[\tilde{p}(\bm{\theta})\;||\;q(\bm{\theta})],\quad
i=1,\ldots,n.$ (3)

This establishes that the measure of contribution the exact likelihood term
makes on both exact posterior distribution and the tilted posterior
distribution is the same. The KL minimization often reduces to moment matching
when the approximate distribution is in the exponential family. Finally, we
update each approximate likelihood term by dividing the new approximate
posterior distribution by the old approximate posterior distribution.

In practice, EP often performs well due to this local updating but it can
become intractable again in high dimensional problems due to the need to
compute normalizing factors caused by marginalizing the parameters.

Many attempts have been made to solve instability problem in Expectation
Propagation (EP) technique [(Vehtari et al.,, 2019); (Gelman et al.,, 2014)].
Several research have focused on the problem of the energy function in
expectation propagation. EP may have stability issues even with one-
dimensional tilted distributions and moment computations are more challenging
if the tilted distribution is multimodal or has long tails (Jylanki et al.,,
2011). Matthias and Bangti, (2014) studies expectation propagation method and
demonstrates its significant potentials on the nonlinear inverse problem. They
propose the coupling of the EP with an iterative linearization strategy. A
direct integral was used for computing the marginal likelihood. However, this
is only tractable in a low dimensional space as the work assumed. Jose and
Ryan, (2015) works on a hierarchical problem as a layered Bayesian Neural
Network. However, their approach can lead to underestimation of variance of
the approximate posterior. Gelman et al., (2014) revisit expectation
propagation as a prototype for scalable algorithm that partition big datasets
into many parts and analyze each part in parallel to perform inference of
shared parameters. Minka, (2004) proposes an extension of EP called a
fractional EP to improve the robustness of EP algorithm when the approximation
family is not enough (Minka,, 2005) or when the propagation of information is
difficult due to vague prior information (Seeger,, 2008). Yingzhen et al.,
(2018) presents an extension of expectation propagation called a Stochastic
Expectation Propagation that maintains a global posterior approximation like
variational inference but updates the posterior in a local way. The work used
the expectation propagation strategy but differs only at the updating stage.
Graves, (2011) proposes a Monte Carlo approximation to compute the lower bound
of the variational inference optimized by stochastic gradient descent (SDG).
John et al., (2011) focuses on the problem of intractable lower bound on the
marginal likelihood in variational inference at the updates of the parameters.
They propose an alternative approach based on stochastic optimization that
allows for direct optimization of the variational lower bound.

Image reconstruction problems are common to medical imaging area such as fast
MRI and low dose CT are generally mathematically ill-posed inverse problems.
Often times, linear imaging system are considered with a forward operator $G$,
e.g. a Fourier transform for MRI and X-ray transform for CT. The measurement
$y$ is given as $y=Gx$, where $x$ is the underlying image in the noise-free
state. The linear operator $G$ is ill-posed for most applications; therefore,
some statistical priors are necessary to make these problems invertible.
Consequently, intractability becomes inevitable due to high-dimensional data.
Moreover, when the number of hidden variables to be estimated grows bigger and
results in multiple integration, intractability poses a problem in Bayesian
inference. In this paper, we generalize EP to the mixture of Gaussian and non-
Gaussian distributions. We propose a new hybrid algorithms called Splitting
Expectation Propagation (SEP). We incorporate the Monte Carlo approximation,
Markov Chain Monte Carlo method (MCMC), and Alternating Direction Method of
Multiplier (ADMM) into the EP. SEP has the ability to handle non-Gaussian
factors in hierarchical Bayesian models in high dimensional space.

The paper is organized as follows; In Sect. (2), we give a detailed
explanation of the standard EP technique whereas the proposed techniques is
fully derived in Sect. (3). Furthermore, to illustrate the methods, in Sect.
(4), we present an experiment on synthetic data sets and finally, in Sect.
(5), we draw some future direction and the main conclusions of this work.

##  3 Splitting Expectation Propagation

In this section, we introduce a new hybrid algorithm motivated by EP, called
Splitting Expectation Propagation (SEP). SEP incorporates Monte Carlo (MC),
Markov Chain Monte Carlo (MCMC), and Alternating Direction method of
multiplier (ADMM) at the EP update. The algorithm can be interpreted as a
version of EP that is capable of handling the marginal likelihood with mixture
of exponential family in a stochastic and deterministic views.

###  3.1 Model Specification

We describe a Bayesian model for the data based on an inverse problem. Given
data $\mathcal{D}=\\{\mathcal{X}_{n},\mathcal{Y}_{n}\\}_{n=1}^{N}$, made up of
$D-$dimensional vector $\mathcal{X}_{n}\in\mathcal{R}^{D}$ and corresponding
target variables $\mathcal{Y}_{n}\in\mathcal{R}^{D}$. In the context of
inverse problems, we only have access to a noisy version $\mathcal{Y}$ of the
exact data $\mathcal{X}$. We assume that $\mathcal{Y}$ is obtained as
$\mathcal{Y}=G\mathcal{X}+\xi$, where the
$\xi\in\mathcal{R}^{\mathrm{N}\times\mathrm{D}}$ represents the noise in the
data.

$\lambda$$\alpha_{\lambda}$$\mathcal{Y}$$\mathcal{X}$$\tau$$\mathbf{\alpha_{\tau}}$
Figure 1: Graph networks connecting priors and their hyper-priors

Figure (1) shows graphical representation of the Bayesian hierarchical model
which is mathematically represented as follows;

$p(\mathcal{X},\lambda,\tau|\mathcal{Y})=\frac{p(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda)p(\mathcal{X}|\tau)p(\lambda)p(\tau)}{\int...\int
p(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda)p(\mathcal{X}|\tau)p(\lambda)p(\tau)d\mathcal{X}d\lambda
d\tau},$ (4)

Equation (4) has a large size of the integration which grows with the size of
the reconstructed image $\mathcal{X}$. The standalone or independent variable
is directly connected to both its mean $\mathcal{X}$ and variance $\lambda$.
The priors are disconnected from each other but directly connected to their
respective hyper-priors. We consider a linear equation

$\mathcal{F}(\mathcal{X})=\mathcal{Y}$ (5)

where the map
$\mathcal{F}:\mathcal{R}^{\mathrm{N}\times\mathrm{D}}\rightarrow\mathcal{R}^{\mathrm{N}\times\mathrm{D}}$,
matrices $\mathcal{X}\in\mathcal{R}^{\mathrm{N}\times\mathrm{D}}$ and
$\mathcal{Y}\in\mathcal{R}^{\mathrm{N}\times\mathrm{D}}$ refer to data
formation mechanism, unknown parameter and the given data respectively.

Let $\mathcal{Y}$ be an $N\times D$ matrix and $\mathcal{X}$ be an $N\times D$
matrix also. The likelihood for the model and the noise variance $\lambda$,
with data $\mathcal{D}=(\mathcal{X},\mathcal{Y})$ is then

$p(\mathcal{Y}|G\mathcal{X},\lambda)=\prod_{i}\prod_{j}\mathcal{N}(\mathcal{Y}_{ij}\big{|}\mathcal{G}\mathcal{X}_{ij},\lambda),$
(6)

Also, we specify a Gaussian prior distribution for each entry
$\mathcal{L}\mathcal{X}$ which is used as a matrix with zero mean and
$\tau^{2}$ as its variance. The Gaussian prior for $\mathcal{L}\mathcal{X}$ is
as follows;

$p(\mathcal{L}\mathcal{X}|\tau)=\prod_{i}\prod_{j}\mathcal{N}(\mathcal{L}\mathcal{X}_{ij}\big{|}0,\tau),$
(7)

The prior distributions for $\lambda$ and $\tau$ are chosen to be an
exponential distributions, i.e., $p(\tau)=\exp(\tau|\alpha_{\tau})$ and also
the noise variance to be $p(\lambda)=\exp(\lambda|\alpha_{\lambda})$. The
posterior distribution for the parameters $\mathcal{X},\lambda,\tau$ can then
be obtained by

$p(\mathcal{X},\lambda,\tau|\mathcal{Y})=\frac{p(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda)p(\mathcal{L}\mathcal{X}|\tau)p(\lambda)p(\tau)}{p(\mathcal{Y})},$
(8)

where $\lambda$, $\tau^{2}$ are unknown variances and $\mathcal{G}$ is a known
matrix, $p(\mathcal{Y})$ is the normalizing factor and it will be denoted as
$Z$ throughout this work and in particular

$Z=\int....\int
p(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda)p(\mathcal{L}\mathcal{X}|\tau)p(\lambda)p(\tau)d\mathcal{X}d\lambda
d\tau,$ (9)

where $\mathcal{X}$ is $\mathcal{N}\times D$ dimension and $\mathcal{Y}$ is
$\mathcal{N}\times D$ dimension. Next, we choose an approximating families.
Here, the approximate distributions are mixture of different distributions of
the exponential family. For instance,

$q(\mathcal{X},\lambda,\tau)=\bigg{[}\prod_{i=1}^{m}\prod_{j=1}^{n}\mathcal{N}(\mathcal{X}_{ij}\big{|}mx_{i,j},vx_{i,j})\bigg{]}\exp(\lambda|\alpha_{\lambda})\exp(\tau|\alpha_{\tau}),$
(10)

$\mathcal{X}\sim\mathcal{N}(mx_{i,j},vx_{i,j})$,
$\lambda\sim\exp(\alpha_{\lambda})$ and $\tau\sim\exp(\alpha_{\tau})$.

####  3.1.1 EP-MC Technique

Here, the approximation parameters $mx_{i,j},vx_{i,j},\alpha_{\tau}$ and
$\alpha_{\lambda}$ are determined by applying a stochastic search on the
posterior in Equation (8). Finally, we sequence through and incorporate the
terms $t_{i}$ into the approximate posterior in Equation (10). At each step,
we move from an old $q_{-i}(\mathcal{X},\lambda,\tau)$ to a new
$q(\mathcal{X},\lambda,\tau)$.

The update rule is given in Equation (11) below;

$\alpha_{\tau}^{new}\approx\alpha_{-ij}+v_{\alpha_{-ij}}\zeta_{\alpha_{-ij}},$
(11)

Where the mean $\frac{1}{K}\sum_{k=1}^{K}\tau^{(k)}$ is denoted as
$\zeta_{\alpha_{-ij}}$ and variance as $v_{\alpha_{-ij}}$.

The update rule for $\lambda$ is given in Equation (12) below;

$\alpha_{\lambda}^{new}\approx\alpha^{\lambda}_{-ij}+v_{\alpha^{\lambda}_{-ij}}\frac{1}{K}\sum_{k=1}^{K}\lambda^{(k)},$
(12)

We Take the derivative of the normalizing factor with respect to $\omega$
concerning $\mathcal{X}$ as follows;

$\nabla_{\omega}\log{Z}=\frac{1}{Z}\int_{\mathcal{X}_{i,j},\tau}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}q_{i,j}(\mathcal{X}_{i,j},\tau)\nabla_{\omega}\log{q_{i,j}(\mathcal{X}_{i,j},\tau)}d\mathcal{X}_{i,j}d\tau,$
(13)

At this juncture, we proceed with the use of Monte Carlo Integration

$\frac{\delta_{1}}{\delta_{0}}=\sum_{n=1}^{N}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}\nabla_{\omega}\log{q_{i,j}(\mathcal{X}_{i,j},\tau)}\bigg{/}\sum_{n=1}^{N}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)},$
(14)

when taking $\mathcal{X}$ into consideration then we fix $\tau$ and vice
versa. So we update the parameter by taking the gradient step

$\omega^{t+1}=\omega^{t}+lr*\frac{\delta_{1}}{\delta_{0}},$ (15)

One of the peculiarities of EP is the ability to minimize Kullback-Leibler
divergence between the tilted posterior and the approximate posterior. The
full derivatives are provided in the Appendix.

####  3.1.2 EP-ADMM technique

We introduce the Alternating direction method of multiplier (ADMM) algorithm
to update the approximate posterior parameters.

$\displaystyle\begin{split}\text{minimize}\hskip
14.45377pt\text{KL}(\tilde{p}(x_{ij}||q(x)_{ij}))\end{split}$ (16)
$\displaystyle\ \text{subject to}:m_{ij}\geq a;v_{ij}\geq b$

where $a$ and $b$ are constants. Then updating according to Minka, (2001) is
as follows;

$m_{x}=\nabla_{m_{-ij}}\log Z_{x}+\alpha+\rho(m_{-ij}-a)$ (17)

and

$v_{x}=\nabla_{v_{-ij}}\log Z_{x}+\beta+\rho(v_{-ij}-b)$ (18)

where

$q(x)=\mathcal{N}(\mathcal{X};m_{x},v_{x})$ (19)

then the tilted distribution is

$\tilde{p}(x_{ij})=\frac{t_{ij}(x)q_{-ij}(x)}{\int_{\mathcal{X}_{ij}}t_{ij}(x)q_{-ij}(x)dx_{ij}}$
(20)

Now the normalizing factor is

$Z_{x}=\int_{x_{ij}}\mathcal{N}(\mathcal{Y}_{ij};G\mathcal{X}_{ij},\lambda)\mathcal{N}(\mathcal{X}_{ij};m_{-ij},v_{-ij})dx_{ij}$
(21)

The update rule for the parameters $m_{x}$ and $v_{x}$ of $q(\mathcal{X})$
according to Minka, (2001) with a method of multipliers are

$m^{new}_{x}=m_{-ij}+v_{-ij}\frac{y_{ij}-gm_{-ij}}{v_{-ij}g^{2}+\lambda}g+\alpha+\rho(m_{i}-a)$
(22)
$v_{x}^{new}=\frac{v_{-ij}\lambda}{v_{-ij}g^{2}+\lambda}+\beta+\rho(v_{i}-b)$
(23) $\alpha^{new}=\alpha^{k}+\rho(m_{x}^{new}-a)\hskip
14.45377pt\text{and}\hskip
14.45377pt\beta^{new}=\beta^{k}+\rho(v_{x}^{new}-b)$ (24)

for full derivatives, Appendix provides details.

####  3.1.3 EP-MCMC Technique

The EP-MCMC method provides an approach when the model is complex or the
number of parameters is large and deterministic gradient-based techniques are
infeasible and most importantly when the posteriors are not Gaussian. The
Kullback-Leibler divergence between the tilted posterior and the approximated
posterior can be addressed using the MCMC technique. In brief, the transitions
in the Markov chain are designed so that an equilibrium distribution exists
and is equal to the target distribution (tilted posterior). According to EP’s
strategy of moment matching, this is also possible but instead we update the
moments of the approximated posterior with the moment of the tilted posterior
from the Markov chain. We update the approximate term and compute the new
cavity distribution from the new approximated posterior and thereafter compute
the new tilted posterior. With this cycle in mind, the advantage of EP-MCMC
over the MCMC is that at each iteration in EP-MCMC the tilted posterior is
updated with new cavity distribution thereby generating samples that are
independent at each step of the algorithm. Unlike EP-MCMC, the major problem
of MCMC is the dependence among samples generated at different iterations.
This is due to a static target posterior. The EP setup for the parameters
$\tau$ and $\lambda$ is approached via Markov Chain Monte Carlo method.

Algorithm 1 The General EP-MCMC Algorithm

Set an initial value for $\tau$ and $\lambda$

  1. 1.

Initialize all of the approximating factors $\tilde{t_{i}}(\tau)$ and
$\tilde{t_{i}}(\lambda)$.

  2. 2.

Compute the initial approximation $q(\tau)$ and $q(\lambda)$ from the product
of the approximating factors:

$q(\tau)=\frac{\prod_{i}\tilde{t}_{i}}{\int\prod_{i}\tilde{t}_{i}\hskip
3.61371ptd\tau}\hskip 14.45377pt\text{and}\hskip
14.45377ptq(\lambda)=\frac{\prod_{i}\tilde{t}_{i}}{\int\prod_{i}\tilde{t}_{i}\hskip
3.61371ptd\lambda}$

  3. 3.

Until all $\tilde{t}_{i}$ converge:

    1. (a)

Choose a $\tilde{t}_{i}$ to refine

    2. (b)

Remove $\tilde{t}_{i}$ from the approximation $q(\tau)$ and $q(\lambda)$ by
division:

$q_{-i}(\tau)\propto\frac{q(\tau)}{\tilde{t}_{i}(\tau)}\hskip
14.45377pt\text{and}\hskip
14.45377ptq_{-i}(\lambda)\propto\frac{q(\lambda)}{\tilde{t}_{i}(\lambda)}$

    3. (c)

Compute the tilted distribution $\tilde{p}_{i}$ from $q_{-i}$ and the exact
factor $t_{i}$

$\tilde{p}_{i}(\tau)\propto t_{i}(\tau)\hskip 2.84544ptq_{-i}(\tau)\hskip
14.45377pt\text{and}\hskip 14.45377pt\tilde{p}_{i}(\lambda)\propto
t_{i}(\lambda)\hskip 2.84544ptq_{-i}(\lambda)$

      1. i.

Generate $\epsilon_{\tau}$ and $\epsilon_{\lambda}$ from a Gaussian
distribution

      2. ii.

Generate a proposed new value $\tau^{*}=\tau+\epsilon_{\tau}$ and
$\lambda^{*}=\lambda+\epsilon_{\lambda}$

      3. iii.

Evaluate

$\alpha_{\tau}=\frac{\pi(\tau^{*}|\mathcal{X})}{\pi(\tau|\mathcal{X})}\hskip
14.45377pt\text{and}\hskip
14.45377pt\alpha_{\lambda}=\frac{\pi(\lambda^{*}|\mathcal{Y})}{\pi(\lambda|\mathcal{Y})}$

      4. iv.

Generate $u_{\tau}$ and $u_{\lambda}$ from a uniform distribution $U(0,1)$

      5. v.

If $\alpha_{\tau}>u_{\tau}$ and $\alpha_{\lambda}>u_{\lambda}$ then accept the
proposals $\tau^{*}$ and $\lambda^{*}$ and set $\tau=\tau^{*}$ and
$\lambda=\lambda^{*}$

    4. (d)

Update the approximate posterior from the inference of tilted posterior, this
is similar to a moment matching in EP. $q^{new}(\tau)$ and $q^{new}(\lambda)$

    5. (e)

Compute the new approximate term:

$\tilde{t}^{new}_{i}(\tau)\propto\frac{q^{new}(\tau)}{q_{-i}(\tau)}\hskip
14.45377pt\text{and}\hskip
14.45377pt\tilde{t}^{new}_{i}(\lambda)\propto\frac{q^{new}(\lambda)}{q_{-i}(\lambda)}$

  4. 4.

Evaluate the approximation to the model evidence:

$p(\mathcal{\mathcal{Y}})\approxeq\int\prod_{i}\tilde{t}_{i}(\tau)\hskip
3.61371ptd\tau\hskip 14.45377pt\text{and}\hskip
14.45377ptp(\mathcal{\mathcal{Y}})\approxeq\int\prod_{i}\tilde{t}_{i}(\lambda)\hskip
3.61371ptd\lambda$ (25)

##  4 Implementation details

After incorporating the factors in Equation (20) for the first time, we
constrained only mean and variance parameters of approximate posterior $q$ of
$\mathcal{X}$ for the purpose of this work using Alternating direction method
of multipliers. Because of the structure of the data presented in Figure (4)
which is sparse, EP often breaks down when the algorithm produces variance
parameter that is negative after incorporating one likelihood factor,
especially at the update step of approximate term which involves division. A
similar operation was reported in Minka, (2001) when negative variances arise
in Gaussian approximation factors.

###  4.1 Monitoring Convergence of EP-MCMC

To monitor convergence of the embedded Metropolis Hastings in EP, we run EP-
MCMC $m=10$ replications and $n=1000$ to compute the Brook-Gelman statistic;
Brooks and Gelman, (1998) which compare the between-chain variance and within-
chain variance. However, EP is generally known to be deterministic which is
still preserved in EP-ADMM for the prior on $\mathcal{X}$ but the priors on
$\lambda$ and $\tau$ are randomly tuned. The reconstruction images presented
are mean of the approximate posterior $m_{ij}$. To know if a convergence has
been reached, Gelman and Rubin, (1992) suggests comparing $m$ inferences
computed from the $m$ chains to the inferences computed by mixing together the
$mn$ draws from all the sequences. To compute the between-chain and within-
chain variances, let $s_{jt}$ where $j=1,...,m$ and $t=1,...,n$ be the $t$th
iterations of $s$ in sequence $j$, the between-chain variance $B/n$ and
within-chain variance are given by

$B/n=\frac{1}{m-1}\sum_{j=1}^{m}(\bar{s}_{j.}-\bar{s}_{..})^{2}$ (26)

and

$W=\frac{1}{m(n-1)}\sum_{j=1}^{m}\sum_{t=1}^{n}(s_{jt}-\bar{s}_{j.})^{2}$ (27)

We also computed the potential scale reduction factor PSRF which is the
variance of the pooled and within-chain inferences

$R=\frac{\hat{V}}{\sigma^{2}}$

R will be estimated by $\hat{R}$ because its denominator is unknown

$\hat{R}=\frac{\hat{V}}{W}=\frac{m+1}{m}\frac{\hat{\sigma}^{2}_{+}}{W}-\frac{n-1}{mn}$
(28)

where the pooled variance is given by

$\hat{V}=\hat{\sigma}^{2}_{+}+B/(mn)$ (29)

and the estimated variance of the $\sigma^{2}$ given by the weighted average
of $B$ and $W$

$\hat{\sigma}^{2}_{+}=\frac{(n-1)}{n}W+\frac{B}{n}$ (30)

Equation (30) is the unbiased estimate of the true variance $\sigma^{2}$. And
$\hat{V}$ accounts for the sampling variability of the estimator which yields
a pooled posterior variance estimate given in Equation (30). The PSRF is
expected to be close to $1$ for convergence to be guaranteed and each of the
$m$ sets of $n$ simulated observations is close to the target distribution,
however large value of $\hat{R}$ indicates no convergence or divergence and
further simulations can be taken into consideration or the proposal
distribution should be critically looked into.

###  4.2 Splitting EP on Clutter Problem

The clutter problem has been addressed in the work by Minka, (2001). This work
replicates this process to compare the proposed algorithm with the EP
algorithm for a one-dimensional space. Suppose we have observations from a
Gaussian distribution embedded in a sea of unrelated clutter where $w$ is the
clutter ratio, so that the density observation is a mixture of two Gaussians:

$p(\textbf{y}|\theta)=(1-w)\mathcal{N}(\textbf{y};\theta,\textbf{I})+w\mathcal{N}(\textbf{y};\textbf{0},10\textbf{I})$
(31)

Let the d-dimensional vector $\vec{\theta}$ with a Gaussian prior
distribution:

$p(\theta)\sim\mathcal{N}(\textbf{0},100\textbf{I})$ (32)

The joint distribution of $\vec{\theta}$ and $n$ observation
$D=\\{\bf{y}_{1},...,\bf{y}_{n}\\}$

$p(D,\theta)=p(\theta)\prod_{i}p(\textbf{y}_{i}|\theta)$ (33)

Figure 2: A complex posterior in the clutter problem produced by EP with ADMM
and EP compared with the exact posterior.

Figure 3: A complex posterior in the clutter problem produced by Stochastic
search EP with Monte Carlo integration compared with EP.

The goal is to approximate the posterior distribution

$p(\vec{\theta}|D)=\frac{p(D,\theta)}{\int p(D,\theta)\hskip 1.42271ptd\hskip
1.42271pt\theta}$ (34)

The first component in Equation 31 contains the parameter to be estimated
which is $\theta$ and the second component describes the clutter where $w$ is
the known ratio of the clutter. The Bayesian network for this problem is
$\theta$ point to the $y_{i}$ [Check Minka, (2001) for more detail].

The strength of Splitting EP over expectation propagation is the use of
stochastic approximation and method of multiplier within EP. The learning rate
(lr) determines how stable the Splitting EP algorithm is, unlike expectation
propagation which fails at the refining step. Learning rate controls the
erratic nature of the expectation propagation algorithm, and this seems to be
a viable improvement over EP. Moreover, it forces the variance not to reduce
to negative value which often leads to EP’s stability issue. EP was compared
to four algorithms for approximate inference in Minka, (2001) such as the
Laplace’s method, variational Bayes, Importance sampling where the prior is
used as the importance sampling, and Gibbs sampling by introducing the hidden
variables that determines if a data point is clutter or not. According to
Minka, (2001), EP competed well with other deterministic algorithms by
approximating the posterior with a Gaussian. However, their performance
improved substantially with more data i.e. the posterior is more Gaussian with
more data. Figure (3) shows true posterior with multi-modal shape (Black
line). Here, EP-ADMM (Red line) captures at least two modes while EP captures
only one of the modes. This justifies that EP is a unimodal algorithm. Figure
(3) shows how Stochastic search expectation propagation competes well with EP.

###  4.3 Splitting EP on Synthetic Data

A cylindrical image produced by a gamma camera is presented in Figure (4). The
true image to be reconstructed is presented in Figure (4)(c) mixed with a
noise. We evaluate EP-Monte Carlo in inverse problem with hierarchical
Bayesian models, with data of cylindrical image from gamma camera. In EP-Monte
Carlo, we retain the originality of the EP. However, we improve on the
intractability of the normalizing factor by following the idea of John et al.,
(2011). Due to the randomness introduced to EP, we ran EP-MC $1000$ iterations
and output the mean of the posterior. The true parameters are the following;
the precision is $100$ and the true standard deviation is $0.1$. We hereby
present the reconstruction results of the EP-MC and MCMC. Figure (6) shows the
final output of EP-MC after $1000$ iterations while Figure (6) shows the
outcome of the MCMC. The reconstruction result of EP-Monte Carlo outperforms
that of MCMC in terms of sharpness of the image.

Figure 4: $(a)$ True data: $\mathcal{X}$, $(b)$ the mean:
$\mathcal{G}\mathcal{X}$ and $(c)$ Noisy data:$\mathcal{Y}$

Figure 5: The reconstruction of the cylindrical image produced by EP-MC

Figure 6: (a): The reconstruction image of the cylindrical image produced by
MCMC.(b): relative error (c): The estimates of precision converge at about
$85$, (d): The estimates of standard deviation converges at $0.085$.

Figure 7: The precision is calculated from the plot as $\tau^{-2}=100$. Since
the variance is $\tau^{2}=0.01$

Figure 8: The variance output from the EP-MC is $0.01$ then the estimated
standard deviation of EP-MC is $0.1$.

Table 1: The Accuracy of the Standard deviation and precision produced
by EP-MC and MCMC for prior $\tau$ and $\lambda$ compared to the True values provided Parameter |  | True-Value |  | EP-MC |  | MCMC |  |
---|---|---|---|---|---|---|---|---
$\tau^{-2}$ |  | $100$ |  | $100$ |  | $85$ |  |
$\lambda$ |  | $0.1$ |  | $0.1$ |  | $0.085$ |  |

Compared to the true image in Figure (4)(a), we observe that the posterior of
$\mathcal{X}$ is well approximated by the EP-MC method. In terms of the
convergent time, i.e. time for both methods to reach convergence,

Figure 9: (a): Mouse injected with $[^{99m}Tc]MDP$ at $4$ h pi (15-min)
acquisition time. (b): A mouse injected with ${}^{99m}Tc$ DMSA at $1h$ (10-min
scan), 3 h (10-min scan), 5 h (10-min scan), 6h (10-min scan), and 24 h pi
(30-min scan), (c): A static image of the mouse injected with the $[^{99m}Tc]$
MIBI at different intervals of time.

EP-MC is faster than MCMC. Unlike MCMC that needs about $500$ chains as a
burn-in period, EP-MC does not require any thin-in or burn-in period. We
observed a clear pattern in the residual plot shown in Figure (6)(b). This
non-randomness may be due to the choice of the starting value of
hyperparameters. The estimated precision for both EP-MC and MCMC are provided
in Figure (8) and Figure (6)(c) respectively. It is not surprising to observe
that EP-MC produces the exact precision parameter value because of the high
accuracy in the approximation of $\mathcal{X}$. Similarly, the estimate of
variance parameter value are provided in Figure (8) and Figure (6)(d) for EP-
MC and MCMC respectively. Table (1) shows the comparison of the estimate of
precision and variance parameter values for EP-MC and MCMC with the true
parameters. We observe that EP-MC provides a closer estimate of both precision
and variance to the true parameters.

###  4.4 Splitting EP on Real Data

The full description of the study on the $\gamma$-eye system evaluated in a
proof-of-concept animal study using normal Webster Swiss Albino mice with
average weights of $25g$ can be found in Maria et al., (2016). The author
conducted the study on three different clinical radio-pharmaceuticals which
was radio-labeled with Tc-99m and injected the mice via their tail vein. The
first mouse presented in Figure (9)(a) was injected with $100\mu 1/7.5$ MBq
$[^{99m}Tc]$ MDP, which is a suitable agent for bone imaging and static images
were obtained at $1h$ post-injection $(pi)$, $2h$ pi, $3h$ pi, and $4h$ pi.
The second mouse was injected by $100\mu 1/7.5$ MBq $[^{99m}Tc]$ DMSA, which
is a tracer for imaging the anatomical structure and functional process of the
kidneys. The author performed dynamic imaging for the first hour after
injection and then static images for $10$ min were acquired every $1h$ up to
$24h$ pi. The image presented in Figure (9)(b) is a 1-h dynamic study of a
mouse injected with $100\hskip 0.56905pt\mu l$ of $[^{99m}Tc]$ DMSA radio-
tracer. Then, static images were acquired at different time intervals. Figure
(9)(c) shows a third mouse injected with $100\mu 1/5.6$ MBq $[^{99mTc}]$ MIBI
for heart perfusion. Dynamic imaging was performed for $2h$ pi and $10$ min
static images were obtained up to $5$-h pi.

####  4.4.1 Parameter Setup for Splitting EP

The Splitting EP algorithm was initialized with the following parameters;
$\tau=\lambda=0.01$ as a starting point for Metropolis Hasting algorithm
embedded in EP to solve the intractability of the normalizing factor. We fixed
the hyper-priors to $\tilde{\alpha}_{ij}^{\tau}=2$,
$\tilde{\alpha}_{ij}^{\lambda}=1$ for the approximate term distribution and
$\alpha^{\tau}_{ij}=\alpha^{\lambda}_{ij}=10$ for the approximate posterior.
According to the update of splitting EP, we computed the cavity distribution
parameters $\alpha_{-ij}^{\tau}$ and $\alpha_{-ij}^{\lambda}$ just once to
initialize the target distribution $p(\tau)$ and $p(\lambda)$. The outer loop
consists of the EP-MCMC, i.e. tilted posteriors of $\tau$ and $\lambda$ are
hereby assumed to be intractable and thereby approximated using MCMC
technique. We note that EP-MCMC has some advantages over ordinary MCMC in
terms of speed to convergence, and accuracy. Also, we observed that due to
different values of cavity distribution computed at different iterations of
the algorithm, EP-MCMC does not need either burn-in period or thinning to
achieve stationarity. To illustrate the accuracy of EP-ADMM for the image
reconstruction and EP-MCMC for parameter recovery, we present also the results
by a standard Metropolis Hasting MCMC with a sample of $1000$ and burn-in
period of $500$.

Figure 10: (a): 3h static_mibi_results from EP-ADMM, (b): relative error (c):
estimates of $\tau$ converges at $158.48$ from EP-MCMC, (d): estimates of
$\lambda$ converges at $0.028$ from EP-MCMC.

Figure 11: (a): 3h static_mibi_results from MCMC, (b)&(c): error & residual,
(d): the estimate of $\sigma_{x}$ converges at $159$, (e): the estimates of
$\sigma_{\epsilon}$ converges to about $0.9$.

Figure 12: The Autocorrelation plot of the mouse injected with Mibi reagent
for sigma; the EP-MCMC algorithm was runs $10$ time each of length $1000$. The
plot shows independence of chains at each replication.

Figure 13: For the $\tau$, the plot shows a strong dependence of chains at
first four replications but at the fifth and ninth there is no correlation and
the last replication shows a decline in correlation of the chains.

The random walk stepsize is chosen so that the acceptance rate is close to
$0.234$, which is widely known to be optimal (Gelman et al.,, 1996). It is
well known that the convergence of MCMC is very crucial.

####  4.4.2 Reconstruction result for Mibi

First, we present the reconstruction result for the mouse injected with
$100\mu 1/5.6$ MBq $[^{99mTc}]$ MIBI for heart perfusion. The estimate based
on the Gaussian prior is shown in Figure (11). It can be seen that a well
reconstruction is obtained with a substantive reduction in the noise. The
result produced from ordinary MCMC is presented in Figure (11)(a). The
reconstructions produced by both algorithms have a clear distinctions. First
off, the reconstruction by EP-ADMM has a white background and smooth edges
while the reconstruction produced by MCMC has a gray background and rough
edges around the image reconstructed. This is a very important distinction
because expectation propagation is well known for its fast computational
speed. In contrast, MCMC is known for its slow convergence and it requires
more computational time to reach a stationary distribution. Figure (11)(b)
shows the error comparing the estimates with the data $\mathcal{Y}$. Comparing
the MCMC results with the reconstruction produced in EP-MCMC, the parameter
estimates of $\sigma_{x}$ in MCMC is $159$ while $\tau$ in EP-MCMC is about
$158.48$, similarly the parameter estimates of $\sigma_{\epsilon}$ is about
$0.9$ while that of EP-MCMC is $0.028$. The residual error shown in Figure
(11)(b) comparing the estimates with the data shows a foggy pattern. It is
possible to produce similar random residual as shown in Figure (11)(b&c) by
adjusting the value of $\alpha_{\lambda}$ and $\alpha_{\tau}$ in the
hyperparameter $\lambda$ and $\tau$ respectively. This may be achieved by
tweaking the values of the hyperparameters $\lambda$ and $\tau$. However,
given our primary aim to produce a constant reconstruction, our focus is on
the appropriate values of hyperparameters that produce smoother
reconstruction.

Figure 14: Iterative PSRF Plot for $\lambda$ in Mibi image data (from $m=10$
parallel sequence and $n=1000$). The convergence starts at about $300$
iterations till the end of the iteration.

Figure 15: Iterative PSRF Plot for $\tau$ in Mibi image data (from $m=10$
parallel sequence and $n=1000$). The convergence starts at About $800$
iteration till the end of the iteration.

Figure 16: Marginal posterior distribution of the parameter $\lambda$ for Mibi
data. The tick horizontal bar signifies the posterior $95\%$ credible interval
and the maximum value at $0.028$ is the posterior mean.

Figure 17: Marginal posterior distribution of the parameter $\tau$ for Mibi
data. The tick horizontal bar signifies the posterior $95\%$ credible interval
and the maximum value at about $160$ is the posterior mean.

Table 2: The mean estimates of the ten EP-MCMC chains for priors $\tau$ and
$\lambda$
of Mibi data parameter | Posterior-Mean | $\hat{R}$ | $\hat{V}$ | $W$ | $B/n$ | $\sigma_{+}^{2}$ |
---|---|---|---|---|---|---|---
$\tau$ | $158.48$ | $1.06$ | $173.58$ | $98.53$ | $68.32$ | $166.75$ |
$\lambda$ | $0.028$ | $1.00$ | $2.98e^{-08}$ | $2.49e^{-08}$ | $2.48e^{-08}$ | $2.73e^{-08}$ |

These parameter estimates are very significant in producing the clearer and
smooth image. Given our primary aim to produce a clear and smoother image,
which means larger $\tau$ in EP-MCMC and $\sigma_{x}$ in MCMC would be
preferred and lower values of $\sigma_{\epsilon}$ in MCMC and $\lambda$ in EP-
MCMC would also be preferred. The autocorrelations of the estimates produced
by MCMC are nested in their respective plot and it can be seen that the
autocorrelations have been controlled by the well set-up of MCMC such as burn-
in of $500$ samples and thin-in of $1$-of-$10$ samples. However, there seems
to be a high correlation in Figure (11)(e) which might be reduced over time
with large iteration. The computed potential scale reduction factors $\hat{R}$
for $\tau$ is $1.06$, with upper C.I to be $1.12$ and for $\lambda$ is $1$
with upper C.I to be $1$ as shown in Table (2). Here, the potential scale
reduction factor for $\lambda$ is exactly $1$ which indicates strong
convergence across the $10$ runs. According to Gelman and Rubin, (1992), large
$\hat{R}$ can be interpreted as a need for more simulations to further reduce
the estimate of the variance $\hat{\sigma}^{2}$ or to increase the $W$. The
closeness to $1$ indicates that each $m$ sets of $n$ simulated samples is
close to the target distribution. However, potential scale reduction factor
for $\tau$ is a bit more than one which is in tandem with further plots below.
The simulation result agrees with the result from MCMC. Figure (13) shows the
correlation between the chains of $\lambda$ for each run. It can be seen that
there is no autocorrection in the chains for

Figure 18: (a): mdp_1h_mouse_results from EP-ADMM, (b): relative error (c):
estimates of $\tau$ which converges at $165.443.$, (d): estimates of $\lambda$
converges at $0.019$.

Figure 19: (a): mdp_1h_mouse_results from MCMC (b)&(c): error & residual, (d):
the estimate $\sigma_{x}$ converges at $160$, (e): the estimates
$\sigma_{\epsilon}$ converges to about $0.9$.

each of the runs, this indicates that there is a convergence at every run
which also contributes to the potential scale reduction factor of $1$.
Similarly, Figure (13) shows an autocorrelation plot of $\tau$ for the $10$
runs. There is a high correlation among the chain of the first four runs, but
at the fifth and ninth runs there is a sharp reduction in the correlation
while the tenth run shows a continuous decline in the correlation. Figure (17)
and Figure (17) show the density of the marginal posterior distributions for
both parameter $\lambda$ and $\tau$ respectively. In each plot, the density
summarizes the posterior sample. The thick horizontal bar at the bottom shows
the posterior $95\%$ credible interval and the maximum value is the posterior
mean. The posterior estimate for $\lambda$ is $\hat{\lambda}=0.028$, with the
credible interval of $(0.0275,0.0285)$. The posterior estimate for $\tau$ is
$\hat{\tau}=158.48$ with the credible interval of $(150,170)$.

####  4.4.3 Reconstruction result for Mouse Data

Now, the reconstruction result for the mouse injected with $100\mu 1/7.5$ MBq
$[^{99m}Tc]$ MDP, to X-ray the bone imaging is presented in Figure (19) and
Figure (19) for splitting EP and MCMC respectively. Figure (19)(a) shows the
estimate of the true image reconstructed from the noise by EP-ADMM. It can be
seen that the image has a clear background with slightly rough edges on the
image. The residual error shown in Figure (19)(b) comparing the estimates with
the data, shows no randomness. Similar comment can be made about the residual
shown in Figure (19)(c&d). It is possible to reduce the patterns in the
residual by adjusting the value of $\alpha_{\lambda}$ in the hyperparameter
$\lambda$. This may lead to introducing another noise into the posterior mean
estimate. Figure (19)(c) and Figure (19)(d) show the estimates of the
parameter which converges to around $165.44$ and $0.019$ produced by EP-MCMC.
Figure (19) on the other hand, shows a result produced by MCMC which is seen
to have slightly similar results with EP-ADMM with a difference of blurry
background. This might be as a result of long chain needed to reach
convergence. Also, in Figure (19)(e), it can be seen that the samples are
highly correlated even at $1500$ iteration. To reduce the correlation, the
number of $n$ samples to thin might be increased that is, the number of
samples to be discarded for each kept sample. Figure (19)(c) shows the trace
plot for the parameter $\tau$. It can be seen that the acceptance rate is
quite low, this might be due to the choice of proposal distribution. Also, in
Figure (19)(e) the $\sigma_{\epsilon}$ chains produced by MCMC have high
correlation which doesn’t show any sign of reduction. On the contrary, the
observations $\sigma_{x}$ of MCMC exhibit independence in the chains. However,
all the parameters from EP-MCMC and MCMC converge to a very similar value with
slight differences. We further establish the convergence of the $\tau$ and
$\lambda$ for EP-MCMC.

Figure 20: The Autocorrelation plot of the mouse injected with Mdp reagent;
the EP-MCMC algorithm was run $10$ time each of length $1000$. The plot shows
independence of chains at each replication.

Figure 21: The plot shows a reduction in dependence of chains at first three
replications but at the fifth and ninth there is no correlation and the last
replication shows a decline in correlation.

Figure 22: Iterative PSRF Plot for $\lambda$ in Mouse image data (from $m=10$
parallel sequence and $n=1000$). About $800$ iterations, the convergence
starts till the end of the iteration.

Figure 23: Iterative PSRF Plot for $\tau$ in Mouse image data (from $m=10$
parallel sequence and $n=1000$). About $600$ iterations, the convergence
starts till the end of the iteration

Table 3: The mean estimate of the ten EP-MCMC chains for prior $\tau$ and $\lambda$ of mouse data parameter | Posterior-Mean | $\hat{R}$ | $\hat{V}$ | $W$ | $B/n$ | $\sigma_{+}^{2}$ |
---|---|---|---|---|---|---|---
$\tau$ | $165.44$ | $1.02$ | $56.84$ | $56.58$ | $0.00031$ | $56.83$ |
$\lambda$ | $0.019$ | $1.00$ | $1.21e^{-08}$ | $1.21e^{-08}$ | $1.14e^{-11}$ | $1.21e^{-08}$ |

Figure (23) and Figure (23) give a visual interpretation of the values
presented. The shrinking factor for $\lambda$ in Figure (23) was above $1$ at
the initial stage of the iteration, but on getting to around $800$ iterations
the equality between the between-chain and within-chain becomes so evident.
This indicates that $1000$ iterations is sufficient to provide convergence for
parameter $\lambda$. Likewise, the shrinking factor for the $\tau$ at the
first $100$ iterations is $10$ which means the chains still need a lot of
iterations to bring reduction in the shrinking factor. At about $600$
iterations, the shrinking factor is on the acceptable line of $1$. Table (3)
shows some statistical inference to establish the convergence properties of
the EP-MCMC parameters. The potential scale reduction faction of $\tau$ is
$1.02$ which is still within the acceptable range of $1$. This can be
interpreted as the chain represents the target distribution of the parameter
$\tau$. Also, the PSRF for $\lambda$ is exactly $1$. The within-chain,
between-chain variances, and pooled variance for $\tau$ are $W=56.58$,
$B/n=3.1\times 10^{-4}$, and $\hat{V}=56.84$ respectively. Similarly, the
within-chain, between-chain variances, and pooled variance for $\lambda$ are
$W=1.21\times 10^{-8}$, $B/n=1.14\times 10^{-11}$, and $\hat{V}=1.21\times
10^{-8}$ respectively. It can be seen that according to the interpretation by
Gelman and Rubin, (1992), the equality of pooled variance $\hat{V}$ and $W$
for both parameters indicates the existence of convergence.

Figure 24: Marginal posterior distribution of the parameter $\lambda$ for
Mouse data. The tick horizontal bar signifies the posterior $95\%$ credible
interval and the maximum value at $0.019$ is the posterior mean.

Figure 25: Marginal posterior distribution of the parameter $\tau$ for Mouse
data. The tick horizontal bar signifies the posterior $95\%$ credible interval
and the maximum value at about $166$ is the posterior mean.

Figure 26: (a): dmsa_1h_mouse_results from SEP, (b): relative error (c):
estimate $\tau$ is $58.9$, (d): estimate $\lambda$ is $0.179$

Figure 27: (a): MCMC Dmsa (b)&(c): error & residual, (d)&(e): the estimate of
$\sigma_{x}$ & $\sigma_{\epsilon}$ converges at $10$ and $3.0$

Figure (25) and Figure (25) show the density of the marginal posterior
distributions for both parameter $\lambda$ and $\tau$ respectively. In each
plot, the density summarizes the posterior sample. The thick horizontal bar at
the bottom shows the posterior $95\%$ credible interval and the maximum value
is the posterior mean. The posterior estimate for $\lambda$ is
$\hat{\lambda}=0.019$, with the credible interval of $(0.01885,0.0198)$. The
posterior estimate for tau is $\hat{\tau}=165.44$ with the credible interval
of $(100,250)$.

####  4.4.4 Reconstruction result for DMSA Data

Now, moving to the reconstruction of the mouse injected with DMSA. The results
produced by splitting EP and MCMC are quite similar in Figure (27)(a), Figure
(27)(b), Figure (27)(a), and Figure (27)(b) respectively. Given that we
require clear and smooth image reconstruction, EP-ADMM produces a smooth image
reconstruction with a clear and white background. In contrast, MCMC produces a
white and clear background with sharp surface coupled with slightly rough
edges. The residual errors shown in Figure (27)(b) and Figure (27)b comparing
the estimates with the data, show similar pattern. Adjusting the value of
$\alpha_{\lambda}$ and $\alpha_{\tau}$ in the hyperparameters $\lambda$ and
$\tau$ respectively may produce random errors. In Figure (27)(d), the mean
estimate of the posterior for $\lambda$ across the $10$ replication is
$0.179$, which can be seen in the trace-plot. Figure (27)(c) does not show any
pattern of convergence even after $1000$ iterations. This however may be due
to the choice of proposal distribution. Figure (31) and Figure (31) show the
density of the marginal posterior distributions for both parameter $\lambda$
and $\tau$ respectively. In each plot, the density summarizes the posterior
sample.

Figure 28: Iterative PSRF Plot for $\lambda$ in Mouse image data (from $m=10$
parallel sequence and $n=1000$). After about $800$ iteration convergence
starts till the end of the iteration.

Figure 29: Iterative PSRF Plot for $\tau$ in Mouse image data (from $m=10$
parallel sequence and $n=1000$). About $200$ iteration convergence starts till
the end of the iteration.

Figure 30: Marginal posterior distribution of the parameter $\lambda$ for DMSA
data. The tick horizontal bar signifies the posterior $95\%$ credible interval
and the maximum value at $0.179$ is the posterior mean.

Figure 31: Marginal posterior distribution of the parameter $\tau$ for DMSA
data. The tick horizontal bar signifies the posterior $95\%$ credible interval
and the maximum value at about $59$ is the posterior mean.

Table 4: The mean estimate of the ten EP-MCMC chains for prior $\tau$ and
$\lambda$
of DMSA data parameter | Posterior-Mean | $\hat{R}$ | $\hat{V}$ | $W$ | $B/n$ | $\sigma_{+}^{2}$ |
---|---|---|---|---|---|---|---
$\tau$ | $58.82$ | $1.00$ | $59.89$ | $59.87$ | $0.75$ | $59.89$ |
$\lambda$ | $0.179$ | $1.00$ | $7.60e^{-07}$ | $7.60e^{-07}$ | $8.93e^{-10}$ | $7.60e^{-07}$ |

The thick horizontal bar at the bottom shows the posterior $95\%$ credible
interval and the maximum value is the posterior mean. The posterior estimate
for $\lambda$ is $\hat{\lambda}=0.179$, with the credible interval of
$(0.176,0.182)$. The posterior estimate for tau is $\hat{\tau}=58.82$ with the
credible interval of $(1,102)$.

####  4.4.5 Reconstruction result of four Circles Data

Finally, Figure (32) shows the original image to be reconstructed. Figure
(34)(a) presents the reconstruction result for the image of four circles. It
can be seen that EP-ADMM produces a clear and sharp reconstruction image,
while MCMC produces a scaly and gray background as shown in Figure (34)(a).
Comparing the parameter estimates of MCMC with EP-MCMC, the parameter
estimates of $\sigma_{x}$ in MCMC is $160$ while $\tau$ in EP-MCMC is about
$200$, similarly the parameter estimates of $\sigma_{\epsilon}$ is about $0.5$
while that of EP-MCMC is $0.023$. The residual error shown in Figure (34)(b)
comparing the estimates with the data, shows a very clear pattern. Random
residual as shown in Figure (11)(b&c) can be produced by using an appropriate
values of the hyperparameters $\lambda$ and $\tau$. However, given our primary
aim to produce a constant reconstruction, this remains a little concern and it
can be further checked in the future. The computed potential scale reduction
factors $\hat{R}$ for $\tau$ is $1.1$, upper C.I is $1.12$, and $\lambda$ is
$1$ with upper C.I as $1$ are shown in Table (5).

Figure 32: $(a)$ Original image of circles to be reconstructed: $\mathcal{X}$,
$(b)$ the mean: $\mathcal{G}\mathcal{X}$ and $(c)$ Noisy data:$\mathcal{Y}$

Figure 33: (a): The reconstruction image of four circles by EP-ADMM, (b):
relative error (c): estimates $\tau$ converges at $200$ from EP-MCMC, (d):
estimates $\lambda$ converges at $0.023$ from EP-MCMC

Figure 34: (a): The reconstruction image of four circles produced by MCMC
(b)&(c): error & residual, (d): the estimate of $\sigma_{x}$ converges at
$150$, (e): the estimates of $\sigma_{\epsilon}$ converges to about $0.5$

Figure 35: The Autocorrelation plot of the image of four circles; the EP-MCMC
algorithm was run $10$ time each of length $1000$. The plot shows independence
of chains at each replication.

Figure 36: The plot shows a strong dependence of chains at first four
replications but at the fifth and ninth there is no correlation and the last
replication shows a decline in correlation

Here, the potential scale reduction factor for $\lambda$ is exactly $1$ which
indicates strong convergence across the $10$ runs. Figure (36) shows the
correlation between the chains of $\lambda$ for each run. It can be seen that
there is no autocorrection in the chains for each of the runs, this indicates
that there is a convergence at every run which also contributes to potential
scale reduction factor of $1$. Similarly, Figure (36) shows an autocorrelation
plot of $\tau$ for the $10$ runs at the initial stage of the runs. There is a
high correlation among the chain of the first four runs, but at the fifth and
ninth runs there is sharp reduction in the correlation while the tenth run
shows a continuous decline in the correlation.

Figure 37: Iterative PSFR Plot for $\lambda$ in image of four circle (from
$m=10$ parallel sequence and $n=1000$). After about $200$ iteration
convergence starts with a bump at the middle but got stabilized around $300$,
till the end of the iteration.

Figure 38: Iterative PSFR Plot for $\tau$ in image of four circles (from
$m=10$ parallel sequence and $n=1000$). About $200$ iteration convergence
started and diverged in the middle but later converged at about $800$ till the
end of the iteration.

Table 5: The mean estimate of the ten EP-MCMC chains for prior $\tau$ and
$\lambda$
of Circles parameter | Posterior-Mean | $\hat{R}$ | $\hat{V}$ | $W$ | $B/n$ | $\sigma_{+}^{2}$ |  |
---|---|---|---|---|---|---|---|---
$\tau$ | $200$ | $1.10$ | $112.06$ | $110$ | $2.16$ | $112.06$ |  |
$\lambda$ | $0.023$ | $1.00$ | $9.84e^{-09}$ | $9e^{-09}$ | $1.41e^{-11}$ | $9.84e^{-09}$ |  |

##  5 Conclusion and Future Direction

This paper presented a hybrid algorithm called Splitting Expectation
Propagation for approximating hierarchical Bayesian models. Our algorithm
focused on the problems of instability and intractability in EP especially at
the refinement stage where it often fails. Also, this paper generalizes EP to
handle hierarchical Bayesian models which poses a problem due to rigorous
moment matching especially when the prior and likelihood family distributions
differ. The modification was employed from stochastic and deterministic
standpoints.

We adopted the technique by John et al., (2011) which uses Monte Carlo
integration for the intractable integration of the EP. SEP was used for the
image reconstruction and it competed well with the MCMC in that SEP produced
well approximated parameter estimates. EP-ADMM was used to handle the image
reconstruction and EP-MCMC was used for parameter estimation. According to the
results presented EP-ADMM produced very clear, sharp and white background. EP-
MCMC achieved low variances when compared to the ordinary MCMC. Splitting EP
in these viewpoints is the umbrella over both EP-MC, EP-ADMM, and EP-MCMC.
Splitting EP was able to solve the problems of intractability and
inflexibility of EP to hierarchical Bayesian models. Splitting EP can
therefore be seen as an alternative EP for image analysis in high dimensional
space.

The main limitation of EP-MC is that it required large number of samples to
fulfill the central limit theorem due to the Monte Carlo integration. As a
result of the large samples, EP-MC might be slow to converge. Moreover, the
limitation of the EP-MCMC differs from the ordinary MCMC in that the inner
loop EP-ADMM that contributes to the EP-MCMC for the hyperparameter is
efficiently fast. On the contrary, ordinary MCMC for the prior $\mathcal{X}$
will be slow resulting in slow convergence in total. As a future direction, we
look forward to handling the hyper-parameters with Variational Inference. This
method would make perfect sense as the approximate posterior from the prior
$\mathcal{X}$ will be used as the auxiliary posterior distribution factorized
over both the hyper-parameters $\lambda$ and $\tau$ which can be seen as a
Coordinate Ascent Variational Inference (CAVI) David et al., (2018).

## Appendix

### Full derivation of SEP

#### Incorporating the priors into $q$

We first incorporated priors on $\tau$ and $\lambda$ as factors. These factors
have the same functional form as

$q(\mathcal{X},\lambda,\tau)=\bigg{[}\prod_{i=1}^{m}\prod_{j=1}^{n}\mathcal{N}(\mathcal{X}_{ij}\big{|}mx_{i,j},vx_{i,j})\bigg{]}\exp(\lambda|\alpha_{\lambda})\exp(\tau|\alpha_{\tau})$
(35)

The first update rules $\alpha_{\lambda}^{new}$ and $\alpha_{\tau}^{new}$ for
$q$ is obtained by SSEP method which will be discussed in the algorithm below.
The rest of the factors are sequentially incorporated into Equation (35) which
is also updated in a similar manner. One difficulty encountered when applying
the update rules in Minka, (2001) is that the normalizer does not have a
closed form. This brings about a uniqueness in SSEP which treats $Z$ as
follows;

$Z=\int\mathcal{N}(\mathcal{L}\mathcal{X}_{ij}\hskip 1.42271pt|\hskip
1.42271pt0,\hskip 1.42271pt\tau)\hskip
1.42271ptq(\mathcal{X},\lambda,\tau)\hskip 1.42271ptd\mathcal{X}\hskip
1.42271ptd\lambda\hskip 1.42271ptd\tau$
$=\int\mathcal{N}(\mathcal{L}\mathcal{X}_{ij}\hskip 1.42271pt|0,\hskip
1.42271pt\tau)\hskip 1.42271pt\exp(\tau\hskip 1.42271pt|\hskip
1.42271pt\alpha_{\tau})\hskip 1.42271ptd\mathcal{X}\hskip 1.42271ptd\tau$ (36)

In Equation (36), the integral involves $\tau$. We adopted the method by John
et al., (2011). The normalizing factor $Z$ in Equation (36) is solved for
$\tau$ as follows,

$Z_{\tau}=\frac{1}{K}\sum_{k=1}^{K}\mathcal{N}(\mathcal{L}\mathcal{X}_{ij}\hskip
1.42271pt|\hskip 1.42271pt0,\hskip 1.42271pt\tau^{(k)})\hskip
1.42271pt\exp(\tau^{(k)}\hskip 1.42271pt|\hskip 1.42271pt\alpha_{\tau}),$

$\tau$ is generated from exponential distribution while $\mathcal{X}$ is
initialized by $\mathcal{Y}$ to compute $\mathcal{LX}$ which is in general the
goal of this work. So the Monte Carlo integration of Equation (36) with
respect to $\mathcal{X}$ becomes

$Z_{x}=\frac{1}{K}\sum_{k=1}^{K}\mathcal{N}(\mathcal{Y}|\mathcal{GX}_{ij},\hskip
1.42271pt\lambda)\hskip
1.42271pt\mathcal{N}(\mathcal{X}^{(k)},mx_{ij},vx_{ij})\hskip 1.42271pt$

Here, $\mathcal{X}$ is generated from any distribution student-t, Normal
distribution, or Uniform distribution from which $\mathcal{LX}$ is recomputed.
Here, we write the exact likelihood that needs to be processed multiple times
i.e $N\times D$ times, where $N$ and $D$ represent the number of rows and
columns respectively.

$p(L\mathcal{X}|\tau)=\prod_{i=1}^{m}\prod_{j=1}^{n}\mathcal{N}(L\mathcal{X}_{ij}\big{|}0,\tau)$
(37)

To approximate the posterior with SSEP algorithm, the posterior distribution
can be factorized as follows

$q(\mathcal{X},\lambda,\tau)=\bigg{[}\prod_{i=1}^{m}\prod_{j=1}^{n}\mathcal{N}(\mathcal{X}_{ij}\big{|}mx_{i,j},vx_{i,j})\bigg{]}\exp(\lambda|\alpha_{\lambda})\exp(\tau|\alpha_{\tau})$

The exact likelihood term is a joint distribution and is denoted as

$t_{i,j}(\tau)=\mathcal{N}(L\mathcal{X}_{ij}\big{|}0,\tau)$

Then, we approximate the exact likelihood term by choosing approximating term

$\tilde{t}_{i,j}(\tau)=\exp(\tau|\tilde{\alpha}_{i,j})$

We choose term to refine by removing $\tilde{t}_{i,j}$ from $q$ to compute the
cavity distribution on $\mathcal{X}_{i,j}$ and $\tau$

$q_{-i,j}(\tau)=\exp(\tau|\alpha_{-i,j})$
$\alpha_{-i,j}=\alpha_{i,j}-\tilde{\alpha}_{i,j}$

After this, we perform the moment matching between $q(\tau)$ and a normalizing
version of $t_{i,j}(\tau)q_{-i,j}(\tau)$. It is computed as follows

$\operatorname*{argmin}_{q}KL(\tilde{p}_{i,j}\big{|}\big{|}q_{i,j})$

where

$\tilde{p}_{i,j}(\tau)=\frac{t_{i,j}(\tau)q_{-i,j}(\tau)}{\int
t_{i,j}(\tau)q_{-i,j}(\tau)d\tau}$ (38)

Following Equation (38), the integral contains an intractable terms resulting
from $\nabla_{\omega}\log{Z}$ where

$Z=\int_{\tau}t_{i,j}(\tau)q_{-i,j}(\tau)d\tau$ (39)

Our goal is to make a stochastic approximation of this gradient.

$\nabla_{\omega}\log{Z}=\frac{1}{Z}\nabla_{\omega}\int_{\tau}t_{i,j}(\tau)q_{-i,j}(\tau)d\tau$
(40)
$\nabla_{\omega}\log{Z}=\frac{1}{Z}\int_{\tau}t_{i,j}(\tau)\nabla_{\omega}q_{-i,j}(\tau)d\tau$
(41)

Equation (41) can be written as

$\nabla_{\omega}\log{Z}=\frac{1}{Z}\int_{\tau}t_{i,j}(\tau)q_{-i,j}(\tau)\nabla_{\omega}\log
q_{-i,j}(\tau)d\tau$ (42)

referring to Equation (38) as the posterior, we have Equation (42) as an
expectation

$\nabla_{\omega}\log{Z}=E_{\hat{p}(\tau)}\big{[}\nabla_{\omega}\log
q_{-i,j}(\tau)\big{]}$ (43)

we use the identity
$\nabla_{\omega}q_{-i,j}(\tau)=q_{-i,j}(\tau)\nabla_{\omega}\log
q_{-i,j}(\tau)$. We can stochastically approximate the expectation using Monte
Carlo integration as follows;

$\nabla_{\omega}\log{Z}=\frac{1}{K}\sum_{k=1}^{K}\nabla_{\omega}\log
q_{-i,j}(\tau)$ (44)

Now, we substitute for $q_{-i,j}(\tau)$

$\nabla_{\alpha}\log{Z}=\frac{1}{K}\sum_{k=1}^{K}\nabla_{\alpha}\log\big{(}\exp(\tau|\alpha_{-i,j})\big{)}$
(45)

we compute the updates as follows;

$\nabla_{\alpha_{-ij}}\log{Z}=\frac{1}{K}\sum_{k=1}^{K}\big{(}\tau^{(k)}\big{)}$
(46)

The update rule is given in equation (47) below;

$\alpha_{\tau}^{new}\approx\alpha_{-ij}+v_{\alpha_{-ij}}\zeta_{\alpha_{-ij}}$
(47)

Where the mean $\frac{1}{K}\sum_{k=1}^{K}\tau^{(k)}$ is denoted as
$\zeta_{\alpha_{-ij}}$ and variance as $v_{\alpha_{-ij}}$.

#### Incorporating the likelihood factors into $q$

Now, the $N\times D$ factors are sequentially incorporated for the likelihood
in Equation (37). However, approaching this with conventional update rule in
EP is difficult to compute because it requires integration of each likelihood
factor. After incorporating all the factors in Equation (37), SSEP
sequentially incorporates the $N\times D$ factors for the likelihood Equation
(37) provided below in Equation (48)

$Z_{\lambda}=\int\mathcal{N}(\mathcal{Y}\hskip 1.42271pt|\hskip
1.42271pt\mathcal{G}\mathcal{X},\hskip 1.42271pt\lambda)\hskip
1.42271ptq(\mathcal{X},\lambda,\tau)\hskip 1.42271ptd\mathcal{X}\hskip
1.42271ptd\lambda\hskip 1.42271ptd\tau$ $=\int\mathcal{N}(\mathcal{Y}\hskip
1.42271pt|\hskip 1.42271pt\mathcal{G}\mathcal{X},\hskip
1.42271pt\lambda)\hskip 1.42271pt\exp(\lambda\hskip 1.42271pt|\hskip
1.42271pt\alpha_{\lambda})d\lambda.$ (48)

Then the Monte Carlo approximation for Equation (48) is

$\frac{1}{K}\sum_{k=1}^{K}\mathcal{N}(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda^{(k)})\exp(\lambda^{(k)}|\alpha_{\lambda}).$
(49)

Now we compute the gradient of the normalizing factor as follows;

$\nabla_{\alpha_{\lambda}}\log{Z}=\frac{1}{Z}\int_{\lambda}t_{i,j}(\lambda)\nabla_{\alpha_{\lambda}}q_{-i,j}(\lambda)d\lambda$
$\nabla_{\alpha_{\lambda}}\log{Z}=\frac{1}{Z}\int_{\lambda}t_{i,j}(\lambda)q_{-i,j}(\lambda)\nabla_{\alpha_{\lambda}}\log
q_{-i,j}(\lambda)d\lambda$
$\nabla_{\alpha_{\lambda}}\log{Z}=\frac{1}{Z}\int_{\lambda}\mathcal{N}(\mathcal{Y}|\mathcal{G}\mathcal{X},\lambda)q_{-i,j}(\lambda)\nabla_{\alpha_{\lambda}}\log\exp(\lambda|\alpha_{\lambda})d\lambda$
(50)

Equation (50) can be approximated as follows;

$E_{\hat{p}}[\lambda]=\sum_{k=1}^{K}\nabla_{\alpha_{\lambda}}\log\exp(\lambda|\alpha_{\lambda})$

Similarly, for $\lambda$

$\nabla_{\omega}\log{Z}=\frac{1}{K}\sum_{k=1}^{K}\nabla_{\omega}\log\big{(}\exp(\lambda|\alpha^{\lambda}_{-i,j})\big{)}$
$\nabla_{\alpha^{\lambda}_{-ij}}\log{Z}=\frac{1}{K}\sum_{k=1}^{K}\big{(}\lambda^{(k)}\big{)}$

We denote the mean $\frac{1}{K}\sum_{k=1}^{K}\lambda^{(k)}$ as
$\zeta_{\lambda_{-ij}}$ and variance as $v_{\alpha^{\lambda}_{-ij}}$. The
update rule for $\lambda$ is given in Equation (51) below;

$\alpha_{\lambda}^{new}\approx\alpha^{\lambda}_{-ij}+v_{\alpha^{\lambda}_{-ij}}\frac{1}{K}\sum_{k=1}^{K}\lambda^{(k)}$
(51)

The cavity of the prior over $\mathcal{X}$ is

$q_{-ij}(\mathcal{X})=\mathcal{N}(x_{ij};m_{-ij},v_{-ij})$

We choose the approximate posterior for $\mathcal{X}$ as

$q(\mathcal{X})=\mathcal{N}(x_{ij};m_{x},v_{x})$

Then the approximate term is chosen from the Gaussian family

$\tilde{t}_{ij}(\mathcal{X})=\mathcal{N}(\mathcal{X}_{ij};\tilde{m}_{ij},\tilde{v}_{ij})$

In order to update the cavity distribution we remove $\tilde{t}_{ij}$ from $q$
and update its parameters as follows;

$v_{-i,j}^{-1}=vx_{i,j}^{-1}-\tilde{v}_{i,j}^{-1}$
$m_{-i,j}=mx_{i,j}+v_{i,j}\tilde{v}_{i,j}(mx_{i,j}-\tilde{m}_{i,j})$

Our goal is to make a stochastic approximation of this gradient.

$\nabla_{\omega}\log{Z}=\frac{1}{Z}\nabla_{\omega}\int_{\mathcal{X}_{i,j},\tau}t_{i,j}(L\mathcal{X}_{i,j},\tau)q_{-i,j}(\mathcal{X}_{i,j},\tau)d\mathcal{X}_{i,j}d\tau$
(52)

Replace $q_{-i,j}$ in Equation (52) with $q_{i,j}/\tilde{t}_{i,j}$ then it
becomes

$\nabla_{\omega}\log{Z}=\frac{1}{Z}\int_{\mathcal{X}_{i,j},\tau}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}\nabla_{\omega}{q_{i,j}(\mathcal{X}_{i,j},\tau)}d\mathcal{X}_{i,j}d\tau$
$\nabla_{\omega}\log{Z}=\frac{1}{Z}\int_{\mathcal{X}_{i,j},\tau}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}q_{i,j}(\mathcal{X}_{i,j},\tau)\nabla_{\omega}\log{q_{i,j}(\mathcal{X}_{i,j},\tau)}d\mathcal{X}_{i,j}d\tau$

At this juncture, we proceed with the use of Monte Carlo Integration

$\frac{\delta_{1}}{\delta_{0}}=\sum_{n=1}^{N}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}\nabla_{\omega}\log{q_{i,j}(\mathcal{X}_{i,j},\tau)}\bigg{/}\sum_{n=1}^{N}\frac{t_{i,j}(L\mathcal{X}_{i,j},\tau)}{\tilde{t}_{i,j}(\mathcal{X}_{i,j},\tau)}$

when taking $\mathcal{X}$ into consideration then we fix $\tau$ and vice
versa. So we update the parameter by taking the gradient step

$\omega^{t+1}=\omega^{t}+lr*\frac{\delta_{1}}{\delta_{0}}$ (53)

### Derivative of EP-ADMM

We introduce the Alternating direction method of multiplier (ADMM) algorithm
to update the approximate posterior parameters.

$\displaystyle\begin{split}\text{minimize}\hskip
14.45377pt\text{KL}(\tilde{p}(x_{ij}||q(x)_{ij}))\end{split}$ $\displaystyle\
\text{subject to}:m_{ij}\geq a;v_{ij}\geq b$

where $a$ and $b$ are constants. Then updating according to Minka, (2001) is
as follows;

$m_{x}=\nabla_{m_{-ij}}\log Z_{x}+\alpha+\rho(m_{-ij}-a)$

and

$v_{x}=\nabla_{v_{-ij}}\log Z_{x}+\beta+\rho(v_{-ij}-b)$

where

$q(x)=\mathcal{N}(\mathcal{X};m_{x},v_{x})$

then the tilted distribution is

$\tilde{p}(x_{ij})=\frac{t_{ij}(x)q_{-ij}(x)}{\int_{\mathcal{X}_{ij}}t_{ij}(x)q_{-ij}(x)dx_{ij}}$

Now the normalizing factor is

$Z_{x}=\int_{x_{ij}}\mathcal{N}(\mathcal{Y}_{ij};G\mathcal{X}_{ij},\lambda)\mathcal{N}(\mathcal{X}_{ij};m_{-ij},v_{-ij})dx_{ij}$

By integrating out $x$, $Z_{x}$ becomes

$\int\exp\biggl{\\{}-\frac{1}{2}\Big{(}y_{ij}-gx_{ij}\Big{)}^{2}\lambda^{-1}+\Big{(}x_{ij}-m_{-ij}\Big{)}^{2}v_{-ij}^{-1}\biggr{\\}}dx_{ij}$

Completing the squares, we have

$-\frac{1}{2}\Bigg{[}x_{ij}-\frac{gy_{ij}\lambda^{-1}+m_{-ij}v_{-ij}^{-1}}{g^{2}\lambda^{-1}+v_{-ij}^{-1}}\Bigg{]}^{2}\bigg{(}g^{2}\lambda^{-1}+v_{-ij}^{-1}\bigg{)}+\frac{1}{2}\frac{\bigg{(}gy_{ij}\lambda^{-1}+m_{-ij}v_{-ij}^{-1}\bigg{)}^{2}}{g^{2}\lambda^{-1}+v_{-ij}^{-1}}$

Then the normalizing factor of $x$ is

$Z_{x}=\mathcal{N}(y_{ij};m_{y},v_{y})$

where

$m_{y}=\frac{g\lambda^{-1}m_{-ij}v_{-ij}^{-1}}{\lambda^{-1}-g^{2}\lambda^{-2}\bigg{(}g^{2}\lambda^{-1}+v_{-ij}^{-1}\bigg{)}^{-1}}=gm_{-ij}$

and

$\frac{1}{\lambda^{-1}-g^{2}\lambda^{-2}\bigg{(}g^{2}\lambda^{-1}+v_{-ij}^{-1}\bigg{)}^{-1}}=v_{-ij}g^{2}+\lambda$

then updating the parameters $m_{x}$ and $v_{x}$ of $q(\mathcal{X})$ as
follows according to Minka, (2001) with a method of multipliers, we have

$m^{new}_{x}=m_{-ij}+v_{-ij}\frac{y_{ij}-gm_{-ij}}{v_{-ij}g^{2}+\lambda}g+\alpha+\rho(m_{i}-a)$
$v_{x}^{new}=\frac{v_{-ij}\lambda}{v_{-ij}g^{2}+\lambda}+\beta+\rho(v_{i}-b)$
$\alpha^{new}=\alpha^{k}+\rho(m_{x}^{new}-a)\hskip 14.45377pt\text{and}\hskip
14.45377pt\beta^{new}=\beta^{k}+\rho(v_{x}^{new}-b)$ (54)

## References

  * Brooks and Gelman, (1998) Brooks, S. P. and Gelman, A. (1998).  General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics.  7:4(434 - 455).
  * David et al., (2018) David, M. B., Kucukelbir, A., and McAuliffe, J. D. (2018).  Variational inference: A review for statisticians arXiv:1601.00670v9 [stat.CO].  (1–41).
  * Gelman and Rubin, (1992) Gelman, A. and Rubin, D. (1992).  Inference from iterative simulation using multiple sequences, Statistical Science.  7(457–511).
  * Gelman et al., (2014) Gelman, A., Vehtari, A., Jylanki, P., Robert, C., and Chopin, N. (2014).  Expectation propagation as a way of life. arXiv:1412.4869v1.  (1–29).
  * Gelman et al., (1996) Gelman, A., F., Bois, Y., and Jiang, J. (1996).  Physiological pharmacokinetic analysis using population modelling and informative prior distributions.  Journal of the American Statistical Association, 91:1400 – 1412\.
  * Graves, (2011) Graves, A. (2011).  Practical variational inference for neural networks. In advances in Neural Information Processing System 24.  (pp. 2348-2356.).
  * John et al., (2011) John, P., David, B., and Michael, J. (2011).  Variational bayesian inference with stochastic search.  (1–8).
  * Jose and Ryan, (2015) Jose, M. H. and Ryan, A. P. (2015).  Probabilistic backpropagation for scalable learning of bayesian neural networks. arXiv:1502.05336v2 [stat.ML].  (Harvard University, Cambridge, MA 02138 USA).
  * Jylanki et al., (2011) Jylanki, P., Vanhatalo, J., and Vehtari, A. (2011).  Robust gaussian process regression with a student-t likelihood. Journal of Machine Learning Research.  (12, 3227–3257).
  * Maria et al., (2016) Maria, G., Eleftherios, F., Fragogeorgi, E., and Loudos, G. (2016).  Characterization of ”y-eye: a low-cost benchtop mouse-sized gamma camera for dynamic and static imaging studies”. DOI: 10.1007/s11307-016-1011-4 [World Molecular Imaging Society].  (398–407).
  * Matthias and Bangti, (2014) Matthias, G. and Bangti, J. (2014).  Expectation propagation for nonlinear inverse problems with an application to electrical impedance tomography. Journal of Computational Physics.  (513 - 535).
  * Minka, (2001) Minka, T. (2001).  Expectation propagation for approximate bayesian inference.  In In proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, ed. J. Breese and D Koller, pages 362–369.
  * Minka, (2004) Minka, T. (2004).  Power expectation propagation. Technical report, Microsoft Research, Cambridge. (1).
  * Minka, (2005) Minka, T. (2005).  Divergence measures and message passing. Technical report, Microsoft Research, Cambridge. (1).
  * Seeger, (2008) Seeger, M. (2008).  Bayesian inference and optimal design for the sparse linear model. Journal of Machine Learning Research.  (9, 759–813).
  * Vehtari et al., (2019) Vehtari, A., Gelman, A., Sivula, T., Jylanki, P., Tran, D., Sahai, S., Blomstedt, P., Cunnigham, J. P., Schminovich, D., and Robert, C. (2019).  Expectation propagation as a way of life: A framework for bayesian inference on partitioned data.  arXiv:1412.4869v5 [stat.CO], pages 1 – 48.
  * Yingzhen et al., (2018) Yingzhen, L., Jose, M., and Richard, E. (2018).  Stochastic expectation propagation. arXiv:1506.04132v2[Stat.ML].