# Variational Inference for
the Latent Shrinkage Position Model

Xian Yao Gwee<EMAIL_ADDRESS>Isobel Claire Gormley
<EMAIL_ADDRESS>Michael Fop<EMAIL_ADDRESS>
Isobel Claire Gormley and Michael Fop have contributed equally to this work.

###### Abstract

The latent position model (LPM) is a popular method used in network data
analysis where nodes are assumed to be positioned in a $p$-dimensional latent
space. The latent shrinkage position model (LSPM) is an extension of the LPM
which automatically determines the number of effective dimensions of the
latent space via a Bayesian nonparametric shrinkage prior. However, the LSPM
reliance on Markov chain Monte Carlo for inference, while rigorous, is
computationally expensive, making it challenging to scale to networks with
large numbers of nodes. We introduce a variational inference approach for the
LSPM, aiming to reduce computational demands while retaining the model’s
ability to intrinsically determine the number of effective latent dimensions.
The performance of the variational LSPM is illustrated through simulation
studies and its application to real-world network data. To promote wider
adoption and ease of implementation, we also provide open-source code.

Keywords— Network data, Latent position model, Bayesian nonparametric
shrinkage priors, Multiplicative gamma process, Variational inference.

##  1 Introduction

In recent years network data have become increasingly prevalent and
complicated, appearing for example as neural connections in neuroscience
(Aliverti and Durante,, 2019), as trade relationships (Lyu et al.,, 2023) and
as the spread of diseases in epidemiology (Chu et al.,, 2021). Such networks
can have hundreds to thousands of nodes, making them challenging to model
using current statistical tools. Thus, there is a need for statistical network
models that can be used for inference on large networks in a computationally
efficient manner.

The latent position model (LPM, Hoff et al.,, 2002) is a popular network model
due to its ease of interpretation. The LPM assumes that nodes have positions
in a $p$-dimensional latent space, and that the distance between the nodes’
latent positions influence the observed edge formation process. The latent
shrinkage position model (LSPM) (Gwee et al., 2023a, ) extended the LPM to
allow for automatic inference on the number of effective dimensions of the
latent space. This is achieved through a Bayesian nonparametric framework
where a multiplicative truncated gamma process (MTGP) prior is assumed on the
variance of the latent positions. The MTGP prior induces shrinkage on the
variance of the latent positions as the number of dimensions increases. This
approach eliminates the requirement to fit multiple LPMs and then choose
between them using model selection criteria to identify the optimal $p$.

In the Bayesian setting, inference for the LPM and its extensions has
typically relied on Markov chain Monte Carlo (MCMC) techniques, for example in
the clustering setting (Handcock et al.,, 2007), in longitudinal network
models (Chu et al.,, 2021), and in multiview network data (D’Angelo et al.,,
2023), among others. These sampling-based methods, while rigorous, have a high
computational cost. As the number of nodes increases, their computational cost
increases quadratically (Salter-Townshend and Murphy,, 2013; Rastelli et al.,,
2018), making the LPMs challenging to scale to large networks. Since the LSPM
uses Metropolis-within-Gibbs sampling, it also suffers from this issue.

Various efforts have been made to improve computational efficiency when
inferring LPMs via MCMC sampling methods, including the adoption of
Hamiltonian Monte Carlo (Spencer et al.,, 2022), noisy MCMC (Rastelli et al.,,
2018), and the use of case-control approximate likelihood (Raftery et al.,,
2012; Sewell and Chen,, 2015), each having varying degrees of success. Instead
of employing such sampling methods, here a variational inference (VI) approach
is adopted (Jordan et al.,, 1999; Blei et al.,, 2017) in which the posterior
distribution is approximated via an optimization process, reducing the
computational cost of inference. The idea stems from literature where
variational inference has been successfully employed for computationally
efficient estimation of LPMs. For example, Salter-Townshend and Murphy, (2013)
introduced a variational Bayes algorithm for the latent position cluster
distance model, while Gollini and Murphy, (2016) introduced a variational
inference algorithm for a latent position distance model for multilayer
networks. In the context of longitudinal networks, both Sewell and Chen,
(2017) and Liu and Chen, (2023) used variational inference, with a focus on
the latent position projection model and the distance model, respectively.
Additionally, Aliverti and Russo, (2021) presented a stratified stochastic
variational inference method for high-dimensional network factor models. Many
variational algorithms have exhibited empirical success, and recent
theoretical contributions by Yang et al., (2020) and Liu and Chen, (2023) have
provided deeper insights into the statistical guarantees of these methods.

Here, we propose the variational LSPM (V-LSPM), a variational inference
approach for the LSPM to perform posterior inference for large scale networks.
Our simulation studies show improved computational cost compared to MCMC
approaches, while the ability to identify the number of effective dimensions
is retained. We also apply V-LSPM to analyse a cat brain connectivity binary
network with a moderate number of nodes and moderate density, and a relatively
large but sparse binary network describing the nervous system of a worm.

In what follows, Section 2 provides an introduction to the LSPM and outlines a
VI algorithm for posterior inference. Sections 3 and 4 illustrate the
performance of the proposed method with simulation studies and real-world
examples, respectively. Section 5 concludes with a discussion. To aid
practitioners in employing the V-LSPM, the necessary R code is available in
the lspm GitLab repository.

##  2 Variational inference for the latent shrinkage position model

###  2.1 The latent shrinkage position model

Consider a binary network with $n$ nodes, encoded in a $n\times n$ adjacency
matrix, $\mathbf{Y}$, where entry $y_{i,j}=1$ if there is an edge between
nodes $i$ and $j$, and $0$ otherwise, for $i,j=\\{1,\ldots,n\\}$. Self-edges
are not considered.

In the LPM (Hoff et al.,, 2002) each node has an unobserved latent position
$\mathbf{z}_{i}$, in a $p$-dimensional Euclidean latent space. Edges are
assumed independent, conditional on the latent positions of the nodes, giving
the sampling distribution:

$\mathbb{P}(\mathbf{Y}\mid\alpha,\mathbf{Z})=\prod_{i\neq
j}\mathbb{P}(y_{i,j}\mid\alpha,\mathbf{z}_{i},\mathbf{z}_{j})$

where $\mathbf{Z}$ is the matrix of latent positions and $\alpha$ is a global
parameter that captures the overall connectivity level in the network. A
logistic model formulation is then employed to characterize the edge-formation
process, where the log odds of an edge between nodes $i$ and $j$ depends on
the squared Euclidean distance (Gollini and Murphy,, 2016; Gwee et al., 2023a,
) between their latent space positions $\mathbf{z}_{i}$ and $\mathbf{z}_{j}$
i.e.,

$\log\frac{q_{i,j}}{1-q_{i,j}}=\alpha-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|^{2}_{2}$
(1)

where $q_{i,j}$ represents the probability of an edge between nodes $i$ and
$j$.

The LSPM (Gwee et al., 2023a, ) extends the LPM through a Bayesian
nonparametric approach which enables automatic inference of $p$, the number of
effective latent space dimensions necessary to fully describe the network. The
latent positions are assumed to lie in an infinite dimensional space and to
have a Gaussian distribution with zero mean and a diagonal precision matrix
$\mathbf{\Omega}$, whose entries $\omega_{\ell}$ denote the precision of the
latent positions in dimension $\ell$, for $\ell=1,\ldots,\infty$. A shrinkage
prior, the multiplicative truncated gamma process (MTGP), is assumed for the
precision parameters, ensuring shrinkage of the variance in higher dimensions.
Each latent dimension, $h=1,\ldots,\infty$, has an associated shrinkage
strength parameter, $\delta_{h}$ and the precision $\omega_{\ell}$ in
dimension $\ell$ is the cumulative product of $\delta_{1}$ to $\delta_{\ell}$,
i.e., for $i=1,\ldots,n$,

$\displaystyle\mathbf{z}_{i}\sim\text{MVN}(\mathbf{0},\mathbf{\Omega}^{-1})\qquad\mathbf{\Omega}=\begin{bmatrix}\omega_{1}^{-1}&\dots&0\\\
\vdots&\ddots&\\\
0&&\omega_{\infty}^{-1}\end{bmatrix}\qquad\omega_{\ell}=\prod_{h=1}^{\ell}\delta_{h}\quad\mbox{
for }\ell=1,\ldots,\infty.$ (2)

Gamma priors are assumed for the shrinkage strength parameters, where for
dimensions $h>1$ a truncated gamma prior with a left truncation point at
$t_{2}=1$ is assumed to ensure shrinkage, leading to the MTGP prior, i.e.

$\displaystyle\delta_{1}\sim\text{Gam}(a_{1},b_{1}=1)\qquad\delta_{h}\sim\text{Gam}^{\text{T}}(a_{2},b_{2}=1,t_{2}=1)\quad\mbox{
for }h>1.$

Under this MTGP prior, the variances of the latent positions in the higher
dimensions will shrink towards zero causing the dimensions with non-negligible
variance to correspond to the effective dimensions necessary to describe the
network.

The joint posterior distribution of the LSPM is

$\displaystyle\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})\propto\mathbb{P}(\mathbf{Y}\mid\alpha,\mathbf{Z})\mathbb{P}(\alpha)\mathbb{P}(\mathbf{Z}\mid\bm{\delta})\mathbb{P}(\bm{\delta}).$
(3)

In Gwee et al., 2023a , inference on the LSPM relies on MCMC which, while
rigorous, is computationally expensive and makes fitting the LSPM to large
networks with thousands of nodes challenging in practice. Here we propose a
variational approach to inference for the LSPM, providing a fast but
approximate solution, making application of the LSPM to large networks
computationally feasible.

###  2.2 Variational inference for the LSPM

Variational inference (VI) is a computational technique used to approximate
complex, intractable posterior distributions with simpler, tractable ones (see
Blei et al., (2017) for a comprehensive review). Instead of directly computing
the posterior, VI reformulates the problem as an optimization task. This is
done by selecting a family of simpler distributions to optimize the parameters
based on the Kullback-Leibler (KL) divergence, such that the constructed
variational distribution is as close as possible to the true posterior. A
widely-adopted technique for VI is the mean-field approach (Consonni and
Marin,, 2007), which represents the posterior distribution using a completely
factorized form. This allows for the simplification of the expectations needed
for the calculation of the KL divergence.

