11institutetext: Computer Science dep., ETH Zurich, Ramistrasse 101, 8092
Zurich, Switzerland 22institutetext: Swiss Data Science Center, Turnerstrasse
1, 8092 Zurich, Switzerland
22email<EMAIL_ADDRESS>
Shirin Goshtasbpour

# Optimization of Annealed Importance Sampling Hyperparameters

Shirin Goshtasbpour
1122 0000-0002-4902-165X Fernando Perez-Cruz 11 2 2 0000-0001-8996-5076

###### Abstract

Annealed Importance Sampling (AIS) is a popular algorithm used to estimates
the intractable marginal likelihood of deep generative models. Although AIS is
guaranteed to provide unbiased estimate for any set of hyperparameters, the
common implementations rely on simple heuristics such as the geometric average
bridging distributions between initial and the target distribution which
affect the estimation performance when the computation budget is limited.
Optimization of fully parametric AIS remains challenging due to the use of
Metropolis-Hasting (MH) correction steps in Markov transitions. We present a
parameteric AIS process with flexible intermediary distributions and optimize
the bridging distributions to use fewer number of steps for sampling. A
reparameterization method that allows us to optimize the distribution sequence
and the parameters of Markov transitions is used which is applicable to a
large class of Markov Kernels with MH correction. We assess the performance of
our optimized AIS for marginal likelihood estimation of deep generative models
and compare it to other estimators.

###### Keywords:

Annealed Importance Sampling Partition function estimation Generative models

##  1 Introduction

Deep generative models, developed over a decade, are now capable of simulating
complex distributions in high dimensional space and synthesizing high quality
samples in various domains such as natural images, text and medical data. Many
of these models are built with the assumption that data $x\in\mathcal{X}$
resides close to a low dimension manifold. In this case, data can be
represented using latent variable $z$ drawn from prior distribution on
$\mathbb{R}^{d}$ with density $p(z)$ and the observation model
$p(x|z)=p(x|f_{\theta}(z))$ is given by a parametric mapping
$f_{\theta}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\to\Gamma$ where
$\Gamma$ is the space of parameters of the likelihood $p(x|z)$.

For complex $f_{\theta}$ mappings, evaluation of log marginal likelihood

$\displaystyle\log p(x)=\log\int p(z)p(x|z)dz$

is intractable and the posterior density $p(z|x)\propto p(x|z)p(z)$ is known
only up to a constant normalization factor111We use normalization factor and
marginal likelihood $p(x)$ interchangeably in this paper. In this case, latent
variable inference, probabilistic model evaluation and model comparison are
performed using variational approximation of the posterior [3, 24, 38] or via
sampling methods like Markov Chain Monte Carlo (MCMC) [4], Nested sampling
[43, 32], Sequential Monte Carlo (SMC) [8] and Annealed Importance Sampling
(AIS) [34]. In Variational Inference (VI) the posterior is approximated with
the most similar probability density from the family of distributions
$\mathcal{Q}=\\{q_{\phi}(z|x)\mathrel{\mathop{\ordinarycolon}}\phi\in\mathbb{\phi}\\}$
indexed by parameter $\phi$. Success of VI depends on sufficient expressivity
of members of $\mathcal{Q}$ and our ability to find the optimal member of this
family by, for instance, maximizing the Evidence Lower BOund (ELBO)

$\displaystyle\mathcal{L}(\phi)=\int\log\left(\frac{p(x,z)}{q_{\phi}(z|x)}\right)q_{\phi}(z|x)dz\leq\log
p(x)$ (1)

where the equality happens only when $q_{\phi}(z|x)=p(z|x)\in\mathcal{Q}$.

In MCMC we use Markov kernels with unique stationary distribution $p(z|x)$ and
sample a Markov chain $\left(z_{k}\right)_{k\in[M]}$ by iterative application
of this kernel on an initial particle $z_{0}\sim q_{\phi}(z|x)$. Contrary to
parametric VI, under mild assumptions, it is theoretically guaranteed that
distribution of $z_{M}$ converges to the target distribution $p(z|x)$ as $M$
goes to infinity. However, if the posterior has multiple modes or heavy tails,
convergence can require large number of iterations and therefore, be
computationally prohibitive. Among the sampling methods, SMC and AIS are of
particular interest as they produce unbiased estimation of marginal likelihood
regardless of the computational budget by assigning importance weights to
samples [40, 50, 31]. In these algorithms an auxiliary sequence of
distributions, $\left(\pi_{k}\right)_{k\in[M]}$, is used to bridge a simple
proposal distribution with density $\pi_{0}(z)=q_{\phi}(z|x)$ and the target
density $\pi_{M}(z)=p(z|x)$. This sequence is defined via unnormalized
densities $\left(\gamma_{k}\right)_{k\in[M]}$ where
$\pi_{k}(z)=\gamma_{k}(z)/Z_{k}$ and $\gamma_{0}(z)=\pi_{0}(z)$ and
$\gamma_{M}(z)=p(x,z)$. The algorithm produces $N$ Markov chain samples
$\left(z_{k}^{j}\right)_{k\in[M]}$ for $j\in[N]$ and their corresponding
importance weights denoted by $w^{j}$ as follows: Initially $z_{0}^{j}$ is
sampled from $\gamma_{0}$. Then $z_{k}^{j}$ is approximately sampled from
$\gamma_{k}$ using a Markov kernel (typically a gaussian distribution around
$z_{k-1}^{j}$ or a transition with invariant distribution $\gamma_{k}$)222SMC
has an additional resampling step to draw exact samples from $\gamma_{k}$ to
reduce the variance of importance weight although it sometimes results in
insufficient sample diversity. The marginal likelihood is approximated with
Monte Carlo method

$\displaystyle
p(x)=\mathbb{E}[w]\approx\hat{p}^{N}(x)\coloneqq\frac{1}{N}\sum_{j\in[N]}w^{j}$
(2)

where the expectation is taken over joint distribution of the Markov chains.
Variance of $\hat{p}^{N}(x)$ depends on the selected density sequence and
hyperparameters especially when computation resources are scarce. However,
little work is available on how to optimize the intermediary distributions. In
this paper we focus on optimization of AIS hyperparameters for more accurate
estimation of log marginal likelihood with only a few annealing distributions.
We optimize the sequence of distributions between proposal and target
distributions as well as the hyperparameters of Markov kernels. We have made
the following contributions:

  * •

We propose the parameterization of a continuous time density path between the
proposal and target distributions which can define general density functions.
Contrary to [24, 25] the densities used do not need to be normalized and we
don’t require their exact samples since we use sampling algorithms to
gradually transition between the intermediary distributions.

  * •

To optimize the bridging distributions we minimize Jefferys and inverse KL
divergences between the distributions of AIS process and its auxiliary reverse
process (defined in Section 2.1). Inspired from recent developments in Monte
Carlo VI (MCVI) methods, we reparameterize particle states with a mixture of
Normalizing Flows and derive a low variance consistent estimation of the
gradient of our optimization objective.

  * •

We show that Jefferys divergence captures both bias and variance of the
estimation while optimization of inverse KL divergence reduces the bias of
estimation. We further evaluate our method on deep generative models with
different training procedures.

The combination of both procedures allows reducing the variance of the
importance sampling weights providing more accurate estimates with fewer
particles and intermediary distributions. The rest of this paper is organized
as follows: Section 2 is dedicated to reviewing vanilla AIS algorithm. In
Section 3 we present our parameterization of AIS process which results in
flexible bridging distributions and motivate out optimization objective. In
Section 4 we illustrate a reparameterization method and derivation of the
objective gradient estimates. Finally, we analyze the accuracy of marginal
likelihood estimation and its variance on synthetic and image datain Section
6.

##  2 Background

In this Section we give a brief introduction to AIS algorithm and its popular
adaptive versions. For the rest of this section, we assume that the
observation $x$ is fixed and $\pi(z)=p(z|x)$ is the target density function
which we can evaluate up to a normalizing constant $\tilde{\pi}(z)=p(z,x)$
where $\pi(z)=\tilde{\pi}(z)/Z$ and $Z=p(x)$. We also define a proposal
distribution with normalized density $\gamma_{0}(z)$ which is easy to sample
and evaluate, such as the variational posterior $\gamma_{0}(z)=q_{\phi}(z|x)$.

###  2.1 Recap of Annealed Importance Sampling

To sample from the target distribution, AIS algorithm requires a sequence of
distributions which change from the proposal distribution to the target. To
match with our parameterization of the bridging distributions, we consider the
generalization of this sequence to a continuously indexed path of density
functions. With some abuse of notation, let
$\gamma\mathrel{\mathop{\ordinarycolon}}[0,1]\times\mathbb{R}^{d}\to\mathbb{R}$
with mapping $(t,z)\mapsto\gamma(t,z)$ denote a path of density functions
indexed by $t\in[0,1]$, with fixed endpoints
$\gamma(0,\cdot)=q_{\phi}(\cdot|x)$ and $\gamma(1,\cdot)=\tilde{\pi}(\cdot)$.
Assume the density functions have non-increasing support (i.e.
$\text{supp}(\gamma(t^{\prime},\cdot))\subseteq\text{supp}(\gamma(t,\cdot))$
for $0\leq t\leq t^{\prime}\leq 1$). Let $t_{k}=k/M$ for $k\in\\{0,...,M\\}$
be the linear schedule in $[0,1]$ with which the path of density functions is
discretized and let $\gamma_{k}(\cdot)=\gamma(t_{k},\cdot)$. We assume that
the $\gamma(\cdot,z)$ curves are flexible enough to absorb any increasing
mapping $\beta\mathrel{\mathop{\ordinarycolon}}[0,1]\to[0,1]$ (representing a
discretization of $[0,1]$ at $\left(\beta(t_{k})\right)_{k\in\\{0,...,M\\}}$).