Although the MTGP prior is nonparametric, it is not implementable in practice
without setting a finite truncation level $p$ on the number of dimensions
fitted. Hence, for fixed $p$, we apply the VI approach to the LSPM posterior
(3) by defining the variational posterior as

$\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})=\mathbb{Q}(\alpha)\prod_{i=1}^{n}\mathbb{Q}(\mathbf{z}_{i})\mathbb{Q}(\delta_{1})\prod_{h=2}^{p}\mathbb{Q}(\delta_{h}),$

where
$\mathbb{Q}(\alpha)=N(\tilde{\mu}_{\alpha},\tilde{\sigma}^{2}_{\alpha}),\mathbb{Q}(\mathbf{z}_{i})=\text{MVN}(\tilde{\mathbf{z}}_{i},\tilde{\bm{\Omega}}^{-1}),\mathbb{Q}(\delta_{1})=\text{Gam}(\tilde{a}_{1},\tilde{b}_{1}),\text{
and
}\mathbb{Q}(\delta_{h})=\text{Gam}^{\top}(\tilde{a}_{2}^{h},\tilde{b}_{2}^{h},1)$.
The variational parameters are distinguished with a tilde for clarity and are
collectively denoted by
$\tilde{\Theta}=\\{\tilde{\mu}_{\alpha},\tilde{\sigma}^{2}_{\alpha},\tilde{\mathbf{z}}_{i},\tilde{\bm{\Omega}}^{-1},\tilde{a}_{1},\tilde{b}_{1},\tilde{a}_{2}^{h},\tilde{b}_{2}^{h}\\}$.
The variational parameter $\tilde{\mathbf{\Omega}}$ has diagonal entries
taking the cumulative products of the expected values of $\delta_{1}$ and
$\delta_{\ell}$ under the variational distributions $\mathbb{Q}(\delta_{1})$
to $\mathbb{Q}(\delta_{\ell})$.

Typically, the KL divergence is used as a measure of closeness to minimise the
distance between the variational posterior
$\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})$ and the true posterior
$\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})$:

$\displaystyle\mathrm{KL}[\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\|\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})]$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})}{\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\mathbf{Y},\alpha,\mathbf{Z},\bm{\delta})}{\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle\qquad+\log\mathbb{P}(\mathbf{Y}).$

Note that minimising the KL divergence is equivalent to maximising the
evidence lower bound (ELBO) with respect to the variational distribution
(Bishop,, 2006).

The KL divergence between the variational posterior and the true posterior for
the LSPM can be written as:

$\displaystyle\mathrm{KL}$
$\displaystyle[\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\|\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})]$
(4) $\displaystyle=$
$\displaystyle\operatorname{KL}\left[\mathbb{Q}(\alpha)\|\mathbb{P}(\alpha)\right]+\sum_{i=1}^{n}\operatorname{KL}\left[\mathbb{Q}\left(\mathbf{z}_{i}\right)\|\mathbb{P}\left(\mathbf{z}_{i}\right)\right]$
$\displaystyle+\operatorname{KL}\left[\mathbb{Q}(\delta_{1})\|\mathbb{P}(\delta_{1})\right]+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]$
$\displaystyle-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left\\{\log\left[\mathbb{P}(\mathbf{Y}\mid\alpha,\mathbf{Z},\bm{\delta})\right]\right\\}+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}+\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle-\frac{np}{2}+\frac{n}{2}\log\operatorname{det}(\mathbf{\Omega}^{-1}\tilde{\mathbf{\Omega}})+\frac{1}{2}\sum_{i=1}^{n}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\frac{n}{2}\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})$
$\displaystyle+\tilde{a}_{1}\left[\psi(\tilde{a}_{1})+\frac{b_{1}}{\tilde{b}_{1}}+1\right]-a_{1}\left[\psi(\tilde{a}_{1})+\log(\frac{b_{1}}{\tilde{b}_{1}})\right]-\log\frac{\Gamma(\tilde{a}_{1})}{\Gamma(a_{1})}$
$\displaystyle+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left\\{\log\left[\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)\right]\right\\}+\log\mathbb{P}(\mathbf{Y}),$

where $\Gamma$ and $\psi$ are the gamma and digamma functions respectively
(see Appendix A for full derivations). As in Gollini and Murphy, (2016), the
expected log-likelihood
$\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\\{\log[\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)]\\}$
is approximated using the Jensen’s inequality:

$\displaystyle\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\\{\log[\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)]\\}$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right]$
$\displaystyle\quad-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log\left(1+\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right)\right]$
$\displaystyle\leq\sum_{i\neq
j}^{n}y_{i,j}\left(\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right]\right)$
$\displaystyle\quad-\log\left(1+\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right]\right)$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad-\log\left\\{1+\frac{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\right.$
$\displaystyle\qquad\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}.$

In contrast to Salter-Townshend and Murphy, (2013), which employs three first-
order Taylor expansions to approximate the log-likelihood, the squared
Euclidean distance formulated in (1) benefits from using just one application
of Jensen’s inequality, yielding a comparatively tighter bound. Jaakkola and
Jordan, (2000) offers an alternative approach to approximate the expected log-
likelihood but, as it requires further approximations, it is more difficult to
compute.

####  2.2.1 Updating the variational parameters

The variational parameters
$\\{\tilde{\mu}_{\alpha},\tilde{\sigma}^{2}_{\alpha},\tilde{\mathbf{z}}_{i}\\}$
are updated by minimising (4) via a bracketing-and-bisection method and a
conjugate gradient routine (see Appendix B for a summary and Appendix C for
the full derivations). However, the variational parameters
$\\{\tilde{a}_{1},\tilde{b}_{1},\tilde{a}_{2}^{h},\tilde{b}_{2}^{h}\\}$ for
$h=2,\ldots,p$ are updated analytically by identifying their optimal
expressions as the gamma (for $h=1$) and truncated gamma (for $h>1$) densities
respectively, with their own shape and rate parameters (see Appendix B for a
summary and Appendix C for the full derivations). The variational inference
algorithm then updates each of the variational parameters in turn, while
holding the others constant.

While Gollini and Murphy, (2016) employed both first and second order Taylor
expansions to obtain an analytical form for updating the variational
parameters, poor performance was observed empirically when the dimension of
the latent space was higher than two. Instead, we adopt the numerical
optimization methods of Salter-Townshend and Murphy, (2013) where the latent
position variational means $\tilde{\mathbf{z}}_{i}$ are optimized via a
conjugate gradient routine (Fletcher,, 1964) while the $\tilde{\mu}_{\alpha}$
and $\tilde{\sigma}^{2}_{\alpha}$ hyperparameters are optimized via a
bracketing-and-bisection algorithm (Burden and Faires,, 1985).

To initialize the estimation algorithm, here, we follow the strategy for
obtaining starting values as taken in Gwee et al., 2023a , with the exception
that the truncation level is set as $p=5$, which was found to work well in
practice for the networks examined here. A well-chosen initial set of
variational parameters is essential for the VI approach due to the potential
presence of multiple local minima in the Kullback–Leibler divergence. To
address the issue of local minima, a range of starting values is used. To do
so, random noise is added to the initial latent positions obtained from
multidimensional scaling, for a pre-specified value of $p$, distributed as
$\text{MVN}(0,r^{2})$ where $r^{2}$ is the five hundredths of the empirical
variance of the multidimensional scaling positions. Ten random initialisations
are considered and the best estimate is selected as the one with the highest
ELBO value.

Convergence of the VI algorithm can be assessed in a number of different ways.
It is typically achieved by monitoring the difference between the ELBO at
successive iterations and stopping the algorithm when this difference is
smaller than a pre-specified threshold. In our implementation, we monitor the
approximate expected log-likelihood across iterations as it often makes the
largest contribution to the evaluation of the ELBO, simplifying the
calculations for the convergence check. In the settings considered here,
convergence was typically achieved after a low number ($<50$) of iterations,
using a threshold of 0.01.

##  3 Simulation studies

The performance of the V-LSPM is assessed on simulated data scenarios with a
computer equipped with an i7-10510U CPU and 16GB RAM. Performance of the
V-LSPM is assessed in 4 ways: computing time, Procrustes correlation (PC), the
area under the receiver operating characteristic curve (AUROC), and the area
under the precision recall curve (AUPR). Time is measured in seconds with
shorter time indicating preferred performance. Procrustes correlation (Mardia
et al.,, 2008) measures the similarity between different latent space
configurations by minimising the sum of squared differences after Procrustes
rotations, with higher values (up to a maximum of 1) indicating better
correspondence. Here the PC is calculated between the estimated latent
positions and the true positions, conditioning on the smaller number of
dimensions between the two configurations compared. The AUROC and AUPR measure
the goodness-of-fit of the model by comparing the estimated edge probabilities
(calculated using the estimated model parameters and latent positions as in
(1)) with the observed data, with value of 1 implying perfect model fit. An
AUROC value of 0.5 implies random predictions, while the AUPR’s baseline value
depends on the density of the network. In networks with low density, the AUPR
will generally be lower than AUROC as it is based on the proportion of present
edges. Comparisons are made between the variational LSPM (V-LSPM) and the
MCMC-based LSPM (MCMC-LSPM), both implemented in the R package available in
the lspm repository.

The simulated data are generated by Bernoulli trials for each node pair based
on the probabilities derived from the distances between the nodes’ latent
positions. The latent positions are simulated according to (2) with the
shrinkage strengths being manually set to explore their effect across
different settings. Hyperparameters are set as
$\mu_{\alpha}=0,\sigma_{\alpha}=3,a_{1}=2,a_{2}=3,b_{1}=b_{2}=1$. A total of
30 networks are simulated in each case.

The simulation studies are structured as follows: Section 3.1 examines the
performance and impact of different initial truncation levels $p$ and Section
3.2 assesses the performance of the V-LSPM under different network sizes $n$.

###  3.1 Study 1: truncation level

Here the effect of different truncation levels of the latent dimension is
explored. Networks with $n=100$ are generated with the true number of latent
dimensions $p^{*}=4$ and shrinkage strengths of $\delta_{1}=0.5$,
$\delta_{2}=1.1$, $\delta_{3}=1.05$, and $\delta_{4}=1.15$. A value $\alpha=6$
is used, leading to networks with moderate density of $\simeq 20\%$. Three
initial truncation levels are considered: $p=\\{2,4,10\\}$, representing
situations where the initial truncation has been underestimated, correctly
specified, and overestimated, respectively.

Table 1 shows that the V-LSPM converges up to two orders of magnitude faster
than the MCMC-LSPM. Following Gwee et al., 2023a , 1,000,000 iterations were
run for the MCMC approach. Although the time per iteration is higher for
V-LSPM due to more complex computations, the substantial reduction in the
number of iterations results in lower computational cost. Table 1 also shows
good correspondence between estimated and true configurations, with PC values
between the true latent positions and the V-LSPM estimated latent positions
greater than 0.7. For both $p=4$ and $p=10$, the Procrustes correlation is
similar, indicating that a higher initial number of dimensions does not affect
the estimated latent positions configuration. The AUROC shows good model fit
with values of more than 0.8 across the different truncation levels but with
$p=\\{4,10\\}$ showing better and similar fit. Similarly, the AUPR shows good
model fit for $p=\\{4,10\\}$, since it takes values much higher than the
network density. These results suggest the use of $p>p^{*}$ is preferable to
$p<p^{*}$.

Table 1: Performance of the LSPM for different truncation levels $p$, assessed through average model run time, Procrustes correlation (PC), the area under the receiver operating characteristic curve (AUROC), and the area under the precision recall curve (AUPR). Standard deviations are given in brackets. $p$ | Time (sec) | PC | AUROC | AUPR
---|---|---|---|---
VI | MCMC
2 | 23 (2) | 2188 (57) | 0.71 (0.13) | 0.825 (0.027) | 0.480 (0.059)
4 | 57 (26) | 2381 (38) | 0.87 (0.08) | 0.918 (0.015) | 0.724 (0.057)
10 | 142 (104) | 2471 (89) | 0.87 (0.08) | 0.934 (0.006) | 0.786 (0.020)

Figure 1 shows that the variational posterior shrinkage strength estimates
have similar levels of accuracy regardless of the truncation level. When
$p>p^{*}=4$, Figure 1 shows that there is a large increase in the shrinkage
strength estimates on the 5th dimension, while, when $p\leq p^{*}$, there is
no sudden increase in shrinkage strength. Since a large shrinkage strength
means relatively smaller variance of the latent positions on that dimension,
this signifies the change from an effective dimension at $p=4$ to a non-
effective dimension at $p=5$ in the $p>p^{*}=4$ case. Although dimensions 7 to
10 tend to have relatively small shrinkage strength compared to dimension 5,
note that the precision of the latent positions is a cumulative product of the
shrinkage strengths from the previous dimensions. Therefore, the dimensions
above 5 will also have small variance and be non-effective. The absence of a
large increase in the posterior mean shrinkage strength across dimensions
indicates that either the truncation level employed is too low or is equal to
the true effective dimension. This behaviour suggests that using a reasonably
high truncation level is advisable in order to obtain a clear indication of
the number of effective dimensions. However, specifying $p$ to be too large
has a negative impact on computational efficiency.

Figure 1: Estimates of the variational posterior shrinkage strength across
different truncation levels $p=\\{2,4,10\\}$. Each dimension has 30 violin
plots, each representing the 10 point estimates from running V-LSPM with
different starting positions, for the 30 simulated networks. Red crosses are
the true parameter values used to simulate the networks.

Both Figure 1 and Figure 2(a) show differences between parameters estimates
and true values, particularly for $\delta_{1}$ and $\alpha$, indicating the
presence of bias regardless of the truncation level. In contrast, Figure 2(b)
shows that under the MCMC approach, except when $p<p^{*}$, there is little
bias when estimating $\alpha$. However, it should be noted that the lack of
bias observed in Figure 2(b) is confounded with the fact that while $p_{0}$
denotes the initial truncation level in the MCMC approach, an adaptive
procedure follows which typically leads to fewer inferred dimensions $p$ and a
better fitting model. In summary, while there is a loss in accuracy, this is
the trade-off for improved computational cost resulting from using an
approximation-based inferential method.

(a)

(b)

Figure 2: (a) V-LSPM estimates of $\alpha$ across different truncation levels
$p=\\{2,4,10\\}$. Each dimension has 30 violin plots, each representing the 10
point estimates from running V-LSPM with different starting positions, for the
30 simulated networks; (b) MCMC-LSPM posterior distributions of $\alpha$ for
different initial truncation levels $p_{0}$ and inferred values of $p$, across
30 simulated networks.

###  3.2 Study 2: network size

Here the focus is on assessing the effect of different network sizes on the
accuracy of the inference. Networks with $n=\\{20,50,100,200\\}$ are generated
with $p^{*}=2$, and shrinkage strengths of $\delta_{1}=0.5$ and
$\delta_{2}=1.1$ giving $\omega_{1}^{-1}=2$ and $\omega_{2}^{-1}=1.82$
respectively. Here, $\alpha=3$ is used meaning moderate network density of
$\simeq 20\%$. The V-LSPM with truncation level $p=5$ is fitted and inference
on the number of effective dimensions is achieved through inspection of the
variational posterior shrinkage strength parameters.

As observed in Section 3.1, Table 2 shows increased computational speed for
the V-LSPM approach compared to the MCMC-LSPM, where, as in Gwee et al., 2023a
, 500,000 iterations are used. Moreover, the V-LSPM demonstrates strong
performance in terms of accuracy as evidenced by the high PC, AUROC, and AUPR
values. As $n$ increases, Table 2 also shows increasingly higher PC, AUROC and
AUPR values.

Table 2: Performance of the LSPM across different network sizes $n$, assessed through average model run time, Procrustes correlation (PC), the area under the receiver operating characteristic curve (AUROC), and the area under the precision recall curve (AUPR). Standard deviations are given in brackets. $n$ | Time (sec) | PC | AUROC | AUPR
---|---|---|---|---
VI | MCMC
20 | 0.55 (0.40) | 304 (15) | 0.79 (0.12) | 0.918 (0.018) | 0.799 (0.058)
50 | 33 (20) | 481 (4) | 0.93 (0.04) | 0.905 (0.011) | 0.781 (0.035)
100 | 85 (24) | 1169 (8) | 0.95 (0.01) | 0.904 (0.007) | 0.788 (0.019)
200 | 437 (408) | 3920 (132) | 0.96 (0.01) | 0.905 (0.004) | 0.797 (0.011)

Figure 3 shows that, with increasing numbers of nodes, the variational
posterior shrinkage strength for dimension 3 becomes increasingly large in
value, signifying the increasingly clearer change from an effective dimension
at $p=2$ to a non-effective dimension at $p=3$.

Figure 3: Estimates of the variational posterior shrinkage strength across
networks of different sizes: $n=\\{20,50,100,200\\}$. Each dimension has 30
violin plots, each representing the 10 point estimates from running V-LSPM
with different starting positions, for the 30 simulated networks. Red crosses
are the true parameter values used to simulate the networks.

##  4 Illustrative network data examples

This section investigates the performance of V-LSPM on illustrative network
data examples. In addition to the variational LSPM (V-LSPM) and the MCMC-based
LSPM (MCMC-LSPM), also the variational-based LPM (V-LPM) and the MCMC-based
LPM (MCMC-LPM) are included in the comparison. The V-LPM is implemented
through the R package lvm4net (Gollini and Murphy,, 2016), while the MCMC-LPM
is implemented using the R package latentnet (Krivitsky and Handcock,, 2008,
2022).

###  4.1 Cat brain connectivity binary network data

The V-LSPM is used to analyse a cat brain connectivity network which
encompasses 65 distinct cortex regions, represented as nodes, and 1139
interregional macroscopic axonal projections, represented as edges (Scannell
et al.,, 1995; de Reus and van den Heuvel,, 2013). The regions within the
cortex are grouped into four primary categories: visual (18 areas), auditory
(10 areas), somatomotor (18 areas), and frontolimbic (19 areas). These
groupings stem from neurophysiological data detailing the specific function of
each brain region. Here, the network is viewed as a binary directed network. A
value of $y_{i,j}=1$ suggests a connection between regions $i$ and $j$, while
$y_{i,j}=0$ denotes the absence of such a connection. The network density is
$27.37\%$.

Fitting the MCMC-LSPM took an average of 10 ($\pm 1$) minutes (from 10
different starting values) via the MCMC approach with 500,000 iterations,
while the variational V-LSPM (from 10 different starting values) with $p=5$
took on average 16 ($\pm 4$) seconds and required 23 ($\pm 3$) iterations
(standard deviations are given in brackets). Figure 4(a) reports the estimated
shrinkage strength parameters and indicates 2 dimensions as the number of
effective dimensions; the MCMC approach suggests 3. Figure 4(a) also shows
that the inference on the number of dimensions is insensitive to the dimension
truncation level. The Procrustes correlation is assessed between the V-LSPM
latent positions and the 2-dimensional latent positions obtained from the
three other inferential methods: comparing with MCMC-LSPM gives PC = 0.69
($\pm 0.01$), with MCMC-LPM gives PC = 0.66 ($\pm 0.01$), and with V-LPM gives
PC = 0.67 ($\pm 0.01$). In terms of model fit, the AUROC and AUPR from the
V-LSPM are 0.891 ($\pm 0.006$) and 0.733 ($\pm 0.012$) respectively,
suggesting good model fit.

(a)

(b)

(c)

Figure 4: Cat brain connectivity data. (a) The V-LSPM estimates of the
variational posterior shrinkage strength parameters from multiple runs of
different starting values coloured according to different truncation levels,
(b) the inferred 2-dimensional variational posterior latent positions under
one of the V-LSPM runs, and (c) the Procrustes errors between the V-LSPM
variational posterior latent positions (arrow tail) and the MCMC-LSPM
posterior mean latent positions on the first two dimensions (arrow head).

Figure 4(b) shows the variational posterior latent positions under the V-LSPM.
Nodes from the same brain region are closely located with nodes from the
visuals and somatomotor regions lying below zero on the first dimension while
auditory and frontolimbic nodes mostly lie above zero. Similarly on dimension
2, nodes from the somatomotor region tend to lie below zero. Figure 4(c) shows
the Procrustes errors between the inferred V-LSPM variational posterior latent
positions and the MCMC-LSPM posterior mean latent positions on the first two
dimensions, indicating the overall structure of the configurations is very
similar with the V-LSPM positions being less dispersed than their MCMC
counterparts.

###  4.2 Worm nervous system binary network data

This binary directed network comprises $n=272$ nodes of neurons derived from
the nervous system of the adult male _Caenorhabditis elegans_ worm. Before
analysis, three nodes without edges were removed. Each of the nodes belongs to
one of 5 cell types: muscle (90 nodes), sensory (72 nodes), interneuron (67
nodes), motorneuron (38 nodes), and other (2 nodes). Each of the 4451 edges
signifies either a chemical or electrical interaction between a pair of nodes.
The data were compiled from serial electron micrograph sections, as documented
by Jarrell et al., (2012). The network density is 6.09%.

Model fitting took an average of 225 ($\pm 28$) minutes (from 10 different
starting values) via the MCMC approach with 1,000,000 iterations. On the other
hand, the V-LSPM (from 10 different starting values) took on average 8 ($\pm
1$) minutes with 43 ($\pm 2$) iterations (standard deviations given in
brackets). Figure 5 indicates 4 effective dimensions, while the MCMC approach
indicated 3. The Procrustes correlation is assessed between the latent
positions on the first 3-dimensions between the V-LSPM and the other methods
considered for comparison: with MCMC-LSPM, PC = 0.752 ($\pm 0.004$), with
MCMC-LPM, PC = 0.671 ($\pm 0.004$) and with the V-LPM, PC = 0.782 ($\pm
0.005$). Also, the AUROC and AUPR metrics for the V-LSPM are 0.853 ($\pm
0.002$) and 0.269 ($\pm 0.006$) respectively, suggesting good model fit. While
the AUPR value is low, it is higher than the baseline corresponding to the
network density, still indicating good model fit in terms of identification of
edges.

(a)

(b)

Figure 5: For the worm nervous system network, the V-LSPM point estimates of
the (a) variational posterior shrinkage strength parameters and (b)
variational posterior variance parameters.

Figure 6 shows the variational posterior latent positions’ estimates on the
first 4 dimensions. Nodes from the muscle cell type have the largest variance
across dimensions but tend to be negative on the first dimension. Nodes
belonging to the rest of the cell types are closely located to each other in
the latent space. On dimension 2, the sensory cells tend to be positive, while
the motorneuron cells are located below zero.

Figure 6: For the worm nervous system network, the latent position estimates
on the first 4 dimensions under the V-LSPM coloured by the cell types.

As shown by the PC values, while there is good correspondence between the
latent configurations under the MCMC and variational approaches, under the
V-LSPM, the latent positions struggle to move from the initial configuration,
initialised as detailed in Section 2.2.1, with a PC between the initial and
estimated configurations of 0.99 $\pm<0.01$. The MCMC approach instead manages
to move away from the initial configuration with PC between the initial and
posterior mean configurations of 0.56 $\pm 0.06$. Considering alternative
initialisation strategies for the V-LSPM may help here. Also, the V-LSPM does
not shrink the variance of the latent positions on the non-effective
dimensions as strongly as its MCMC counterpart. However, inference on the
number of effective dimensions was still clear in the case of the network
considered here. For reference, the initial positions and the MCMC-LSPM
posterior mean positions are given in Appendix D.

##  5 Discussion

In this work a variational inference approach for the LSPM has been developed
and has shown an improvement in computational cost when compared to the MCMC
approach. Although the variational approach is approximate, the V-LSPM
demonstrated good accuracy when inferring the number of effective latent
dimensions. Moreover, for the V-LSPM in the simulation studies, there was good
correspondence between the estimated and true latent positions, despite some
bias in some parameter estimates. Goodness of fit metrics via the AUROC and
AUPR also indicated satisfactory performance in both simulation studies and
applications.

Performance of the V-LSPM was observed to be poor in networks that were very
sparse or very dense, as in such cases $\alpha$ dominates the logistic
function in a way that the distances between the latent positions have little
influence on the edge formation process. This is especially problematic when
the truncation level is different from the true number of effective
dimensions, since biases are further incorporated into the estimate of
$\alpha$. The MCMC-LSPM addresses this problem via an adaptive sampling
procedure which automatically shrinks or grows the number of latent dimensions
as required, therefore reducing parameter bias. However, as the V-LSPM does
not sample from the true posterior distribution, a way to perform such an
adaptive sampling procedure is less intuitive. To obviate this problem, a
different shrinkage prior may be more suitable such as employing the Indian
buffet process (Knowles and Ghahramani,, 2011; Ročková and George,, 2016).

While the main focus of this work was to develop a computationally efficient
approach to inferring the LSPM, the proposed methodology could be easily
adapted to a range of network data types, much like its MCMC counterpart. For
example, developing a variational approach to inference for the LSPM for
networks with count edges as in Gwee et al., 2023a  or to perform clustering
of nodes using the LSPCM as in Gwee et al., 2023b  would be natural
extensions.

## Acknowledgments

This publication has emanated from research conducted with the financial
support of Science Foundation Ireland under grant number 18/CRT/6049. For the
purpose of Open Access, the author has applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission.

The authors are grateful for the discussions with Dr Isabella Gollini that
assisted with this work.

##  Appendix A Derivation of the KL divergence

The derivation of the KL divergence is derived with the help of the matrix
cookbook (Petersen and Pedersen,, 2012). By definition, the KL divergence
between the variational posterior $\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})$
and the true posterior
$\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})$ is:

$\displaystyle\mathrm{KL}[\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\|\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})]$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})}{\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\mathbf{Y},\alpha,\mathbf{Z},\bm{\delta})}{\mathbb{P}(\mathbf{Y})\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\mathbf{Y},\alpha,\mathbf{Z},\bm{\delta})}{\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})+\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{P}(\mathbf{Y})d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\mathbf{Y},\alpha,\mathbf{Z},\bm{\delta})}{\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})}d(\alpha,\mathbf{Z},\bm{\delta})+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\frac{\mathbb{P}(\mathbf{Y}\mid\alpha,\mathbf{Z})\mathbb{P}(\alpha)\prod_{i=1}^{n}\mathbb{P}(\mathbf{z}_{i})\mathbb{P}(\delta_{1})\prod_{h=2}^{p}\mathbb{P}(\delta_{h})}{\mathbb{Q}(\alpha)\prod_{i=1}^{n}\mathbb{Q}(\mathbf{z}_{i})\mathbb{Q}(\delta_{1})\prod_{h=2}^{p}\mathbb{Q}(\delta_{h})}d(\alpha,\mathbf{Z},\bm{\delta})+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle=-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{P}(\mathbf{Y}\mid\alpha,\mathbf{Z})d(\alpha,\mathbf{Z},\bm{\delta})+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle\qquad-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{P}(\alpha)d(\alpha,\mathbf{Z},\bm{\delta})+\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{Q}(\alpha)d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle\qquad-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\prod_{i=1}^{n}\mathbb{P}(\mathbf{z}_{i})d(\alpha,\mathbf{Z},\bm{\delta})+\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\prod_{i=1}^{n}\mathbb{Q}(\mathbf{z}_{i})d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle\qquad-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{P}(\delta_{1})d(\alpha,\mathbf{Z},\bm{\delta})+\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\mathbb{Q}(\delta_{1})d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle\qquad-\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\prod_{h=2}^{p}\mathbb{P}(\delta_{h})d(\alpha,\mathbf{Z},\bm{\delta})+\int\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\log\prod_{h=2}^{p}\mathbb{Q}(\delta_{h})d(\alpha,\mathbf{Z},\bm{\delta})$
$\displaystyle=-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}[\log(\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha))]+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle\qquad-\mathbb{E}_{\mathbb{Q}}[\log\mathbb{P}(\alpha)]+\mathbb{E}_{\mathbb{Q}}[\log\mathbb{Q}(\alpha)]-\sum_{i=1}^{n}\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{P}\left(\mathbf{z}_{i}\right)\right]+\sum_{i=1}^{n}\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{Q}\left(\mathbf{z}_{i}\right)\right]$
$\displaystyle\qquad-\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{P}\left(\delta_{1}\right)\right]+\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{Q}\left(\delta_{1}\right)\right]-\sum_{h=2}^{p}\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{P}\left(\delta_{h}\right)\right]+\sum_{h=2}^{p}\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{Q}\left(\delta_{h}\right)\right]$
$\displaystyle=-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log(\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)\right]+\log\mathbb{P}(\mathbf{Y})+\operatorname{KL}\left[\mathbb{Q}(\alpha)\|\mathbb{P}(\alpha)\right]$
$\displaystyle\qquad+\sum_{i=1}^{n}\operatorname{KL}\left[\mathbb{Q}\left(\mathbf{z}_{i}\right)\|\mathbb{P}\left(\mathbf{z}_{i}\right)\right]+\operatorname{KL}\left[\mathbb{Q}(\delta_{1})\|\mathbb{P}(\delta_{1})\right]+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]$

Here follows the derivation of
$\operatorname{KL}\left[\mathbb{Q}(\alpha)\|\mathbb{P}(\alpha)\right]$:

$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}[\log{\mathbb{Q}}(\alpha)]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{(\alpha-\tilde{\mu}_{\alpha})^{2}}{2\tilde{\sigma}_{\alpha}^{2}}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{(\alpha^{2}-2\alpha\tilde{\mu}_{\alpha}+\tilde{\mu}_{\alpha}^{2})}{2\tilde{\sigma}_{\alpha}^{2}}\right]$
$\displaystyle=-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\alpha^{2}\right]-2\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\alpha\right]\tilde{\mu}_{\alpha}+\tilde{\mu}_{\alpha}^{2}}{2\tilde{\sigma}_{\alpha}^{2}}$
$\displaystyle=-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{\tilde{\mu}_{\alpha}^{2}+\tilde{\sigma}_{\alpha}^{2}-2\tilde{\mu}_{\alpha}\tilde{\mu}_{\alpha}+\tilde{\mu}_{\alpha}^{2}}{2\tilde{\sigma}_{\alpha}^{2}}$
$\displaystyle=-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{1}{2}$
$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}[\log{\mathbb{P}}(\alpha)]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[-\frac{1}{2}\log\left(2\pi\sigma_{\alpha}^{2}\right)-\frac{(\alpha-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[-\frac{1}{2}\log\left(2\pi\sigma_{\alpha}^{2}\right)-\frac{(\alpha^{2}-2\alpha\mu_{\alpha}+\mu_{\alpha}^{2})}{2\sigma_{\alpha}^{2}}\right]$
$\displaystyle=-\frac{1}{2}\log\left(2\pi\sigma_{\alpha}^{2}\right)-\frac{\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\alpha^{2}\right]-2\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\alpha\right]\mu_{\alpha}+\mu_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle=-\frac{1}{2}\log\left(2\pi\sigma_{\alpha}^{2}\right)-\frac{\tilde{\mu}_{\alpha}^{2}+\tilde{\sigma}_{\alpha}^{2}-2\tilde{\mu}_{\alpha}\mu_{\alpha}+\mu_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle\operatorname{KL}\left[\mathbb{Q}(\alpha)\|\mathbb{P}(\alpha)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\log{\mathbb{Q}}(\alpha)\right]-\mathbb{E}_{\mathbb{Q}}\left[\log{\mathbb{P}}(\alpha)\right]$
$\displaystyle=\left\\{-\frac{1}{2}\log\left(2\pi\tilde{\sigma}_{\alpha}^{2}\right)-\frac{1}{2}\right\\}-\left\\{-\frac{1}{2}\log\left(2\pi\sigma_{\alpha}^{2}\right)-\frac{\tilde{\mu}_{\alpha}^{2}+\tilde{\sigma}_{\alpha}^{2}-2\tilde{\mu}_{\alpha}\mu_{\alpha}+\mu_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}\right\\}$
$\displaystyle=-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}+(\tilde{\mu}_{\alpha}^{2}-2\tilde{\mu}_{\alpha}\mu_{\alpha}+\mu_{\alpha}^{2})}{2\sigma_{\alpha}^{2}}$
$\displaystyle=-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}+(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle=-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}+\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}$