AIS algorithm starts by sampling the initial particles
$z_{0}\sim\gamma_{0}(z)$ and evaluating the initial weights $w_{0}=1$. Then,
particles are gradually moved in iteration $k$ using Markov kernel
$\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})$ such that the distribution
of particles $z_{M}$ approaches the target distribution $\pi(z)$. The particle
trajectory distribution in AIS has the joint density
$\overrightarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M})=\gamma_{0}(z_{0})\prod_{k}\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})$.
In each iteration a second Markov kernel
$\overleftarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})$ is used to approximate the
probability density of backward transition $z_{k}\to z_{k-1}$ such that the
ratio
$w_{k}(z_{k-1},z_{k})=\gamma_{k}(z_{k})\overleftarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})/\gamma_{k-1}(z_{k-1})\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})$
is well-defined on $\mathbb{R}^{d}$ and the importance weights are adjusted by
multiplying $w_{k}(z_{k-1},z_{k})$ with the current weights. The final
importance weights can be evaluated as

$\displaystyle
w(z_{0\mathrel{\mathop{\ordinarycolon}}M})=\frac{\tilde{\pi}(z_{M})\overleftarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M-1}|z_{M})}{\overrightarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M})}=\prod_{k=1}^{M}w_{k}(z_{k-1},z_{k})$

where
$\overleftarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M-1}|z_{M})=\prod_{k}\overleftarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})$.

The optimal backward process
$\overleftarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M-1}|z_{M})$
with zero variance estimator $\hat{p}^{N}(x)$ is given by inverse kernels

$\displaystyle\overleftarrow{\mathcal{T}}^{*}_{k}(z_{k-1}|z_{k})=\frac{\overrightarrow{q}_{\phi,k-1}(z_{k-1})\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})}{\overrightarrow{q}_{\phi,k}(z_{k})}$

where $\overrightarrow{q}_{\phi,k}(z_{k})$ is the marginalization of
$\overrightarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M})$ over
all $z_{m}$ except $m\neq k$. However it is intractable to use
$\overleftarrow{\mathcal{T}}^{*}_{k}(z_{k-1}|z_{k})$ for sampling [8].

Alternatively, we can use a heuristically fixed $\left(\gamma_{k}\right)_{k}$
to guide the samples towards the target’s basins of energy. One may choose
$\overrightarrow{\mathcal{T}}_{k}$ from popular MCMC kernels like Random Walk
Metropolis-Hastings (RWMH), Hamiltonian Monte Carlo (HMC), Metropolis Adjusted
Langevin Algorithm (MALA) or Unadjusted Langevin Algorithm (ULA) with
$\pi_{k}$ as their stationary distribution [45, 46]. MH-corrected kernels
typically do not admit a transition probability density. However, due to
detailed balance
$\gamma_{k}(z_{k})\overrightarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})=\gamma_{k}(z_{k-1})\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})$
the kernel $\overrightarrow{\mathcal{T}}_{k}$ is reversible with respect to
$\gamma_{k}$ and with backward kernel set to its reversal
$\gamma_{k}(z_{k-1})\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})/\gamma_{k}(z_{k})$
we can obtain well-defined importance weight updates
$w_{k}(z_{k-1},z_{k})=\gamma_{k}(z_{k-1})/\gamma_{k-1}(z_{k-1})$ [34]. When
$\overrightarrow{\mathcal{T}}_{k}$ is not $\gamma_{k}$-invariant, an
approximate reversal kernel with same assumption may be constructed from
$\overleftarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})=\overrightarrow{\mathcal{T}}_{k}(z_{k-1}|z_{k})$
while the weight updates preserves their original form [46].

In another approach, reparameterizable forward and backward transition kernels
can be used where the parameters are optimized to push samples from one
predetermined bridging distribution to the next, which limits the complexity
of the transition probability density considerably [31, 49, 1].

##  3 Optimized Path Annealed Importance Sampling

Typical implementations of AIS use predefined heuristics for the sequence
$\left(\gamma_{k}\right)_{k}$ and the schedule $\beta$. Popular choices
include the geometric average density path
$\gamma(t,z)=\tilde{\pi}(z)^{\beta(t)}q_{\phi}(z|x)^{\beta(1-t)}$ with linear
$\beta(t)=t$ or geometric $\beta(t)=\alpha^{\log t}$ schedule, for some real
$\alpha>1$. In contrast, we propose to optimize the bridging distribution path
$\gamma(t,\cdot)$ in a class of positive parametric functions. If the class is
sufficiently large the schedule $\beta$ will be implicit in $\gamma$ and
therefore we can fix the schedule to be linear and only focus on
parameterizing the $\gamma(t,\cdot)$. We use Deep Neural Networks (DNN) to
define our parametric density function path. We name our procedure Optimized
Path Annealed Importance Sampling (OP-AIS).

###  3.1 Parameterized AIS

We opt to parameterize the time dependent path $\gamma$ as opposed to
parameterizing the intermediary distributions of each iteration individually.
The main reason for this choice is that it results in a density path that is
continuous in $t$ and allows parameters to be shared by intermediate
distributions. We use a continuous mapping
$u_{\phi}\mathrel{\mathop{\ordinarycolon}}\mathcal{X}\times\mathcal{Z}\times[0,1]\to\mathbb{R}$
with additional terms that adjust the boundary constraints at $t\in\\{0,1\\}$.
We consider the following boundary adjustment which coincides with geometric
average of the proposal and target distributions when $u_{\phi}\equiv 0$.

$\displaystyle\log\gamma(t,z)=$ $\displaystyle u_{\phi}(x,z,t)$
$\displaystyle+(1-t)\left[-u_{\phi}(x,z,0)+\log\gamma_{0}(z)\right]$
$\displaystyle+t\left[-u_{\phi}(x,z,1)+\log\tilde{\pi}(z)\right]$ (3)

If $u_{\phi}(x,z,t)\in(-\infty,\infty)$ is well-defined for all
$x\in\mathcal{X}$, $z\in\mathbb{R}^{d}$ and $t\in(0,1)$ and the proposal is
chosen such that
$\text{supp}(\gamma(1,\cdot))\subseteq\text{supp}(\gamma(0,\cdot))$ then
$\text{supp}(\gamma(1,\cdot))=\text{supp}(\gamma(t,\cdot))$ for all
$t\in[0,1]$ and the importance weights updates are always defined using the
common gaussian or invariant transition kernels. Notably, we can apply this
parameterization to any arbitrary path of density functions by twisting the
path with multiplication of an arbitrary positive function
$\exp(-u_{\phi}(x,z,t))$ with correct boundary behavior.

###  3.2 Optimization for the parametric distribution

Given an unnormalized target distribution, Effective Sample Size (ESS) and its
conditional extension are the de facto measure used to adaptively tune AIS [8,
21, 22, 34]. ESS of $N$ parallel runs of AIS is given by

$\displaystyle\text{ESS}(w^{1},...,w^{N})=\frac{\left(\sum
w^{j}\right)^{2}}{N\sum(w^{j})^{2}}$

and is expected to reflect the estimator variance odds with exact samples from
$\pi(z)$ and importance sampling (2). ESS is sensitive to weight degeneracy
problem that is prominent in AIS and other sequential importance sampling
methods (when only a few particles contribute to the Monte Carlo approximate).
However, high ESS does not determine if the particles have sufficient
dispersion according to the target distribution. See [11] for further
discussion on ESS and its defects.

Grosse et. al. in [16] derived a bound on $\log Z$ by running the AIS
algorithm in forward and backward directions. Original direction of AIS
evaluates $\log\overrightarrow{w}=\log
w(z_{0\mathrel{\mathop{\ordinarycolon}}M})$ which in expectation lower bounds
$\log Z$.

$\displaystyle\mathbb{E}_{\overrightarrow{q}_{\phi,M}}[\log\overrightarrow{w}]$
$\displaystyle\leq\log Z$ (4)

Drawing sample trajectories $z_{0\mathrel{\mathop{\ordinarycolon}}M}^{\prime}$
from the reverse direction of AIS algorithm defined by the ancestral sampling
in joint distribution
$\pi(z_{M}^{\prime})\overleftarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M-1}^{\prime}|z_{M}^{\prime})$
we can compute random variable $\log\overleftarrow{w}=\log
w(z_{0\mathrel{\mathop{\ordinarycolon}}M}^{\prime})$ which upper bounds $\log
Z$ in expectation. The gap between the two expectations in forward and
backward processes is called the BiDirectional Monte Carlo (BDMC) gap and is
equivalent to Jefferys divergence between these distributions.