Here follows the derivation of
$\operatorname{KL}\left[\mathbb{Q}\left(\mathbf{z}_{i}\right)\|\mathbb{P}\left(\mathbf{z}_{i}\right)\right]$:

$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{Q}}\left(\mathbf{z}_{i}\mid\mathbf{\Omega}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log\left\\{\left(\frac{1}{2\pi}\right)^{p/2}\operatorname{det}(\tilde{\mathbf{\Omega}})^{1/2}\exp\left[-\frac{1}{2}\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)^{\top}\tilde{\mathbf{\Omega}}\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)\right]\right\\}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{1}{2}\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)^{\top}\tilde{\mathbf{\Omega}}\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)\right]$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{1}{2}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)^{\top}\tilde{\mathbf{\Omega}}\left(\mathbf{z}_{i}-\tilde{\mathbf{z}}_{i}\right)\right]$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{1}{2}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\mathbf{z}_{i}^{\top}\tilde{\mathbf{\Omega}}\mathbf{z}_{i}-\mathbf{z}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\mathbf{z}_{i}+\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}\right]$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})$
$\displaystyle\qquad-\frac{1}{2}\left\\{\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\mathbf{z}_{i}^{\top}\tilde{\mathbf{\Omega}}\mathbf{z}_{i}\right]-\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\mathbf{z}_{i}^{\top}\right]\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\mathbf{z}_{i}\right]+\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}\right\\}$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})$
$\displaystyle\qquad-\frac{1}{2}\left\\{\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}+\operatorname{tr}(\tilde{\mathbf{\Omega}}\tilde{\mathbf{\Omega}}^{-1})-\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}+\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{\Omega}}\tilde{\mathbf{z}}_{i}\right\\}$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{1}{2}\operatorname{tr}(\tilde{\mathbf{\Omega}}\tilde{\mathbf{\Omega}}^{-1})$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{p}{2}$
$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{P}}\left(\mathbf{z}_{i}\mid\mathbf{\Omega}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log\left\\{\left(\frac{1}{2\pi}\right)^{p/2}\operatorname{det}(\mathbf{\Omega})^{1/2}\exp\left[-\frac{1}{2}\left(\mathbf{z}_{i}-0\right)^{\top}\mathbf{\Omega}\left(\mathbf{z}_{i}-0\right)\right]\right\\}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\sum_{i=1}^{n}\left\\{-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega})-\frac{1}{2}\mathbf{z}_{i}^{\top}\mathbf{\Omega}\mathbf{z}_{i}\right\\}\right]$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega})-\frac{1}{2}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\mathbf{z}_{i}^{\top}\mathbf{\Omega}\mathbf{z}_{i}\right]$
$\displaystyle=-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega})-\frac{1}{2}\left[\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})\right]$

$\displaystyle\operatorname{KL}\left[\mathbb{Q}\left(\mathbf{z}_{i}\right)\|\mathbb{P}\left(\mathbf{z}_{i}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{Q}}\left(\mathbf{z}_{i}\mid\mathbf{\Omega}\right)\right]-\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{P}}\left(\mathbf{z}_{i}\mid\mathbf{\Omega}\right)\right]$
$\displaystyle=\left\\{-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\tilde{\mathbf{\Omega}})-\frac{p}{2}\right\\}$
$\displaystyle\qquad-\left\\{-\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega})-\frac{1}{2}\left[\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})\right]\right\\}$
$\displaystyle=-\frac{p}{2}+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega}^{-1}\tilde{\mathbf{\Omega}})+\frac{1}{2}\left[\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})\right]$
$\displaystyle=-\frac{p}{2}+\frac{1}{2}\log\operatorname{det}(\mathbf{\Omega}^{-1}\tilde{\mathbf{\Omega}})+\frac{1}{2}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\frac{1}{2}\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})$

Here follows the derivation of
$\operatorname{KL}\left[\mathbb{Q}(\delta_{1})\|\mathbb{P}(\delta_{1})\right]$:

$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{Q}}\left(\delta_{1}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log\left\\{\frac{\tilde{b}_{1}^{\tilde{a}_{1}}}{\Gamma\left(\tilde{a}_{1}\right)}\delta_{1}^{\tilde{a}_{1}-1}\exp\left(-\tilde{b}_{1}\delta_{1}\right)\right\\}\right]$
$\displaystyle=\tilde{a}_{1}\log\tilde{b}_{1}-\log\Gamma\left(\tilde{a}_{1}\right)+\left(\tilde{a}_{1}-1\right)\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log\left(\delta_{1}\right)\right]-\tilde{b}_{1}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\delta_{1}\right]$
$\displaystyle=\tilde{a}_{1}\log\tilde{b}_{1}-\log\Gamma\left(\tilde{a}_{1}\right)+\left(\tilde{a}_{1}-1\right)\left[\psi(\tilde{a}_{1})+\log(\tilde{b}_{1})\right]-\tilde{b}_{1}\left(\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right)$
$\displaystyle=\tilde{a}_{1}\log\tilde{b}_{1}-\log\Gamma\left(\tilde{a}_{1}\right)+\left(\tilde{a}_{1}-1\right)\left[\psi(\tilde{a}_{1})+\log(\tilde{b}_{1})\right]-\tilde{a}_{1}$
$\displaystyle\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{P}}\left(\delta_{1}\right)\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log\left\\{\frac{b_{1}^{a_{1}}}{\Gamma\left(a_{1}\right)}\delta_{1}^{a_{1}-1}\exp\left(-b_{1}\delta_{1}\right)\right\\}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[a_{1}\log
b_{1}-\log\Gamma\left(a_{1}\right)+\left(a_{1}-1\right)\log\left(\delta_{1}\right)-b_{1}\delta_{1}\right]$
$\displaystyle=a_{1}\log
b_{1}-\log\Gamma\left(a_{1}\right)+\left(a_{1}-1\right)\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log(\delta_{1})\right]-b_{1}\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\delta_{1}\right]$
$\displaystyle=a_{1}\log
b_{1}-\log\Gamma\left(a_{1}\right)+\left(a_{1}-1\right)\left[\psi(\tilde{a}_{1})+\log(\tilde{b}_{1})\right]-b_{1}\left(\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right)$
$\displaystyle\operatorname{KL}\left[\mathbb{Q}(\delta_{1})\|\mathbb{P}(\delta_{1})\right]$
$\displaystyle=\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{Q}}\left(\delta_{1}\right)\right]-\mathbb{E}_{\mathbb{{\mathbb{Q}}}}\left[\log{\mathbb{P}}\left(\delta_{1}\right)\right]$
$\displaystyle=\left\\{\tilde{a}_{1}\log\tilde{b}_{1}-\log\Gamma\left(\tilde{a}_{1}\right)+\left(\tilde{a}_{1}-1\right)\left[\psi(\tilde{a}_{1})+\log(\tilde{b}_{1})\right]-\tilde{a}_{1}\right\\}$
$\displaystyle\quad-\left\\{a_{1}\log
b_{1}-\log\Gamma\left(a_{1}\right)+\left(a_{1}-1\right)\left[\psi(\tilde{a}_{1})+\log(\tilde{b}_{1})\right]-b_{1}\left(\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right)\right\\}$
$\displaystyle=\tilde{a}_{1}\left[\psi(\tilde{a}_{1})+\frac{b_{1}}{\tilde{b}_{1}}+1\right]-a_{1}\left[\psi(\tilde{a}_{1})+\log(\frac{b_{1}}{\tilde{b}_{1}})\right]-\log\frac{\Gamma(\tilde{a}_{1})}{\Gamma(a_{1})}$

Here follows the derivation of
$\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log(\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)\right]$:

$\displaystyle\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}[\log(\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha))]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log\left\\{\prod_{i\neq
j}^{n}\frac{\exp\left[y_{i,j}(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2})\right]}{1+\exp(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2})}\right\\}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\sum_{i\neq
j}^{n}\left\\{y_{i,j}\alpha-
y_{i,j}\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right.\right.$
$\displaystyle\quad\left.\left.-\log\left[1+\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right]\right\\}\right]$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right]$
$\displaystyle\quad-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log\left(1+\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right)\right]$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\mathbb{E}_{\mathbb{Q}}\left[\alpha\right]-\mathbb{E}_{\mathbb{Q}}\left[\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right]$
$\displaystyle\quad-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log\left(1+\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right)\right]$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log\left(1+\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right)\right]$
since log is concave, by Jensen inequality: $\displaystyle\leq\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad-\log\left(1+\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right]\right)$
$\displaystyle=\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad-\log\left\\{1+\frac{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\right.$
$\displaystyle\qquad\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}.$

$\displaystyle\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\exp\left(\alpha-\left\|\mathbf{z}_{i}-\mathbf{z}_{j}\right\|^{2}\right)\right]=$
$\displaystyle\mathbb{E}_{\mathbb{Q}}\left[\exp(\alpha)\right]\mathbb{E}_{\mathbb{Q}}\left[\exp(-\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\|)^{2}\right]$
$\displaystyle=$
$\displaystyle\left[\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)\right]\mathbb{E}_{\mathbb{Q}}\left[\exp(-\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\|)^{2}\right]$
$\displaystyle=$
$\displaystyle\left[\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)\right]\left\\{\operatorname{det}\left(4\tilde{\mathbf{\Omega}}^{-1}+\mathbf{I}\right)^{-1/2}\right.$
$\displaystyle\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(4\tilde{\mathbf{\Omega}}^{-1}+\mathbf{I}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}$

$\displaystyle\mathbb{Q}\left(\mathbf{z}_{i}\right)\sim\operatorname{MVN}$
$\displaystyle\left(\tilde{\mathbf{z}}_{i},\tilde{\mathbf{\Omega}}^{-1}\right)\qquad\mathbb{Q}\left(\mathbf{z}_{j}\right)\sim\operatorname{MVN}\left(\tilde{\mathbf{z}}_{j},\tilde{\mathbf{\Omega}}^{-1}\right)$
$\displaystyle\therefore\mathbb{Q}\left(\mathbf{z}_{i}-\mathbf{z}_{j}\right)$
$\displaystyle\sim\operatorname{MVN}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j},2\tilde{\mathbf{\Omega}}^{-1}\right)$
$\displaystyle\text{ let
}\mathbf{D}_{i,j}=\mathbf{z}_{i}-\mathbf{z}_{j},\tilde{\mathbf{S}}$
$\displaystyle=2\tilde{\mathbf{\Omega}}^{-1}\qquad\therefore\mathbb{Q}\left(\mathbf{D}_{i,j}\right)\sim\operatorname{MVN}\left(\tilde{\mathbf{D}}_{i,j},\tilde{\mathbf{S}}\right)$

$\displaystyle\mathbb{E}_{\mathbb{Q}}\left[\exp\left(-\mathbf{D}_{i,j}^{\top}\mathbf{D}_{i,j}\right)\right]=$
$\displaystyle\int_{-\infty}^{\infty}\exp\left(-\mathbf{D}_{i,j}^{\top}\mathbf{D}_{i,j}\right)\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}$
$\displaystyle\exp\left[-\frac{1}{2}\left(\mathbf{D}_{i,j}-\tilde{\mathbf{D}}_{i,j}\right)^{\top}\tilde{\mathbf{S}}^{-1}\left(\mathbf{D}_{i,j}-\tilde{\mathbf{D}}_{i,j}\right)\right]d\mathbf{D}_{i,j}$
$\displaystyle=$
$\displaystyle\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}\int_{-\infty}^{\infty}\exp\left[-\mathbf{D}_{i,j}^{\top}\mathbf{D}_{i,j}-\frac{1}{2}\mathbf{D}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\mathbf{D}_{i,j}\right.$
$\displaystyle\left.+\mathbf{D}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]d\mathbf{D}_{i,j}$
$\displaystyle=$
$\displaystyle\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}\exp\left[-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]$
$\displaystyle\int_{-\infty}^{\infty}\exp\left[-\mathbf{D}_{i,j}^{\top}\mathbf{D}_{i,j}-\frac{1}{2}\mathbf{D}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\mathbf{D}_{i,j}+\mathbf{D}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]d\mathbf{D}_{i,j}$
$\displaystyle=$
$\displaystyle\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}\exp\left[-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]$
$\displaystyle\int_{-\infty}^{\infty}\exp\left\\{-\frac{1}{2}\mathbf{D}_{i,j}^{\top}\left[\left(\frac{1}{2}\mathbf{I}\right)^{-1}+\tilde{\mathbf{S}}^{-1}\right]\mathbf{D}_{i,j}+\mathbf{D}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right\\}d\mathbf{D}_{i,j}$

$\displaystyle\mathbb{E}_{\mathbb{Q}}$
$\displaystyle\left[\exp\left(-\mathbf{D}_{i,j}^{\top}\mathbf{D}_{i,j}\right)\right]$
$\displaystyle=$
$\displaystyle\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}\exp\left[-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]\operatorname{det}\left\\{2\pi\left[\left(\frac{1}{2}\mathbf{I}\right)^{-1}+\tilde{\mathbf{S}}^{-1}\right]^{-1}\right\\}^{1/2}$
$\displaystyle\exp\left\\{\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\left[\left(\frac{1}{2}\mathbf{I}\right)^{-1}+\tilde{\mathbf{S}}^{-1}\right]^{-1}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle=$
$\displaystyle\operatorname{det}(2\pi\tilde{\mathbf{S}})^{-1/2}\operatorname{det}\left\\{2\pi\left[\left(\frac{1}{2}\mathbf{I}\right)^{-1}+\tilde{\mathbf{S}}^{-1}\right]^{-1}\right\\}^{1/2}$
$\displaystyle\exp\left[-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right]\exp\left\\{\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\left[\left(\frac{1}{2}\mathbf{I}\right)^{-1}+\tilde{\mathbf{S}}^{-1}\right]^{-1}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle=$
$\displaystyle\operatorname{det}\left\\{\tilde{\mathbf{S}}\left[2\mathbf{I}+\tilde{\mathbf{S}}^{-1}\right]\right\\}^{-1/2}\exp\left\\{-\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{D}}_{i,j}\right.$
$\displaystyle\left.+\frac{1}{2}\tilde{\mathbf{D}}_{i,j}^{\top}\left[2\tilde{\mathbf{S}}\tilde{\mathbf{S}}+\tilde{\mathbf{S}}\tilde{\mathbf{S}}^{-1}\tilde{\mathbf{S}}\right]^{-1}\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle=$
$\displaystyle\operatorname{det}\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1/2}\exp\left\\{\tilde{\mathbf{D}}_{i,j}^{\top}\left(\frac{1}{2}\tilde{\mathbf{S}}^{-1}\right)\left[-\mathbf{I}+\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1}\right]\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle=$
$\displaystyle\operatorname{det}\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1/2}\exp\left\\{\tilde{\mathbf{D}}_{i,j}^{\top}\left(\frac{1}{2}\tilde{\mathbf{S}}^{-1}\right)\left[\left(-\mathbf{I}-2\tilde{\mathbf{S}}+\mathbf{I}\right)\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1}\right]\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle=$
$\displaystyle\operatorname{det}\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1/2}\exp\left\\{-\tilde{\mathbf{D}}_{i,j}^{\top}\left(2\tilde{\mathbf{S}}+\mathbf{I}\right)^{-1}\tilde{\mathbf{D}}_{i,j}\right\\}$
$\displaystyle\therefore\mathbb{E}_{\mathbb{Q}}$
$\displaystyle\left[\exp(-\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\|)^{2}\right]$
$\displaystyle=$
$\displaystyle\operatorname{det}\left(4\tilde{\mathbf{\Omega}}^{-1}+\mathbf{I}\right)^{-1/2}\exp\left\\{-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(4\tilde{\mathbf{\Omega}}^{-1}+\mathbf{I}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right\\}$

To sum up, the KL divergence is:

$\displaystyle\mathrm{KL}$
$\displaystyle[\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\|\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})]$
$\displaystyle=$
$\displaystyle\operatorname{KL}\left[\mathbb{Q}(\alpha)\|\mathbb{P}(\alpha)\right]+\sum_{i=1}^{n}\operatorname{KL}\left[\mathbb{Q}\left(\mathbf{z}_{i}\right)\|\mathbb{P}\left(\mathbf{z}_{i}\right)\right]+\operatorname{KL}\left[\mathbb{Q}(\delta_{1})\|\mathbb{P}(\delta_{1})\right]$
$\displaystyle+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]-\mathbb{E}_{\mathbb{Q}(\mathbf{Z},\alpha\mid\mathbf{Y})}\left[\log(\mathbb{P}(\mathbf{Y}\mid\mathbf{Z},\alpha)\right]+\log\mathbb{P}(\mathbf{Y})$
$\displaystyle\leq$
$\displaystyle-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}+\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle-\frac{np}{2}+\frac{n}{2}\log\operatorname{det}(\mathbf{\Omega}^{-1}\tilde{\mathbf{\Omega}})+\frac{1}{2}\sum_{i=1}^{n}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\frac{n}{2}\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})$
$\displaystyle+\tilde{a}_{1}\left[\psi(\tilde{a}_{1})+\frac{b_{1}}{\tilde{b}_{1}}+1\right]-a_{1}\left[\psi(\tilde{a}_{1})+\log(\frac{b_{1}}{\tilde{b}_{1}})\right]-\log\frac{\Gamma(\tilde{a}_{1})}{\Gamma(a_{1})}$
$\displaystyle+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]$
$\displaystyle-\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad+\log\left\\{1+\frac{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\right.$
$\displaystyle\qquad\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}$
$\displaystyle+\log\mathbb{P}(\mathbf{Y})$

##  Appendix B Summary for updating the variational parameters

The variational parameters are updated in turns, while holding the others
constant as follow:

  1. 1.

Optimize $\tilde{\mu}_{\alpha}$ via bracketing-and-bisection algorithm on the
following expression:

$\displaystyle\frac{\partial\operatorname{KL}}{\partial\tilde{\mu}_{\alpha}}$
$\displaystyle=\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})}{\sigma_{\alpha}^{2}}-\sum_{i\neq
j}^{n}y_{i,j}+\sum_{i\neq
j}^{n}\left[\exp\left(-\tilde{\mu}_{\alpha}\right)A_{i,j}^{-1}+1\right]^{-1}$
(5) $\displaystyle\text{where }A_{i,j}$
$\displaystyle=\exp\left(\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)\operatorname{det}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1/2}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]$

  2. 2.

Optimize $\tilde{\sigma}^{2}_{\alpha}$ via bracketing-and-bisection algorithm
on the following expression:

$\displaystyle\frac{\partial\operatorname{KL}}{\partial\tilde{\sigma}_{\alpha}^{2}}$
$\displaystyle=\frac{1}{2\sigma_{\alpha}^{2}}-\frac{1}{2\tilde{\sigma}_{\alpha}^{2}}+\sum_{i\neq
j}^{n}\frac{1}{2}\left[1+\exp(-\frac{1}{2}\tilde{\sigma}_{\alpha}^{2})B_{i,j}^{-1}\right]^{-1}$
(6) $\displaystyle\text{where }B_{i,j}$
$\displaystyle=\exp\left(\tilde{\mu}_{\alpha}\right)\operatorname{det}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1/2}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]$

  3. 3.

Optimize $\tilde{\mathbf{z}}_{i}$ via conjugate gradient routine on the
following expression:

$\displaystyle\frac{\partial\operatorname{KL}}{\partial\tilde{\mathbf{z}}_{i}}$
$\displaystyle=2\tilde{\mathbf{z}}_{i}\left(\frac{1}{2}\bm{\Omega}+\sum_{i\neq
j}^{n}y_{i,j}\right)-2\left(\sum_{i\neq
j}^{n}y_{i,j}\tilde{\mathbf{z}}_{j}\right)+2G(\tilde{\mathbf{z}}_{i})$ (7)
$\displaystyle\text{where }G\left(\tilde{\mathbf{z}}_{i}\right)$
$\displaystyle=-2\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)\sum_{i\neq
j}^{n}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\left\\{1+\frac{\operatorname{det}\left(\mathbf{I}+\tilde{\bm{\Omega}}^{-1}\right)^{1/2}}{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}\right.$
$\displaystyle\quad\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}^{-1}$

  4. 4.

Updates $\delta_{1}$ analytically with the following hyperparameters:

$\tilde{a}_{1}=\frac{np}{2}+a_{1}\qquad\tilde{b}_{1}=\quad
b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\mathbb{E}(\delta_{m})\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]$

  5. 5.

Updates $\delta_{h}$ for $h=2,\ldots,p$ analytically with the following
hyperparameters:

$\tilde{a}_{2}^{h}=\frac{n(p-h+1)}{2}+a_{2}\qquad\tilde{b}_{2}^{h}=\quad
b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\prod_{m=1,m\neq
h}^{\ell}\mathbb{E}(\delta_{m})\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]$

##  Appendix C Derivations of the optimal variational parameters

The full approximate KL divergence needed for the starting point of the
partial derivatives is as follow:

$\displaystyle\mathrm{KL}\left[\mathbb{Q}(\alpha,\mathbf{Z},\bm{\delta})\|\mathbb{P}(\alpha,\mathbf{Z},\bm{\delta}\mid\mathbf{Y})\right]\leq$
$\displaystyle-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)-\frac{1}{2}+\frac{\tilde{\sigma}_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}+\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}$
$\displaystyle-\frac{np}{2}+\frac{n}{2}\log\operatorname{det}(\mathbf{\Omega}^{-1}\tilde{\mathbf{\Omega}})+\frac{1}{2}\sum_{i=1}^{n}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\frac{n}{2}\operatorname{tr}(\mathbf{\Omega}\tilde{\mathbf{\Omega}}^{-1})$
$\displaystyle+\tilde{a}_{1}\left[\psi(\tilde{a}_{1})+\frac{b_{1}}{\tilde{b}_{1}}+1\right]-a_{1}\left[\psi(\tilde{a}_{1})+\log(\frac{b_{1}}{\tilde{b}_{1}})\right]-\log\frac{\Gamma(\tilde{a}_{1})}{\Gamma(a_{1})}$
$\displaystyle+\sum_{h=2}^{p}\operatorname{KL}\left[\mathbb{Q}\left(\delta_{h}\right)\|\mathbb{P}\left(\delta_{h}\right)\right]$
$\displaystyle-\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mu}_{\alpha}-2\operatorname{tr}(\tilde{\mathbf{\Omega}}^{-1})-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)$
$\displaystyle\quad+\log\left\\{1+\frac{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\right.$
$\displaystyle\qquad\left.\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}+\log\mathbb{P}(\mathbf{Y})$

The partial derivatives for the variational parameter $\tilde{\mu}_{\alpha}$:

$\displaystyle\mathrm{KL}_{\tilde{\mu}_{\alpha}}\leq$
$\displaystyle\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}-\sum_{i\neq
j}^{n}y_{i,j}\tilde{\mu}_{\alpha}+\sum_{i\neq
j}^{n}\log\left[1+\exp\left(\tilde{\mu}_{\alpha}\right)A_{i,j}\right]$
$\displaystyle\text{where }A_{i,j}=$
$\displaystyle\frac{\exp\left(\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]$
$\displaystyle\frac{\partial\mathrm{KL}_{\tilde{\mu}_{\alpha}}}{\partial\tilde{\mu}_{\alpha}}=$
$\displaystyle\frac{2(\tilde{\mu}_{\alpha}-\mu_{\alpha})}{2\sigma_{\alpha}^{2}}-\sum_{i\neq
j}^{n}y_{i,j}+\sum_{i\neq
j}^{n}\left[\frac{\exp\left(\tilde{\mu}_{\alpha}\right)A_{i,j}}{1+\exp\left(\tilde{\mu}_{\alpha}\right)A_{i,j}}\right]$
$\displaystyle=$
$\displaystyle\frac{(\tilde{\mu}_{\alpha}-\mu_{\alpha})}{\sigma_{\alpha}^{2}}-\sum_{i\neq
j}^{n}y_{i,j}+\sum_{i\neq
j}^{n}\left[\exp\left(-\tilde{\mu}_{\alpha}\right)A_{i,j}^{-1}+1\right]^{-1}$

The partial derivatives for the variational parameter
$\tilde{\sigma}_{\alpha}^{2}$:

$\displaystyle\mathrm{KL}_{\tilde{\sigma}_{\alpha}^{2}}\leq$
$\displaystyle-\frac{1}{2}\log\left(\frac{\tilde{\sigma}_{\alpha}^{2}}{\sigma_{\alpha}^{2}}\right)+\frac{\tilde{\sigma}_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}+\sum_{i\neq
j}^{n}\log\left[1+\exp\left(\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)B_{i,j}\right]$
$\displaystyle\text{where }B_{i,j}=$
$\displaystyle\frac{\exp\left(\tilde{\mu}_{\alpha}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]$
$\displaystyle\frac{\partial\mathrm{KL}_{\tilde{\sigma}_{\alpha}^{2}}}{\partial\tilde{\sigma}_{\alpha}^{2}}=$
$\displaystyle-\frac{1}{2\tilde{\sigma}_{\alpha}^{2}}+\frac{1}{2\sigma_{\alpha}^{2}}+\sum_{i\neq
j}^{n}\log\left[\exp\left(-\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)B_{i,j}^{-1}+1\right]^{-1}$

The partial derivatives for the variational parameter
$\tilde{\mathbf{z}}_{i}$:

$\displaystyle\mathrm{KL}_{\tilde{\mathbf{z}}_{i}}\leq$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}-\sum_{i\neq
j}^{n}y_{i,j}\left(-\left\|\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right\|^{2}\right)+G\left(\tilde{\mathbf{z}}_{i}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\tilde{\mathbf{z}}_{i}^{\top}\mathbf{\Omega}\tilde{\mathbf{z}}_{i}+\sum_{i\neq
j}^{n}y_{i,j}\left(\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{z}}_{i}-2\tilde{\mathbf{z}}_{i}^{\top}\tilde{\mathbf{z}}_{j}+\tilde{\mathbf{z}}_{j}^{\top}\tilde{\mathbf{z}}_{j}\right)+G\left(\tilde{\mathbf{z}}_{i}\right)$
$\displaystyle\text{where }G\left(\tilde{\mathbf{z}}_{i}\right)=$
$\displaystyle\log\left\\{1+\frac{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}{\operatorname{det}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{\frac{1}{2}}}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}(\mathbf{I}+4\tilde{\mathbf{\Omega}}^{-1})^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}$
$\displaystyle\text{note that if }F(x)=$
$\displaystyle\log\left[1+a\exp(bX)\right],\text{then: }$ $\displaystyle
F^{\prime}(x)=$
$\displaystyle\left[1+a\exp(bX)\right]^{-1}\left[a\exp(bX)\right]\left[b\right]=\left[1+a^{-1}\exp(-bX)\right]^{-1}\left[b\right]$
$\displaystyle\therefore G^{\prime}\left(\tilde{\mathbf{z}}_{i}\right)=$
$\displaystyle\left\\{1+\frac{\operatorname{det}\left(\mathbf{I}+\tilde{\bm{\Omega}}^{-1}\right)^{1/2}}{\exp\left(\tilde{\mu}_{\alpha}+\frac{1}{2}\tilde{\sigma}_{\alpha}^{2}\right)}\exp\left[-\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)^{\top}\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right]\right\\}^{-1}$
$\displaystyle\left\\{-2\left(\mathbf{I}+4\tilde{\bm{\Omega}}^{-1}\right)^{-1}\sum_{i\neq
j}^{n}\left(\tilde{\mathbf{z}}_{i}-\tilde{\mathbf{z}}_{j}\right)\right\\}$
$\displaystyle\frac{\partial\mathrm{KL}_{\tilde{\mathbf{z}}_{i}}}{\partial\tilde{\mathbf{z}}_{i}}=$
$\displaystyle 2\tilde{\mathbf{z}}_{i}\left(\frac{1}{2}\bm{\Omega}+\sum_{i\neq
j}^{n}y_{i,j}\right)-2\left(\sum_{i\neq
j}^{n}y_{i,j}\tilde{\mathbf{z}}_{j}\right)+G^{\prime}(\tilde{\mathbf{z}}_{i})$

Identifying the gamma density of $\delta_{1}$ to update for the variational
parameters $\tilde{a}_{1}$ and $\tilde{b}_{1}$

$\displaystyle\mathbb{P}(\delta_{1}|-)\propto$
$\displaystyle\delta_{1}^{\frac{np}{2}+a_{1}-1}\exp\left[-\left(b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\delta_{m}z_{i\ell}^{2}\right)\delta_{1}\right]$
$\displaystyle\log\mathbb{P}(\delta_{1}|-)\propto$
$\displaystyle\left(\frac{np}{2}+a_{1}-1\right)\log\left(\delta_{1}\right)-\left(b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\delta_{m}z_{i\ell}^{2}\right)\delta_{1}$
$\displaystyle\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{P}(\delta_{1}|-)\right]\propto$
$\displaystyle\left(\frac{np}{2}+a_{1}-1\right)\log\left(\delta_{1}\right)-\left(b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\mathbb{E}_{\mathbb{Q}}\left[\delta_{m}\right]\mathbb{E}_{\mathbb{Q}}\left[z_{i\ell}^{2}\right]\right)\delta_{1}$
$\displaystyle=$
$\displaystyle\left(\frac{np}{2}+a_{1}-1\right)\log\left(\delta_{1}\right)-\left(b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\left[\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\right]\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]\right)\delta_{1}$
$\displaystyle\delta_{1}$
$\displaystyle\sim\text{Gam}\left(\frac{np}{2}+a_{1},\quad
b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]\right)$
$\displaystyle\therefore\tilde{a}_{1}$
$\displaystyle=\frac{np}{2}+a_{1},\quad\tilde{b}_{1}=b_{1}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=1}^{p}\prod_{m=2}^{\ell}\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]$