$\displaystyle\mathcal{L}^{\text{BDMC}}(x)$
$\displaystyle=\mathbb{E}_{\pi\overleftarrow{q}_{\phi,M}}[\log\overleftarrow{w}]-\mathbb{E}_{\overrightarrow{q}_{\phi,M}}[\log\overrightarrow{w}]$
(5)
$\displaystyle=\text{D}_{\text{KL}}(\pi\overleftarrow{q}_{\phi,M}||\overrightarrow{q}_{\phi,M})+\text{D}_{\text{KL}}(\overrightarrow{q}_{\phi,M}||\pi\overleftarrow{q}_{\phi,M})$

This bound is also frequently used to assess the accuracy of the AIS estimator
[15, 50, 20]. Unbiased estimation of this bound requires exact samples from
the target distribution. Therefore, an approximation is made replacing the
data samples with the synthesized samples by the generative model, where we
have access to an underlying latent variable. This approximation results in
more accurate log-likelihood estimation of samples closer to the model which
affects the objectivity of the test. Assuming a long chain and close to
perfect transitions the final particles from AIS forward process can be used
to approximate this bound instead. However, with low computation budget this
marginal distribution may be far off from the target especially in the initial
optimization iterations. Although, the BDMC gap achieves a natural trade-off
between minimization of the bias in ELBO bound and $\log$ empirical variance
of the importance weights (Figure 2), the two terms in
$\mathcal{L}^{\text{BDMC}}(x)$ have contradicting gradients resulting in
unstable optimization of hyperparameters with gradient based methods.

In order to achieve low variance importance weights we need to match the
distributions of the forward and backward processes. It is common to minimize
the inverse KL-divergence between these distributions similar to maximization
of (1) as this results in reducing the bias in (4).

$\displaystyle\mathcal{L}^{\text{AIS}}(x)=\text{D}_{\text{KL}}(\overrightarrow{q}_{\phi,M}||\tilde{\pi}\overleftarrow{q}_{\phi,M})=-\mathbb{E}_{\overrightarrow{q}_{\phi,M}}[\log\overrightarrow{w}]$

We can evaluate an unbiased estimate of inverse KL divergence using samples
from $\overrightarrow{q}_{\phi,M}$.

As was shown by authors of [9, 28] under some assumptions, the variance of
importance weights can be controlled by the inverse KL-divergence between two
distributions. In particular, with some adjustment to the AIS setting we have

$\displaystyle\text{D}_{\text{KL}}(\overrightarrow{q}_{\phi,M}||\tilde{\pi}\overleftarrow{q}_{\phi,M})\approx\frac{1}{2Z^{2}}\text{Var}_{\overrightarrow{q}_{\phi,M}}[w(z_{0\mathrel{\mathop{\ordinarycolon}}M})]$

Therefore, by minimizing the inverse KL-divergence we minimize the log
marginal likelihood bias and the variance of the importance weights,
simultaneously. The equality holds when
$\pi(z_{k})\overleftarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M-1}|z_{M})=\overrightarrow{q}_{\phi,M}(z_{0\mathrel{\mathop{\ordinarycolon}}M})$
which results in zero importance weight variance and
$\mathcal{L}^{\text{AIS}}(x)=0$.

##  4 Stochastic Optimization for OP-AIS

In order to evaluate the gradients of $\mathcal{L}^{\text{AIS}}(x)$ with
respect to the parameters $\phi$ we use a combination of reparameterization
trick and score function estimation introduced by [42] and [31]. The MH-
corrected Markov kernel $\overrightarrow{\mathcal{T}}_{k}(z_{k}|z_{k-1})$ is
given by

$\displaystyle\overrightarrow{\mathcal{T}}_{k}(z_{k}$
$\displaystyle|z_{k-1})=r(z_{k}|z_{k-1})\alpha_{k}(z_{k-1},z_{k})$
$\displaystyle+\int\left(1-\alpha_{k}(z|z_{k-1})\right)r(z|z_{k-1})dz\delta_{z_{k-1}}(z_{k})$

with $r(.|z_{k-1})$ denoting the conditional probability density of proposed
particle in the $k$th iteration,
$\alpha_{k}(z_{k-1},z)=\left(1\land\frac{\gamma_{k}(z)r(z_{k-1}|z)}{\gamma_{k}(z_{k-1})r(z|z_{k-1})}\right)$
the acceptance rate at this iteration and $\delta_{z}$ denoting the Dirac
measure at $z$. We assume that the proposal $r(z|z_{k-1})$ is
reparameterizable with a auxiliary random variable
$\epsilon_{k}\sim\eta(\epsilon)$ and transformation
$z=T_{k}(\epsilon_{k},z_{k-1})$ which is surjective in second argument. After
transformation the proposal state is accepted with probability
$\alpha_{k}(z_{k-1},T_{k}(\epsilon_{k},z_{k-1}))$ (denoted with $\alpha_{k}$
as a shorthand) and we set $z_{k}=T_{k}(\epsilon_{k},z_{k-1})$. Otherwise, we
keep the old state and set $z_{k}=z_{k-1}$. We determine the proper action by
sampling a binary random variable
$a_{k}\sim\text{Bern}(\alpha_{k}(z_{k-1},T_{k}(\epsilon_{k},z_{k-1})))$.

For instance, for a Random Walk kernel $\overrightarrow{\mathcal{T}}_{k}$, we
assume that $\eta(\epsilon)$ is the probability density function of standard
normal distribution and $T_{k}(\epsilon,z)=z+\Sigma^{1/2}\epsilon$ for a
positive definite covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$. We
refer the readers to the derivation of transformations corresponding to HMC,
MALA and ULA by [45] for further details.

Using the above reparameterization we get

$\displaystyle\mathcal{L}^{\text{AIS}}(x)=-\mathbb{E}_{\gamma_{0}(z_{0})\prod_{k}\eta(\epsilon_{k})\prod_{k}\zeta_{k}(a_{k}|z_{0},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k},a_{1\mathrel{\mathop{\ordinarycolon}}k-1})}\left[\log
w\right]$

for
$\zeta_{k}(a^{k}|z_{0},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k},a_{1\mathrel{\mathop{\ordinarycolon}}k-1})=\alpha_{k}^{a_{k}}\left(1-\alpha_{k}\right)^{a_{k}}$.

The gradient of $\mathcal{L}^{\text{AIS}}(x)$ is given by

$\displaystyle\nabla_{\phi}\mathcal{L}^{\text{AIS}}(x)=$
$\displaystyle-\mathbb{E}\left[\nabla_{\phi}\log w\right]$
$\displaystyle-\sum_{k}\mathbb{E}\left[\log
w\nabla_{\phi}\log\zeta_{k}(a_{k}|z_{0},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k},a_{1\mathrel{\mathop{\ordinarycolon}}k-1})\right]$

We can derive $\mathcal{L}^{\text{AIS}}(x)$ approximation using the particle
and importance weights generated from the AIS .

$\displaystyle\nabla_{\phi}\hat{\mathcal{L}}^{\text{AIS}}(x)=$
$\displaystyle-\frac{1}{N}\sum_{i=1}^{N}\nabla_{\phi}\log w^{i}$
$\displaystyle+\underbrace{-\frac{1}{N}\sum_{i=1}^{N}\sum_{k}\log
w^{i}\nabla_{\phi}\log\zeta(a_{k}^{i}|z_{0}^{i},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k}^{i},a_{1\mathrel{\mathop{\ordinarycolon}}k-1}^{i})}_{\nabla_{\phi}\hat{\mathcal{L}}^{\text{score}}(x)}$

$\nabla_{\phi}\hat{\mathcal{L}}^{\text{AIS}}(x)$ is a strongly consistent
estimator of $\nabla_{\phi}\mathcal{L}^{\text{AIS}}(x)$ [27].

Score function estimator is notorious for its high variance. It is a standard
practice to omit $\nabla_{\phi}\mathcal{L}^{\text{score}}$ in parameter
updates for more stable optimization. However, this results in biased gradient
estimate which adversely effects our estimation in latent variable model
experiements. To reduce the high variance of
$\nabla_{\phi}\hat{\mathcal{L}}^{\text{score}}(x)$ we can employ the variation
reduction techniques [33, 36]. For instance, in our experiments we use the
common particles and weights to estimate the two expectations in
$\mathcal{L}^{\text{AIS}}(x)$ and the one-out control variates proposed by
[30] to substitute $\nabla_{\phi}\hat{\mathcal{L}}^{\text{score}}(x)$ with

$\displaystyle\sum_{i=1}^{N}\left[\frac{NR^{i}-\sum_{j}R^{j}}{N-1}\right]\sum_{k}\nabla_{\phi}\log\zeta_{k}(a_{k}^{i}|z_{0}^{i},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k}^{i},a_{1\mathrel{\mathop{\ordinarycolon}}k-1}^{i})$

where $R^{i}=-\log w^{i}/N$.

1

input : Target $\tilde{\pi}$, proposal density $\gamma_{0}$, and $\eta$

output : Parameters $\phi$

2 for  _$K$ epochs_ do

3 Set $\Delta\leftarrow 0$

4 for  _$i=1$ to $N$_ do

5 Set $\text{logw}^{i}\leftarrow 0$, $\text{log}\zeta^{i}\leftarrow 0$

6 Draw $z^{i}\sim q_{0}(z)$

7 Draw
$\epsilon_{1\mathrel{\mathop{\ordinarycolon}}M}^{i}\sim\prod_{k}\eta(\epsilon_{k})$

8

9 for  _$k=1$ to $M$_ do

10 for  _$i=1$ to $N$_ do

11 Set
$\text{logw}^{i}\leftarrow\text{logw}^{i}+\log\frac{\gamma_{k}(z^{i})}{\gamma_{k-1}(z^{i})}$

12 Draw $a_{k}^{(n)}\sim\text{Bern}(\alpha_{k})$

13 Set
$\text{log}\zeta^{i}\leftarrow\text{log}\zeta^{i}+\log\zeta(a_{k}^{i}|z_{0}^{i},\epsilon_{1\mathrel{\mathop{\ordinarycolon}}k}^{i},a_{1\mathrel{\mathop{\ordinarycolon}}k-1}^{i})$

14 Set $z^{i}\leftarrow
T_{k}(\epsilon_{k}^{i},z^{i})^{a_{k}^{i}}(z^{i})^{1-a_{k}^{i}}$

15

16

17 Get $R^{i}$ for $i\in[N]$

18 $\phi\leftarrow\phi-\kappa\nabla_{\phi}\hat{\mathcal{L}}^{\text{AIS}}(x)$

Algorithm 1 Optimization of parameterized AIS

Pseudocode of the optimization algorithm is given in Algorithm (1).

##  5 Literature Review

We can interpret AIS as the discretization of the path sampling algorithm
where the thermodynamics integral identity holds

$\displaystyle\log Z$
$\displaystyle=\int_{0}^{1}\mathbb{E}_{\gamma(t,z)/Z_{t}}\left[\frac{d}{dt}\log\gamma(t,z)\right]dt$

for $Z_{t}=\int\gamma(t,z)dz$. Gelman et. al. in [13] derived the optimal path
and scaling of the bridging distributions for minimum variance estimation.
Evaluation of optimal path is intractable in general. In turn, authors of [6,
2] provided closed form optimal path for special cases. Viewing the importance
weight updates as the finite difference approximations, the authors in [17]
derive the asymptotic bias of $\log w$ estimator and propose moment average
path between members of exponential family. A variational representation of
the annealing distributions minimizing weighted $\alpha$-divergences between
proposal and target distributions are derived in [17, 29]. In [26], the
authors minimize the asymptotic variance of estimator $\log w$ for a given
parametric path of distributions.

The analysis in these works are derived based on the assumption of reaching
equilibrium in each iteration and independence of particles. This is not the
case when transition kernels are not invertable. On the other hand, exhaustive
heuristic optimization to tune AIS hyperparameters is laborious [40]. Other
methods to adapt AIS hyperparameters include monitoring the acceptance rate in
Markov transitions to adjust its parameters [21, 46] and designing problem
specific or heuristic intermediary distribution sequences [21, 40, 29, 17]. On
the contrary, while our method requires training before evaluation, we can
amortize the overhead by parameter sharing and reduce the computation cost by
using only a small subset of validation samples. As a consequence, we find an
efficient and shorter annealing process for sampling. Our aim is to enable AIS
to achieve high accuracy despite limited computation resources.

Recently, Variational Inference and Monte Carlo methods have been combined to
increase the flexibility of parametric variational family $\mathcal{Q}$ and
benefit from the convergence properties of sampling methods. Authors of [42,
48, 7, 19, 45] combine MCMC and VI by extending the space of latent variables
$\left(z_{k}\right)_{k\in[M]}$. A framework with parametric transitions for
SMC algorithms and a surrogate lower bound ELBO are provided in [18, 47, 31].
For example, in [28], the authors derive an approximation of the asymptotic
variance of importance weights in Variational Filtering Objectives, and in
[49, 1], flows are used to overcome energy density barriers. A differentiable
alternative to MH-correcting step in HMC transitions are proposed [52], while
in [12] the authors derive the importance weights resulting from unadjusted
HMC process and its reversal.

Our work in particular bares more similarity to MCVAE proposed by [46] among
others. In MCVAE the authors use AIS with heuristic density path for low
variance estimation of the log marginal likelihood during training of a latent
variable model. Our goal here is to evaluate an already trained model by
optimizing the intermediary distributions and we argue that the learnt AIS
parameters by MCVAE are not efficient since the model is constantly changing
during training and the posterior approximator is known to overfit the model
[50].

##  6 Experiments

Here, we evaluate the proposed algorithm on complex synthetic 2d distributions
which are often used to benchmark generative models [37]. Target distributions
are illustrated in Figure 1.

Figure 1: Target distributions and particles of samplers with optimized path
(trained with KL and Jefferys divergence) with RWMH, HMC and MALA transitions
and $M=8$. Particles with larger weights are coded to show comparable scale to
maximum weight in each plot. Larger dots with warmer colors show particles
with larger weights, with red for maximum weight in each set and dark blue for
particles with 0.1 of the maximum weight

We compare our algorithm to vanilla AIS [34] with geometric average path and
geometric and linear schedule. We use RWMH, HMC and MALA transitions and
normal proposal with learnable mean and diagonal covariance matrix. $M$
between 2-128 bridging distributions and $256$ samples are used for training.
$u_{\phi}$ is implemented with a DNN with one hidden layer with 4 dimensions
and LeakyReLU(0.01) nonlinearty except at the output layer.

The parametric AIS is trained with Adam optimizer [23] with learning rate of
0.03 and betas = (0.5, 0.999) for 100 epochs. The RWMH kernels have local
normal proposal steps with learnable diagonal covariance matrix and the step
size in all transition kernels are trained along with other parameters. We use
cross validation to choose the step size in vanilla AIS from
$\\{0.01,0.05,0.1,0.5,1.\\}$. The normalization constants are estimated using
$N=4096$ samples. To act as the reference normalization constant we use
vanilla AIS with $M=1024$, 10 MCMC steps per iteration with 3 different step
sizes. All the code is written with Pytorch [35] and experiments are run on a
GeForce GTX 1080 Ti. Code is available here:
https://github.com/shgoshtasb/op_ais.

In Figure 1 we illustrate particles generated from the samplers trained on
each distribution with inverse KL divergence (PKL) and Jefferys divergence
(PJ) color coded by their weight for $M=8$. Warmer colors show particles with
higher weights which are main contributors to the partition function
estimation, while dark blue particles are less effective in the estimation. We
observe that training with KL divergence results in smaller bias, while using
Jefferys divergence objective results in more effective particles and lower
variance estimation in Target 3-6.

Figure 2: From left to right: $\log Z$ estimation with Eq. 4, log empirical
variance of importance weights and the Jefferys divergence of forward and
backward processes vs $M$. (Vg, Vl for vanilla AIS with geometric and linear
schedules, PJ and PKL for OP-AIS with Jefferys and KL-divergence optimization)

Figure 2 shows the log partition function estimation and the empirical
variance of importance weights for the mentioned samplers. We index the target
distributions in the same order they appear in the Figure 1. In most of the
experiments, for small number of bridging distributions ($M\in\\{2,4,8\\}$)
our method is able to improve over vanilla AIS algorithm with linear schedule
and achieves slightly tighter lower bound in comparison to geometric schedule
version. In comparison to training with inverse KL divergence, Jefferys
divergence improves the variance of importance weights. we speculate this to
be related to higher correlation of Jefferys divergence with weight variances.

Interestingly, in multimodal and heavy tail distributions, empirical variance
for small M remains below the variance of exhaustive AIS (Vl M=1024 in dashed
Turquoise) and grows with $M$ With only one MCMC step used in the experiments,
particles ultimately cover a small area and don’t reach the high density
regions in the target distribution, resulting in low empirical weight variance
despite the fact that the actual variance is much larger than our observation.
This is specially problematic with the geometric schedule as the few bridging
distributions are placed very close to the proposal and the particles are not
encouraged to move far. In the low $M$ setting as it is unlikely to observe
high weight particles, it is misleading to use ESS or even empirical
$\chi^{2}$-divergence to tune the hyperparameters, whereas the inverse KL
divergence provides better guidance for tuning. Increasing $M$ results in more
dispersed particles and higher variance. For large enough $M$ as the particles
can reach high density regions in the target, we observe that this variance
decreases with $M$. Variance in parametric samplers is more steady and has a
small growth in comparison due to the optimized choice of bridging
distributions for each M.

###  6.1 MNIST Dataset

Figure 3: VAE sample reconstruction with posterior samples

In this section, we assess the effectiveness of our method on a generative
models based on similar DNN decoder architecture trained on MNIST dataset. We
use the binarized MNIST dataset with Bernoulli observation model and data
dequantization, rescaled to [0,1] [40]. The decoder model we consider is the
"Wide" architecture used in [50] with 50-1024-1024-1024-784 Fully-Connected
(FC) layers, Tanh activations for hidden layers and Sigmoid for the output.
This model is trained as a Variational AutoEncoder (VAE) [24], Importance
Weighted AutoEncoder (IWAE) [5], Generative Adversarial Network (GAN) [14] and
Bidirectional GAN (BiGAN) [10]. We compare OP-AIS with $M\in\\{8,16\\}$ and
$N=16$ particles per data sample, HMC transition kernels with 1 leapfrog step
to MCMC and Vanilla AIS with $M=1024$ and HMC transitions with 16 leapfrog
steps. OP-AIS is trained with Adam optimizer with $0.01$ learning rate and
betas=(0.5, 0.999) for 100 epochs on 1400 samples and validated on 600 samples
for highest log likelihood estimation. Reported log-likelihood values are over
remaining 8000 test samples. Here, we use a symmetric encoder with Tanh
activation and Dropout(0.2) at each hidden layer to transform $x$, concatenat
with $z$, and input to $u_{\phi}$ with 32-hidden units and LeakyReLU(0.01)
activation. The encoder network is similar to the one used in [50] for the
"Wide" decoder setup. The training and evaluation time of OP-AIS (2h38 for
$M=16$ on average) is approximately similar to the sampling time required by
the MCMC and Vanilla AIS algorithms (2h54 on average).