Identifying the truncated gamma density of $\delta_{h}$ to update for the
variational parameters $\tilde{a}_{1}$ and $\tilde{b}_{1}$

$\displaystyle\mathbb{P}(\delta_{h}|-)\propto$
$\displaystyle\delta_{h}^{\frac{n(p-h+1)}{2}+a_{2}-1}\exp\left[-\left(b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\prod_{m=1,m\neq
h}^{\ell}\delta_{m}z_{i\ell}^{2}\right)\delta_{h}\right]$
$\displaystyle\log\mathbb{P}(\delta_{h}|-)\propto$
$\displaystyle\left(\frac{n(p-h+1)}{2}+a_{2}-1\right)\log\left(\delta_{h}\right)-\left(b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\prod_{m=1,m\neq
h}^{\ell}\delta_{m}z_{i\ell}^{2}\right)\delta_{h}$
$\displaystyle\mathbb{E}_{\mathbb{Q}}\left[\log\mathbb{P}(\delta_{h}|-)\right]\propto$
$\displaystyle\left(\frac{n(p-h+1)}{2}+a_{2}-1\right)\log\left(\delta_{h}\right)$
$\displaystyle-\left(b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\prod_{m=1,m\neq
h}^{\ell}\mathbb{E}_{\mathbb{Q}}\left[\delta_{m}\right]\mathbb{E}_{\mathbb{Q}}\left[z_{i\ell}^{2}\right]\right)\delta_{h}$
$\displaystyle=$
$\displaystyle\left(\frac{n(p-h+1)}{2}+a_{2}-1\right)\log\left(\delta_{h}\right)$
$\displaystyle-\left(b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\left[\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right]\prod_{m=2,m\neq
h}^{\ell}\left[\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\right]\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]\right)\delta_{h}$
since $\delta_{h}$ is bounded between $[1,\infty)$, $\displaystyle\delta_{h}$
$\displaystyle\sim\text{Gam}^{\top}\left(\frac{n(p-h+1)}{2}+a_{2},b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\left[\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right]\prod_{m=2,m\neq
h}^{\ell}\left[\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\right]\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]\right)$
$\displaystyle\therefore\tilde{a}_{2}^{(h)}$
$\displaystyle=\frac{n(p-h+1)}{2}+a_{2},\quad\tilde{b}_{2}^{(h)}=b_{2}+\frac{1}{2}\sum_{i=1}^{n}\sum_{\ell=h}^{p}\left[\frac{\tilde{a}_{1}}{\tilde{b}_{1}}\right]\prod_{m=2,m\neq
h}^{\ell}\left[\frac{\tilde{a}_{2}^{m}}{\tilde{b}_{2}^{m}}\right]\left[\tilde{z}_{i\ell}^{2}+\tilde{\omega}_{\ell}^{-1}\right]$

##  Appendix D Latent positions of the worm network

Figure 7 shows the initial latent positions for the first 3 dimensions under
the multidimensional scaling for $p=5$ for the worm nervous system network
while Figure 8 shows the posterior mean latent positions under the MCMC-LSPM
with $p=3$ for the worm nervous system network.

Figure 7: Initial latent positions estimates for the first 4 dimensions under
the multidimensional scaling for $p=5$ for the worm nervous system network.
Figure 8: Posterior mean latent positions under the MCMC-LSPM with $p=4$ for
the worm nervous system network.

## References

  * Aliverti and Durante, (2019) Aliverti, E. and Durante, D. (2019).  Spatial modeling of brain connectivity data via latent distance models with nodes clustering.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 12:185–196.
  * Aliverti and Russo, (2021) Aliverti, E. and Russo, M. (2021).  Stratified stochastic variational inference for high-dimensional network factor model.  Journal of Computational and Graphical Statistics, 31:502–511.
  * Bishop, (2006) Bishop, C. M. (2006).  Pattern recognition and machine learning.  Springer.
  * Blei et al., (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017).  Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112:859–877.
  * Burden and Faires, (1985) Burden, R. L. and Faires, J. D. (1985).  Numerical analysis. Prindle, Weber & Schmidt, 3rd ed edition.
  * Chu et al., (2021) Chu, A. M. Y., Chan, T. W. C., So, M. K. P., and Wong, W.-K. (2021).  Dynamic network analysis of covid-19 with a latent pandemic space model.  International Journal of Environmental Research and Public Health, 18(6):3195–3195.
  * Consonni and Marin, (2007) Consonni, G. and Marin, J.-M. (2007).  Mean-field variational approximate bayesian inference for latent variable models.  Computational Statistics & Data Analysis, 52:790–798.
  * de Reus and van den Heuvel, (2013) de Reus, M. A. and van den Heuvel, M. P. (2013).  Rich club organization and intermodule communication in the cat connectome.  Journal of Neuroscience, 33:12929–12939.
  * D’Angelo et al., (2023) D’Angelo, S., Alfò, M., and Fop, M. (2023).  Model-based clustering for multidimensional social networks.  Journal of the Royal Statistical Society Series A: Statistics in Society, 186:481–507.
  * Fletcher, (1964) Fletcher, R. (1964).  Function minimization by conjugate gradients.  The Computer Journal, 7:149–154.
  * Gollini and Murphy, (2016) Gollini, I. and Murphy, T. B. (2016).  Joint modeling of multiple network views.  Journal of Computational and Graphical Statistics, 25:246–265.
  * (12) Gwee, X. Y., Gormley, I. C., and Fop, M. (2023a).  A latent shrinkage position model for binary and count network data.  Bayesian Analysis, pages 1 – 29.
  * (13) Gwee, X. Y., Gormley, I. C., and Fop, M. (2023b).  Model-based clustering for network data via a latent shrinkage position cluster model.  arXiv.
  * Handcock et al., (2007) Handcock, M. S., Raftery, A. E., and Tantrum, J. M. (2007).  Model-based clustering for social networks.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 170:301–354.
  * Hoff et al., (2002) Hoff, P. D., Raftery, A. E., and Handcock, M. S. (2002).  Latent space approaches to social network analysis.  Journal of the American Statistical Association, 97:1090–1098.
  * Jaakkola and Jordan, (2000) Jaakkola, T. S. and Jordan, M. I. (2000).  Bayesian parameter estimation via variational methods.  Statistics and Computing, 10:25–37.
  * Jarrell et al., (2012) Jarrell, T. A., Wang, Y., Bloniarz, A. E., Brittin, C. A., Xu, M., Thomson, J. N., Albertson, D. G., Hall, D. H., and Emmons, S. W. (2012).  The connectome of a decision-making neural network.  Science, 337:437–444.
  * Jordan et al., (1999) Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999).  An introduction to variational methods for graphical models.  Machine Learning, 37:183–233.
  * Knowles and Ghahramani, (2011) Knowles, D. and Ghahramani, Z. (2011).  Nonparametric bayesian sparse factor models with application to gene expression modeling.  The Annals of Applied Statistics, 5.
  * Krivitsky and Handcock, (2008) Krivitsky, P. N. and Handcock, M. S. (2008).  Fitting position latent cluster models for social networks withlatentnet.  Journal of Statistical Software, 24.
  * Krivitsky and Handcock, (2022) Krivitsky, P. N. and Handcock, M. S. (2022).  latentnet: Latent Position and Cluster Models for Statistical Networks.  The Statnet Project (https://statnet.org).  R package version 2.10.6.
  * Liu and Chen, (2023) Liu, Y. and Chen, Y. (2023).  Variational inference for latent space models for dynamic networks.  Statistica Sinica, 32.
  * Lyu et al., (2023) Lyu, Z., Dong, X., and Zhang, Y. (2023).  Latent space model for higher-order networks and generalized tensor decomposition.  Journal of Computational and Graphical Statistics, pages 1–17.
  * Mardia et al., (2008) Mardia, K. V., Kent, J. T., and Bibby, J. M. (2008).  Multivariate analysis.  Academic Press.
  * Petersen and Pedersen, (2012) Petersen, K. B. and Pedersen, M. S. (2012).  The matrix cookbook.  Technical University of Denmark.  Version 20121115.
  * Raftery et al., (2012) Raftery, A. E., Niu, X., Hoff, P. D., and Yeung, K. Y. (2012).  Fast inference for the latent space network model using a case-control approximate likelihood.  Journal of Computational and Graphical Statistics, 21:901–919.
  * Rastelli et al., (2018) Rastelli, R., Maire, F., and Friel, N. (2018).  Computationally efficient inference for latent position network models.  arXiv.
  * Ročková and George, (2016) Ročková, V. and George, E. I. (2016).  Fast bayesian factor analysis via automatic rotations to sparsity.  Journal of the American Statistical Association, 111:1608–1622.
  * Salter-Townshend and Murphy, (2013) Salter-Townshend, M. and Murphy, T. B. (2013).  Variational bayesian inference for the latent position cluster model for network data.  Computational Statistics & Data Analysis, 57:661–671.
  * Scannell et al., (1995) Scannell, J., Blakemore, C., and Young, M. (1995).  Analysis of connectivity in the cat cerebral cortex.  The Journal of Neuroscience, 15:1463–1483.
  * Sewell and Chen, (2015) Sewell, D. K. and Chen, Y. (2015).  Latent space models for dynamic networks.  Journal of the American Statistical Association, 110:1646–1657.
  * Sewell and Chen, (2017) Sewell, D. K. and Chen, Y. (2017).  Latent space approaches to community detection in dynamic networks.  Bayesian Analysis, 12:351–377.
  * Spencer et al., (2022) Spencer, N. A., Junker, B. W., and Sweet, T. M. (2022).  Faster mcmc for gaussian latent position network models.  Network Science, 10:20–45.
  * Yang et al., (2020) Yang, Y., Pati, D., and Bhattacharya, A. (2020).  $\alpha$-variational inference with statistical guarantees.  The Annals of Statistics, 48.