Table 1: Negative Log likelihood estimation of generative models on binarized
MNIST

Model | Encoder | OP-AIS | AIS
---|---|---|---
$M\times$steps |  | $8\times 1$ | $16\times 1$ | $1024\times 16$ | $1024\times 10^{\text{Tuned}}$
VAE | $95.7$ | $216.7\pm 29.8$ | 193$\pm$31.9 | $215.6\pm 35.4$ |  87.8$\pm$5.8
IWAE | $88.5$ | 213.9$\pm$27.3 |  221.2$\pm$34.1 | $235.4\pm 35.8$ |  85.8$\pm$5.6
GAN |  | $805.9\pm 94.7$ | 785.4$\pm$67.0 | $902.9\pm 115.8$ |  415.3 $\pm$48.8
BiGAN |  | $900.5\pm 124.6$ | $885.3\pm 138.3$ | 864.2$\pm$98.1 |  373.2$\pm$44.6

Table 1 shows the negative log marginal likelihood estimation using posterior
samples from each sampling algorithm. The ground truth values are approximated
using a more exhaustive version of AIS with 10 tuning iterations for each
transition kernel on the 1400 separated data samples (AIS $1024\times
10^{\text{Tuned}}$) (8h12 average tuning + sampling time). Although the values
of aquired by ground truth and encoder network (when available) are far from
sampler estimations, estimates are more accurate than the sampling counterpart
with similar computation time. We also observe that sample reconstruction with
OP-AIS achieves similar samples to the test at least in style and are
comparable to the other baselines in Figure 3.

##  7 Conclusion

AIS yields an unbiased estimation for any annealing path, schedule and
transition operator. However, the variance of the Monte Carlo estimator
changes considerably with different choices. We propose an optimization
procedure to find the best set of intermediary distributions that matches the
distribution of the AIS particle trajectory to joint distribution of backward
process with target distribution as its marginal. We formulate AIS algorithm
with MH-corrections as a mixture of NFs and minimize the the divergence
between the two distributions.

It is important to mention that log-likelihood has to be used with caution for
model comparison, since a model with very good or bad generative sample
quality or memorization may have high log likelihood [44]. It is recommended
to use other sample quality measures (e.g. [51, 41]) and model comparison
methods (e.g. [20, 39]) along with AIS.

####  7.0.1 Acknowledgements

This work is supported by funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 813999 for this project.

## References

  * [1] Arbel, M., Matthews, A.G., Doucet, A.: Annealed flow transport Monte Carlo. arXiv preprint arXiv:2102.07501 (2021)
  * [2] Behrens, G., Friel, N., Hurn, M.: Tuning tempered transitions. Statistics and computing 22(1), 65–78 (2012)
  * [3] Blei, D.M., Kucukelbir, A., McAuliffe, J.D.: Variational inference: A review for statisticians. Journal of the American statistical Association 112(518), 859–877 (2017)
  * [4] Brooks, S., Gelman, A., Jones, G., Meng, X.L.: Handbook of Markov chain Monte Carlo. CRC press (2011)
  * [5] Burda, Y., Grosse, R., Salakhutdinov, R.: Importance weighted autoencoders. arXiv preprint arXiv:1509.00519 (2015)
  * [6] Calderhead, B., Girolami, M.: Estimating bayes factors via thermodynamic integration and population mcmc. Computational Statistics & Data Analysis 53(12), 4028–4045 (2009)
  * [7] Caterini, A.L., Doucet, A., Sejdinovic, D.: Hamiltonian variational auto-encoder. arXiv preprint arXiv:1805.11328 (2018)
  * [8] Del Moral, P., Doucet, A., Jasra, A.: Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(3), 411–436 (2006)
  * [9] Domke, J., Sheldon, D.: Importance weighting and variational inference. arXiv preprint arXiv:1808.09034 (2018)
  * [10] Donahue, J., Krähenbühl, P., Darrell, T.: Adversarial feature learning. arXiv preprint arXiv:1605.09782 (2016)
  * [11] Elvira, V., Martino, L., Robert, C.P.: Rethinking the effective sample size. arXiv preprint arXiv:1809.04129 (2018)
  * [12] Geffner, T., Domke, J.: MCMC variational inference via uncorrected Hamiltonian annealing. arXiv preprint arXiv:2107.04150 (2021)
  * [13] Gelman, A., Meng, X.L.: Simulating normalizing constants: From importance sampling to bridge sampling to path sampling. Statistical science pp. 163–185 (1998)
  * [14] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. Advances in neural information processing systems 27 (2014)
  * [15] Grosse, R.B., Ancha, S., Roy, D.M.: Measuring the reliability of mcmc inference with bidirectional Monte Carlo. arXiv preprint arXiv:1606.02275 (2016)
  * [16] Grosse, R.B., Ghahramani, Z., Adams, R.P.: Sandwiching the marginal likelihood using bidirectional Monte Carlo. CoRR abs/1511.02543 (2015), http://arxiv.org/abs/1511.02543
  * [17] Grosse, R.B., Maddison, C.J., Salakhutdinov, R.: Annealing between distributions by averaging moments. In: NIPS. pp. 2769–2777. Citeseer (2013)
  * [18] Gu, S., Ghahramani, Z., Turner, R.E.: Neural adaptive sequential Monte Carlo. arXiv preprint arXiv:1506.03338 (2015)
  * [19] Hoffman, M.D.: Learning deep latent gaussian models with Markov chain Monte Carlo. In: International Conference on Machine Learning. pp. 1510–1519. PMLR (2017)
  * [20] Huang, S., Makhzani, A., Cao, Y., Grosse, R.: Evaluating lossy compression rates of deep generative models. In: International Conference on Machine Learning. pp. 4444–4454. PMLR (2020)
  * [21] Jasra, A., Stephens, D.A., Doucet, A., Tsagaris, T.: Inference for lévy-driven stochastic volatility models via adaptive sequential Monte Carlo. Scandinavian Journal of Statistics 38(1), 1–22 (2011)
  * [22] Johansen, A.M., Aston, J.A., Zhou, Y.: Towards automatic model comparison: An adaptive sequential Monte Carlo approach (2015)
  * [23] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. In: Bengio, Y., LeCun, Y. (eds.) 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015), http://arxiv.org/abs/1412.6980
  * [24] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
  * [25] Kingma, D.P., Salimans, T., Jozefowicz, R., Chen, X., Sutskever, I., Welling, M.: Improved variational inference with inverse autoregressive flow. In: Advances in Neural Information Processing Systems. vol. 29 (2016)
  * [26] Kiwaki, T.: Variational optimization of annealing schedules. arXiv preprint arXiv:1502.05313 (2015)
  * [27] Laubenfels, R.: Feynman–kac formulae: genealogical and interacting particle systems with applications (2005)
  * [28] Maddison, C.J., Lawson, D., Tucker, G., Heess, N., Norouzi, M., Mnih, A., Doucet, A., Teh, Y.W.: Filtering variational objectives. arXiv preprint arXiv:1705.09279 (2017)
  * [29] Masrani, V., Brekelmans, R., Bui, T., Nielsen, F., Galstyan, A., Steeg, G.V., Wood, F.: q-paths: Generalizing the geometric annealing path using power means. arXiv preprint arXiv:2107.00745 (2021)
  * [30] Mnih, A., Rezende, D.: Variational inference for Monte Carlo objectives. In: International Conference on Machine Learning. pp. 2188–2196. PMLR (2016)
  * [31] Naesseth, C., Linderman, S., Ranganath, R., Blei, D.: Variational sequential Monte Carlo. In: International Conference on Artificial Intelligence and Statistics. pp. 968–977. PMLR (2018)
  * [32] Naesseth, C., Lindsten, F., Schon, T.: Nested sequential Monte Carlo methods. In: International Conference on Machine Learning. pp. 1292–1301. PMLR (2015)
  * [33] Naesseth, C., Ruiz, F., Linderman, S., Blei, D.: Reparameterization gradients through acceptance-rejection sampling algorithms. In: Artificial Intelligence and Statistics. pp. 489–498. PMLR (2017)
  * [34] Neal, R.M.: Annealed importance sampling. Statistics and computing 11(2), 125–139 (2001)
  * [35] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch (2017)
  * [36] Ranganath, R., Gerrish, S., Blei, D.: Black box variational inference. In: Artificial intelligence and statistics. pp. 814–822. PMLR (2014)
  * [37] Rezende, D., Mohamed, S.: Variational inference with normalizing flows. In: International Conference on Machine Learning. pp. 1530–1538. PMLR (2015)
  * [38] Rezende, D.J., Mohamed, S., Wierstra, D.: Stochastic backpropagation and approximate inference in deep generative models. In: International conference on machine learning. pp. 1278–1286. PMLR (2014)
  * [39] Sajjadi, M.S., Bachem, O., Lucic, M., Bousquet, O., Gelly, S.: Assessing generative models via precision and recall. Advances in Neural Information Processing Systems 31 (2018)
  * [40] Salakhutdinov, R., Murray, I.: On the quantitative analysis of deep belief networks. In: Proceedings of the 25th international conference on Machine learning. pp. 872–879 (2008)
  * [41] Salimans, T., Goodfellow, I., Zaremba, W., Cheung, V., Radford, A., Chen, X.: Improved techniques for training gans. Advances in neural information processing systems 29 (2016)
  * [42] Salimans, T., Kingma, D., Welling, M.: Markov chain Monte Carlo and variational inference: Bridging the gap. In: International Conference on Machine Learning. pp. 1218–1226. PMLR (2015)
  * [43] Skilling, J., et al.: Nested sampling for general bayesian computation. Bayesian analysis 1(4), 833–859 (2006)
  * [44] Theis, L., Oord, A.v.d., Bethge, M.: A note on the evaluation of generative models. arXiv preprint arXiv:1511.01844 (2015)
  * [45] Thin, A., Kotelevskii, N., Denain, J.S., Grinsztajn, L., Durmus, A., Panov, M., Moulines, E.: Metflow: A new efficient method for bridging the gap between Markov chain Monte Carlo and variational inference. arXiv preprint arXiv:2002.12253 (2020)
  * [46] Thin, A., Kotelevskii, N., Doucet, A., Durmus, A., Moulines, E., Panov, M.: Monte Carlo variational auto-encoders. In: International Conference on Machine Learning. pp. 10247–10257. PMLR (2021)
  * [47] Tom, T.A.L.M.I., Wood, J.T.R.F.: Auto-encoding sequential Monte Carlo. stat 1050, 29 (2017)
  * [48] Wolf, C., Karl, M., van der Smagt, P.: Variational inference with Hamiltonian Monte Carlo. arXiv preprint arXiv:1609.08203 (2016)
  * [49] Wu, H., Köhler, J., Noé, F.: Stochastic normalizing flows. arXiv preprint arXiv:2002.06707 (2020)
  * [50] Wu, Y., Burda, Y., Salakhutdinov, R., Grosse, R.: On the quantitative analysis of decoder-based generative models. CoRR arXiv preprint arXiv:1611.04273 (2016)
  * [51] Yu, Y., Zhang, W., Deng, Y.: Frechet inception distance (fid) for evaluating gans
  * [52] Zhang, G., Hsu, K., Li, J., Finn, C., Grosse, R.: Differentiable annealed importance sampling and the perils of gradient noise. arXiv preprint arXiv:2107.10211 (2021)


11institutetext:  King’s College London, School of Biomedical Engineering &
Imaging Sciences, London, SE1 7EU, UK 11email<EMAIL_ADDRESS>
22institutetext: Ultromics Ltd., Oxford, OX4 2SU, UK 22email:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
# Efficient Pix2Vox++ for 3D Cardiac Reconstruction from 2D echo views

David Stojanovski  This work was supported by the Wellcome/EPSRC Centre for
Medical Engineering [WT203148/Z/16/Z], by the British Heart Foundation
[TG/17/3/33406] and the National Institute for Health Research (NIHR)
Biomedical Research Centre at Guy’s and St Thomas’ NHS Foundation Trust and
King’s College London. Pablo Lamata holds a Wellcome Trust Senior Research
Fellowship [209450/Z/17/Z]. The views expressed are those of the author(s) and
not necessarily those of the NHS, the NIHR or the Department of Health. 11
Uxio Hermida 11 Marica Muffoletto 11 Pablo Lamata 11 Arian Beqiri 1122 Alberto
Gomez 1122

###### Abstract

Accurate geometric quantification of the human heart is a key step in the
diagnosis of numerous cardiac diseases, and in the management of cardiac
patients. Ultrasound imaging is the primary modality for cardiac imaging,
however acquisition requires high operator skill, and its interpretation and
analysis is difficult due to artifacts. Reconstructing cardiac anatomy in 3D
can enable discovery of new biomarkers and make imaging less dependent on
operator expertise, however most ultrasound systems only have 2D imaging
capabilities. We propose both a simple alteration to the Pix2Vox++ networks
for a sizeable reduction in memory usage and computational complexity, and a
pipeline to perform reconstruction of 3D anatomy from 2D standard cardiac
views, effectively enabling 3D anatomical reconstruction from limited 2D data.
We evaluate our pipeline using synthetically generated data achieving accurate
3D whole-heart reconstructions (peak intersection over union score $>0.88$)
from just two standard anatomical 2D views of the heart. We also show
preliminary results using real echo images.

###### Keywords:

Convolutional Neural Networks 2D to 3D reconstruction Deep learning Ultrasound

##  1 Introduction

###  1.1 Motivation and background

The most common imaging modality for cardiac assessment worldwide is
ultrasound (US) [2] because it is more affordable and safer, and it has higher
temporal resolution when compared to other modalities. The reliance on
operator skill to acquire high quality standard views is the main limitation
of US, particularly for cardiac applications (echocardiography).

Most US systems only have 2D capabilities, and most echo protocols are limited
to 2D modes [13]. 3D data is desirable both for anatomical and functional
understanding of the heart. The heart is a complex 3D structure, and
subsequently most features cannot be fully captured within 2D planes.
Ventricular and atrial walls are curved surfaces, valve hinges do not sit on a
2D circle but rather a 3D saddle shape, and as a result current protocols
require many standard planes to be able to assess most features, while a high
quality 3D reconstruction would capture them all. Cardiac function is
quantified primarily via tissue motion and blood flow. Blood quantity is
volumetric, but usually 2D approximations are used, e.g. 2D left ventricular
delineation on a 4 chamber view to compute ejection fraction. This has been
commonly used both in classical clinical disease detection and deep learning
methods [15, 10]. Tissue motion is analysed using 2D components of the motion
e.g. longitudinal, radial, etc., however cardiac motion is complex and
combines torsion, vertical motion and displacement, which are fundamentally 3D
phenomena. Capturing only 2D views can impede further analysis if the quality
of the views is suboptimal from factors such as foreshortening and artifacts.

3D volumes would allow for more accurate quantification of the cardiac
geometries, and in turn allow clinicians to reduce the rate of disease
misdiagnosis as a result of the incorrect quantification arising from using 2D
data. 3D ultrasound probes are available for clinical use. They are however
very prone to artifacting, particularly in moving anatomies such as the heart,
and as a result are very rarely used in practice [9].

The ability to quantify cardiac function using standard 2D US systems, without
the need for external trackers to define a spatial geometry, can be of great
clinical use, but is limited by lack of paired US and native 3D ground truth
anatomy. The aim of this work is to develop a method for reconstructing a full
3D representation of the heart from 2D image views.

###  1.2 Related works

The advent of Deep Learning (DL) has led to development of 3D image
reconstruction methods that have been shown to successfully learn shape and
structure from partial observations, achieving start-of-the-art (SOTA)
performance in natural and medical imaging.

####  1.2.1 2D to 3D reconstruction in natural images

Current SOTA methods using DL allow accurate 2D to 3D reconstruction of
objects from one or more RGB images, without the aid of camera positioning
calibration information. A popular standardized dataset used to perform
comparisons of reconstruction techniques is ShapeNetCore [4]. ShapeNetCore
covers 55 object categories with over 51,300 unique 3D models. In 2019,
ShapeNetCore reconstruction SOTA accuracy was achieved by Xie et al. who were
able to reach 0.706 Intersection over Union (IoU) using the Pix2Vox++ (PiVox)
network [16]. Two network variations, Pix2Vox++/Fast (PiVox/Fast) and
Pix2Vox++/Accurate, or (PiVox/Accurate), were proposed with the PiVox/Fast
network being a lighter weight, albeit lower performing, version of the full
PiVox/Accurate network. However, PiVox/Fast is still an expensive model to
train and perform inference, thus severely limiting the possible
reconstruction resolution on currently available single Graphics Processing
Units (GPUs).

####  1.2.2 Ultrasound specific 2D to 3D reconstruction

A number of US-specific 3D reconstruction algorithms have been developed
showing effective performance but have not been applied to cardiac
reconstruction. Cerrolaza et al. [1] performed fetal skull reconstruction
using conditional hierarchical generative network Variational Auto-Encoders
(VAEs) and achieved a Dice Coefficient (DC) of 0.91 using three orthogonal US
views. In contrast to our goal of performing cardiac reconstructions, their
reconstruction target had a generally regular shape which does not deform
during imaging.

Prevost et al. [12] combine DL with an Inertial Measurement Unit (IMU), using
two consecutive frames with optical flow as a channel to the network along
with the Euler angles provided by the IMU. On a phantom they achieved a
minimum, medium and maximum drift of 1.70, 18.30 and 36.90 mm respectively.

###  1.3 Contributions of this study

No previous work has addressed 3D cardiac shape reconstructions from untracked
2D echo views. As a first step towards this goal, we present an exploratory
investigation using synthetic echo data: 1) A pipeline to synthesize 2D echo
views with corresponding 3D ground truth; 2) A demonstration that PiVox can
produce accurate 2D to 3D reconstructions of realistic synthetic hearts; 3) A
simple modification to both PiVox networks which vastly reduces memory and
computational expense, while still achieving high reconstruction accuracy.

##  2 Methods

To explore reconstruction of a 3D heart from 2D data, we used two types of
synthetic data: 1) A segmentation dataset, containing binary tissue masks
(segmentations) on 2D standard views. Masks are simulated by slicing 3D
computational models (which allows using 3D ground truths for evaluation); 2)
A synthetic US dataset, containing synthetic 2D standard US views, generated
from the tissue masks. We used both for model training and testing. An
overview of the proposed pipeline is shown in Figure 1.

Figure 1: Proposed pipeline for extracting 2D standard US views from a mesh
and performing a reconstruction. A) Initial heart model with tissue labels, B)
Mesh slicing with a plane and extracted segmentation-like images and C) Input
images being trained/inferenced in the PiVox networks to produce 3D
reconstructions.

###  2.1 Synthetic 2D segmentations from 3D heart models

We obtain synthetic 2D segmentations by slicing 3D mesh models of the heart at
standard echocardiographic views. These will be used to asses both how the 2D
to 3D reconstruction networks are affected by the large amounts of variation
that will be present in real ultrasound images, and also to provide insight
into the feasibility of performing training on synthetic data and testing on
real ultrasound data. If it can be shown that training using just synthetic
data and testing on real data is feasible then it would allow generation of
large training datasets and greatly reduce the amount of paired 3D CT/2D
ultrasound data required to validate methodologies.

Our slice extraction technique utilizes a combination of the Visualization
Toolkit (VTK) [8] and the PyVista Python packages [3]. The mesh data used was
a set of 1000 synthetic cardiac meshes created by Rodero et al. [14].

The first step in extracting the standard plane cardiac views requires
defining either 3 landmarks for each view to which a plane can be fit, or 2
landmarks and a projection along a vector. We defined these automatically-
calculated landmarks based on their relation to various cardiac structures
present in the meshes/segmentations (cf. Table 1). A trained sonographer
examined 10 sets of segmentations to confirm the chosen landmarks resulted in
suitably realistic extracted slices from our synthetic meshes.

The left ventricular apex (LVA) was found by performing a ray casting from the
centre of mass of the mitral valve to all mesh cell faces in the LV mesh. If
there was a pair of intersection points along the ray line (i.e. the ray
penetrates both the endo- and epicardial wall) the distance between these
intersection points was calculated. The shortest distance was chosen as the
LVA, to find the thinnest location of the cardiac wall and to minimize apical
foreshortening.

Table 1: Ultrasound views and the corresponding landmarks to define these
planes.

Ultrasound View Landmarks Right Ventricular Inflow Right Ventricle, Right
Atrium, Pulmonary Valve Left Ventricular Parasternal Long-Axis Left Ventricle
Apex, Mitral Valve, Aortic Valve Parasternal Short-Axis Aortic Valve Level
Left Atrium, Right Atrium, Aortic Valve Parasternal Short-Axis Mitral Valve
Level Left Ventricle Apex, Whole Heart Parasternal Short-Axis Papillary Muscle
Level Left Ventricle Apex, Whole Heart Parasternal Short-Axis Apex Level Left
Ventricle Apex, Whole Heart Apical 4 Chamber Left Atrium, Right Atrium, Left
Ventricle Apex Apical 5 Chamber Left Atrium, Aortic Valve, Left Ventricle Apex
Apical 2 Chamber Left Ventricle Apex, Mitral Valve, Right Ventricle

###  2.2 Synthetic Ultrasound images

Synthetic Ultrasound images (i.e. with real appearance) are generated from the
same 3D heart models to study feasibility in realistic echo data while having
an exact ground truth. We adapted the technique proposed by Gilbert et al. [5]
using the CAMUS dataset [7]. In brief, pseudo-images, i.e. tissue masks
corresponding to 2 and 4 chamber images, with noise and Gaussian blurring are
used as input to a CycleGAN network [17] with unpaired US images to create the
final synthetic images, as exemplified in Fig. 2.

Figure 2: Example case of data generated at each stage in the process of
creating a synthetic ultrasound image from a segmentation mask. From left to
right, a) Initial segmentation image b) generated pseudo image c) generated
synthetic US d) overlay of input mask on generated image.

Producing high quality synthetic US images required small variations in 1)
noise parameters 2) sequence of additive noise and blurring, and 3) size of
the Gaussian blur kernel. These parameters were changed due to a different
resolution input image, and our requirement to not perform any geometric
transforms, e.g. changing of input anatomy dimensions.

###  2.3 Efficient Pix2Vox++

The general Pix2Vox++ (PiVox) network architecture is composed of a series of
parallel encoder and decoder branches for each input view, before being passed
into the fusion and refiner modules. The decoder, in particular, is memory and
computationally expensive due to the 3D deconvolutional kernels. This greatly
limits the resolution of the 3D volumes that can be used, thus necessitating a
more efficient network with high reconstruction accuracy. It also contains a
very large number of learnable parameters, in turn requiring a very large
training set.

Xie et al. proposed a lighter weight, albeit lower performing, variant of the
PiVox/Accurate network, referred to as PiVox/Fast, for reducing memory and
computational complexity by 1) Using ResNet-18 instead of ResNet-50 2)
Decreased de/convolution kernel sizes 3) Removal of the refiner module.

We propose a simple adaptation of the PiVox networks, referred to as E-PiVox,
to greatly decrease both memory usage and computational expense, with minimal
impact to performance. We achieve this by adding a 3D convolution at the end
of the encoder module along the input image dimension to reduce the decoder to
a single branch for any number of input images. All information is propagated
through the parallel decoder branches in the PiVox networks passes through a
compressed latent space, in theory allowing for minimal decreases in
performance from having just a single decoder branch rather than multiple
branches.

###  2.4 2D to 3D reconstruction of synthetic hearts

####  2.4.1 Implementation details

The dataset of 1000 cardiac meshes was divided into a 70/15/15$\%$
train/validation/test split and remained unchanged for each of the different
training runs. The models were trained with binary segmentation masks, the
label maps converted to realistic echo-like images using the CycleGAN network
and also on the ShapeNet dataset. All training was performed using Pytorch
1.9.1 [11] on an Nvidia RTX 3090 for 200 epochs. The Adam optimizer was used
with a $\beta_{1}$ of 0.9 and a $\beta_{2}$ of 0.999 [6]. The code is
available at https://github.com/david-stojanovski/E-Pix2Vox-reconstruction

##  3 Experiments and Results

Reconstruction accuracy was assessed using the thresholded Intersection over
Union (IoU) [16].

The final reconstruction results of the ShapeNet and cardiac training runs are
shown in Tables 2 and 3 respectively. Table 2 shows PiVox/Accurate was
consistently the best performing network, however E-PiVox/Fast and
E-PiVox/Accurate were consistently within $1.2\%$ and $2\%$ respectively of
PiVox/Fast and PiVox/Accurate. This was achieved while being much more memory
and computationally efficient.

Table 2: IoU Comparison of multi-view 2D to 3D object reconstruction on the ShapeNet data set for a $32^{3}$ resolution. Method | 1 View | 2 Views | 4 Views | 8 Views
---|---|---|---|---
PiVox/Fast | 0.645 | 0.669 | 0.682 | 0.690
E-PiVox/Fast | 0.641 | 0.663 | 0.675 | 0.682
Percentage diff.(%) | -0.620 | -0.897 | -1.026 | -1.16
PiVox/Accurate | 0.670 | 0.695 | 0.708 | 0.715
E-PiVox/Accurate | 0.6687 | 0.6836 | 0.6949 | 0.701
Percentage diff.(%) | -0.19 | -1.64 | -1.85 | -1.96

It can be seen in Table 2 that our E-PiVox networks was able to come very
close to the performance of the much more computationally expensive PiVox
networks. The relationship between computational expense and memory usage in
regard to number of input views is shown in Figure 3 and highlights the
importance of efficient network architectures when using a larger number of
views, such as in video recordings, or higher resolutions e.g. full resolution
CT.

Table 3: IoU Comparison of multi-view 2D to 3D object reconstruction on the cardiac data set for a $64^{3}$ resolution. A2C and A4C refer to atrial 2 and 4 chamber image views respectively, and 9 views denotes all 9 of the previously described standard ultrasound views. | Binary | CycleGAN inferenced
---|---|---
method | 9 views | A2C | A4C | A2C & A4C | A2C | A4C | A2C & A4C
PiVox/Fast | 0.859 | 0.828 | 0.841 | 0.818 | 0.644 | 0.652 | 0.7
E-PiVox/Fast | 0.86 | 0.837 | 0.85 | 0.822 | 0.651 | 0.66 | 0.712
PiVox/Accurate | 0.895 | 0.85 | 0.867 | 0.873 | 0.673 | 0.732 | 0.741
E-PiVox/Accurate | 0.903 | 0.868 | 0.878 | 0.881 | 0.674 | 0.734 | 0.74

Table 3 shows the results when using both binary and realistic ultrasound
CycleGAN inferenced images. As shown in Table 3 both E-PiVox networks
generally outperformed the PiVox networks in the relative comparisons.The
reconstruction accuracy dropped most in areas where the walls of the heart
were quite thin, generally in the atrial region. Small errors appearing over
the entire structure appear to arise due to errors in the exact location
discretization of the voxels.

Figure 3: Comparison of the 4 networks and their theoretical total
Multiply/Add operations for each number of input views at a $64^{3}$
resolution.

We present examples of 3D heart reconstructions in Figure 4, including
preliminary results using real 2D echo images from the CAMUS dataset.
Synthetically trained reconstructions show close correspondence with the
ground truth anatomy (black wireframe). However, reduced number of input views
resulted in decreased accuracy and clear non-physiological holes. Such
phenomena were exacerbated when using US-like images (see Figure 4C). Larger
differences were observed in the LVA and valve planes for all models.

Figure 4: A) 9 binary views; B) Binary Apical 2 and 4 chamber views; C)
CycleGAN inferenced Apical 2 and 4 chamber views; D) Example case from the
CAMUS dataset. The same case is shown in A, B and C. The E-PiVox output column
shows the ground truth (White) and the network predictions (Orange). The 3D
Heart Anatomy column shows the absolute error in relation to the ground truth
mesh (Black wireframe).

##  4 Discussion and Conclusions

This work provides a proof of concept that a complex geometry like the human
heart can be reconstructed with reasonable accuracy from a limited number of
standard 2D anatomical views. And that this can be achieved with efficient
training and inference.

Results in Table 3 show that the PiVox network can successfully reconstruct a
full heart with a peak IoU of 0.903 for 9 segmented views, and a minimal
decrease in performance to 0.881 for just 2 views. The 0.741 IoU achieved with
just 2 synthetic US images represents an acceptable performance given the much
more challenging input image data that emulates real-world acquisitions.

It is important to note that real US images present great variability in
contrast and appearance. The preliminary results using real data (panel D in
Figure 4) show reconstructions are possible, albeit needing further
improvement. Given a suitably large enough set of real US images, accurate 3D
reconstructions could be possible from a small number of standard clinical 2D
views.

The keystone for this contribution is the synthesis of training data (i.e. 2D
views) with idealized ground truth 3D reconstructions. This approach addresses
the limitations in medical data availability, data privacy, and cost of expert
annotation. Preliminary results show promise in the use on real data using
models trained with synthetic data, but further research and evidence is
needed.

The main limitation of our study is that it is based on the anatomical
variability of a synthetic cohort of healthy subjects. As such, the
performance could dramatically drop in the presence of disease.

In conclusion, this work demonstrates, in a synthetic workbench, the
feasibility of 3D cardiac reconstruction from standard 2D views.

## References

  * [1] B, J.J.C., Li, Y., Biffi, C., Gomez, A., Sinclair, M., Matthew, J., Knight, C., Kainz, B., Rueckert, D.: 3D Fetal Skull Reconstruction from 2DUS via Deep Conditional Generative Networks, vol. 1. Springer International Publishing (2018). https://doi.org/10.1007/978-3-030-00928-1, http://dx.doi.org/10.1007/978-3-030-00928-1˙44
  * [2] Braga, J.R., Leong-Poi, H., Rac, V.E., Austin, P.C., Ross, H.J., Lee, D.S.: Trends in the Use of Cardiac Imaging for Patients with Heart Failure in Canada. JAMA Network Open 2(8), 1–13 (2019). https://doi.org/10.1001/jamanetworkopen.2019.8766
  * [3] Castro, D.d.l.I., Mader, K., Hanchen, Pandu-rao, Sullivan, B., Cole, J., Wallbaum, M., Mitchell, N., Sebastian, Movva, M., McCartney, P., Jadhav, S., de Miguel, Y.B., Milan, A., Mitzkus, B., Bazazian, D., Gandenberger, G., Buchanan, J., Trevisani, M.D., Quach, M., N-McA, Joshi, R.P., Tweedie, R., Koonce, B., Iindovina, Hadfield, J., Joskaaaa: daavoo/pyntcloud: v0.1.6 (jan 2022). https://doi.org/10.5281/ZENODO.5841822, https://zenodo.org/record/5841822
  * [4] Chang, A.X., Funkhouser, T., Guibas, L., Hanrahan, P., Huang, Q., Li, Z., Savarese, S., Savva, M., Song, S., Su, H., Xiao, J., Yi, L., Yu, F.: ShapeNet: An Information-Rich 3D Model Repository (2015), http://arxiv.org/abs/1512.03012
  * [5] Gilbert, A., Marciniak, M., Rodero, C., Lamata, P., Samset, E., McLeod, K.: Generating Synthetic Labeled Data from Existing Anatomical Models: An Example with Echocardiography Segmentation. IEEE Transactions on Medical Imaging 40(10), 2783–2794 (2021). https://doi.org/10.1109/TMI.2021.3051806
  * [6] Kingma, D.P., Ba, J.L.: Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings pp. 1–15 (2015)
  * [7] Leclerc, S., Smistad, E., Pedrosa, J., Ostvik, A., Cervenansky, F., Espinosa, F., Espeland, T., Berg, E.A.R., Jodoin, P.M., Grenier, T., Lartizien, C., Dhooge, J., Lovstakken, L., Bernard, O.: Deep Learning for Segmentation Using an Open Large-Scale Dataset in 2D Echocardiography. IEEE transactions on medical imaging 38(9), 2198–2210 (2019). https://doi.org/10.1109/TMI.2019.2900516
  * [8] Lowekamp, B.C., Chen, D.T., Ibáñez, L., Blezek, D.: The design of simpleITK. Frontiers in Neuroinformatics 7(DEC), 1–14 (2013). https://doi.org/10.3389/fninf.2013.00045
  * [9] Nelson, T.R., Pretorius, D.H., Hull, A., Riccabona, M., Sklansky, M.S., James, G.: Sources and impact of artifacts on clinical three-dimensional ultrasound imaging. Ultrasound in Obstetrics and Gynecology 16(4), 374–383 (2000). https://doi.org/10.1046/j.1469-0705.2000.00180.x
  * [10] Ouyang, D., He, B., Ghorbani, A., Yuan, N., Ebinger, J., Langlotz, C.P., Heidenreich, P.A., Harrington, R.A., Liang, D.H., Ashley, E.A., Zou, J.Y.: Video-based AI for beat-to-beat assessment of cardiac function. Nature 580(7802), 252–256 (2020). https://doi.org/10.1038/s41586-020-2145-8, http://dx.doi.org/10.1038/s41586-020-2145-8
  * [11] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.: PyTorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems 32(NeurIPS) (2019)
  * [12] Prevost, R., Salehi, M., Jagoda, S., Kumar, N., Sprung, J., Ladikos, A., Bauer, R., Zettinig, O., Wein, W.: 3D freehand ultrasound without external tracking using deep learning. Medical Image Analysis 48, 187–202 (2018). https://doi.org/10.1016/j.media.2018.06.003, https://doi.org/10.1016/j.media.2018.06.003
  * [13] Robinson, S., Rana, B., Oxborough, D., Steeds, R., Monaghan, M., Stout, M., Pearce, K., Harkness, A., Ring, L., Paton, M., Akhtar, W., Bedair, R., Bhattacharyya, S., Collins, K., Oxley, C., Sandoval, J., Schofield, R., Siva, A., Parker, K., Willis, J., Augustine, D.X.: A practical guideline for performing a comprehensive transthoracic echocardiogram in adults: The British Society of Echocardiography minimum dataset. Echo Research and Practice 7(4), G59–G93 (2020). https://doi.org/10.1530/ERP-20-0026
  * [14] Rodero, C., Strocchi, M., Marciniak, M., Longobardi, S., Whitaker, J., O’Neill, M.D., Gillette, K., Augustin, C., Plank, G., Vigmond, E.J., Lamata, P., Niederer, S.A.: Linking statistical shape models and simulated function in the healthy adult human heart. PLoS Computational Biology 17(4), 1–28 (2021). https://doi.org/10.1371/journal.pcbi.1008851, http://dx.doi.org/10.1371/journal.pcbi.1008851
  * [15] Upton, R., Mumith, A., Beqiri, A., Parker, A., Hawkes, W., Gao, S., Porumb, M., Sarwar, R., Marques, P., Markham, D., Kenworthy, J., O’Driscoll, J.M., Hassanali, N., Groves, K., Dockerill, C., Woodward, W., Alsharqi, M., McCourt, A., Wilkes, E.H., Heitner, S.B., Yadava, M., Stojanovski, D., Lamata, P., Woodward, G., Leeson, P.: Automated Echocardiographic Detection of Severe Coronary Artery Disease Using Artificial Intelligence. JACC: Cardiovascular Imaging pp. 1–13 (2022). https://doi.org/10.1016/j.jcmg.2021.10.013
  * [16] Xie, H., Yao, H., Sun, X., Zhou, S., Zhang, S.: Pix2Vox: Context-aware 3D reconstruction from single and multi-view images. Proceedings of the IEEE International Conference on Computer Vision 2019-Octob, 2690–2698 (2019). https://doi.org/10.1109/ICCV.2019.00278
  * [17] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks (2017). https://doi.org/10.48550/ARXIV.1703.10593, https://arxiv.org/abs/1703.10593