# Improving Variational Autoencoder Estimation from
Incomplete Data with Mixture Variational Families

Vaidotas Simkus<EMAIL_ADDRESS>
Michael U. Gutmann<EMAIL_ADDRESS>
School of Informatics
University of Edinburgh

###### Abstract

We consider the task of estimating variational autoencoders (VAEs) when the
training data is incomplete. We show that missing data increases the
complexity of the model’s posterior distribution over the latent variables
compared to the fully-observed case. The increased complexity may adversely
affect the fit of the model due to a mismatch between the variational and
model posterior distributions. We introduce two strategies based on (i) finite
variational-mixture and (ii) imputation-based variational-mixture
distributions to address the increased posterior complexity. Through a
comprehensive evaluation of the proposed approaches, we show that variational
mixtures are effective at improving the accuracy of VAE estimation from
incomplete data.

##  1 Introduction

Deep latent variable models, as introduced by Kingma & Welling (2013); Rezende
et al. (2014); Goodfellow et al. (2014); Sohl-Dickstein et al. (2015);
Krishnan et al. (2016); Dinh et al. (2017), have emerged as a predominant
approach to model real-world data. The models excel in capturing the intricate
nature of data by representing it within a well-structured latent space.
_However, they typically require large amounts of fully-observed data at
training time, while practitioners in many domains often only have access to
incomplete data sets_.

In this paper we focus on the class of variational autoencoders (VAEs, Kingma
& Welling, 2013; Rezende et al., 2014) and investigate the implications of
incomplete training data on model estimation. Our contributions are as
follows:

  * •

We show that data missingness can add significant complexity to the model
posterior of the latent variables, hence requiring more flexible variational
families compared to scenarios with fully-observed data (section 3).

  * •

We propose finite variational-mixture approaches to deal with the increased
complexity due to missingness for both standard and importance-weighted ELBOs
(section 4.1).

  * •

We further propose an imputation-based variational-mixture approach, which
decouples model estimation from data missingness problems, and as a result,
improves model estimation when using the standard ELBO (section 4.2).

  * •

We evaluate the proposed methods for VAE estimation on synthetic and realistic
data sets with missing data (section 6).

The proposed methods achieve better or similar estimation performance compared
to existing methods that do not use variational mixtures. Moreover, the
mixtures are formed by the variational families that are used in the fully-
observed case, which allows us to seamlessly re-use the inductive biases from
the well-studied scenarios with fully-observed data (see e.g. Miao et al.,
2022, for the importance of inductive biases in VAEs).

##  2 Background: Standard approach for VAEs estimation from incomplete data

We consider the situation where some part of the training data-points might be
missing. We denote the observed and missing parts of the $i$-th data-point
${\bm{x}}^{i}$ by ${\bm{x}}_{\mathrm{obs}}^{i}$ and
${\bm{x}}_{\mathrm{mis}}^{i}$, respectively. This split into observed and
missing components corresponds to a missingness pattern
${\bm{m}}^{i}\in\\{0,1\\}^{D}$, which can be different for each $i$, and is
itself a random variable that follows a typically unknown missingness
distribution ${p^{*}}({\bm{m}}^{i}\mid{\bm{x}}^{i})$. We make the common
assumption that the missingness distribution does not depend on the missing
variables ${\bm{x}}_{\mathrm{mis}}^{i}$, that is,
${p^{*}}({\bm{m}}^{i}\mid{\bm{x}}^{i})={p^{*}}({\bm{m}}^{i}\mid{\bm{x}}_{\mathrm{obs}}^{i})$,
which is known as the ignorable missingness or missing-at-random assumption
(MAR, e.g. Little & Rubin, 2002, Section 1.3). The MAR assumption allows us to
ignore the missingness pattern ${\bm{m}}^{i}$ when fitting a model
${p_{\bm{\theta}}}({\bm{x}})$ of the true distribution ${p^{*}}({\bm{x}})$
from incomplete data.

The VAE model with parameters ${\bm{\theta}}$ is typically specified using a
decoder distribution ${p_{\bm{\theta}}}({\bm{x}}\mid{\bm{z}})$, parametrised
using a neural network, and a prior ${p_{\bm{\theta}}}({\bm{z}})$ over the
latents ${\bm{z}}$ that can either be fixed or learnt. A principled approach
to handling incomplete training data is then to marginalise the missing
variables from the likelihood ${p_{\bm{\theta}}}({\bm{x}})$, which yields the
marginal likelihood

$\displaystyle{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}}^{i})$
$\displaystyle=\int{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{i})\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}^{i}=\iint{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{i}\mid{\bm{z}}){p_{\bm{\theta}}}({\bm{z}})\mathop{}\\!\mathrm{d}{\bm{z}}\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}^{i}=\int{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}}^{i}\mid{\bm{z}}){p_{\bm{\theta}}}({\bm{z}})\mathop{}\\!\mathrm{d}{\bm{z}},$
(1)

where the inner integral
$\int{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}}\mid{\bm{z}})\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}$
is often computationally tractable in VAEs due to standard assumptions, such
as the conditional independence of ${\bm{x}}$ given ${\bm{z}}$ or the use of
the Gaussian family for the decoder ${p_{\bm{\theta}}}({\bm{x}}\mid{\bm{z}})$.
However, the marginal likelihood above remains intractable to compute as a
consequence of the integral over the latents ${\bm{z}}$.

Due to the intractable integral, VAEs are typically fitted via a variational
evidence lower-bound (ELBO)

$\displaystyle\log{p_{\bm{\theta}}}({\bm{y}})$
$\displaystyle\geq\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{y}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{y}}\mid{\bm{z}}){p_{\bm{\theta}}}({\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{y}})}\right]=\log{p_{\bm{\theta}}}({\bm{y}})-D_{\text{KL}}({q_{\bm{\phi}}}({\bm{z}}\mid{\bm{y}})\;||\;{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{y}})),$
(2)

where ${\bm{y}}$ refers to ${\bm{x}}^{i}$ in the fully-observed case, and to
${\bm{x}}_{\mathrm{obs}}^{i}$ in the incomplete-data case, and
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{y}})$ is an (amortised) variational
distribution with parameters ${\bm{\phi}}$ that is shared for all data-points
in the data set (Gershman & Goodman, 2014). The amortised distribution is
parametrised using a neural network (the encoder), which takes the data-point
${\bm{y}}$ as the input and predicts the distributional parameters of the
variational family. Moreover, when the data is incomplete, i.e.
${\bm{y}}={\bm{x}}_{\mathrm{obs}}^{i}$, sharing of the encoder for any pattern
of missingness is often achieved by fixing the input dimensionality of the
encoder to twice the size of ${\bm{x}}$ and providing
$\gamma({\bm{x}}_{\mathrm{obs}}^{i})$ and ${\bm{m}}^{i}$ as the
inputs,111Alternative encoder architectures, such as, permutation-invariant
networks (Ma et al., 2019) are also used. where $\gamma(\cdot)$ is a function
that takes the incomplete data-point ${\bm{x}}_{\mathrm{obs}}$ and produces a
vector of length $D$ with the missing dimensions set to zero222Equivalent to
setting the missing dimensions to the empirical mean for zero-centered data.
(Nazábal et al., 2020; Mattei & Frellsen, 2019).

In eq. 2 we show that the training objective for incomplete and fully-observed
data has the same form, and therefore it may seem that fitting VAEs from
incomplete data would be similarly difficult to the fully-observed case.
However, as we will see next, data missingness can make model estimation much
harder than in the complete data case.

##  3 Implications of incomplete data for VAE estimation

Figure 1: _Illustration of the posterior complexity due to missing data._ Each
colour represents a different data-point ${\bm{x}}^{i}$. First: the model
posterior ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ under complete data
${\bm{x}}$. Second: the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ under incomplete data
${\bm{x}}_{\mathrm{obs}}$. Third: variational approximation
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ of the complete-data posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$. Fourth: an imputation-mixture
variational approximation
$\mathbb{E}_{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]$
of the incomplete posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$.

The decomposition of the ELBO in eq. 2 emphasises that accurate estimation of
the VAE model requires us to accurately approximate the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ with the variational
distribution ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$. While it
might appear that the marginalisation of the missing variables in eq. 1 comes
at no cost since the ELBO maintains the same form as in the complete case, we
here illustrate that his is not the case.

In the two left-most columns of fig. 1 we illustrate the model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid\cdot)$ under fully-observed data ${\bm{x}}$
and partially-observed data ${\bm{x}}_{\mathrm{obs}}$.333In fig. 1 we use a
VAE with Gaussian variational, prior, and decoder distributions fitted on
complete data. We discover that the model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$, which exhibited a certain
regularity in the complete-data scenario, have become irregular multimodal
distributions ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ when
evaluated with incomplete data.444A related phenomenon, called posterior
inconsistency, has been recently reported in concurrent work by Sudak &
Tschiatschek (2023), relating
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}\setminus u})$, where $u$
is a subset of the observed dimensions (see section 5). Hence, accurate
estimation of VAEs from incomplete data may require more flexible variational
families than in the fully-observed case: while a Gaussian family may
sufficiently well approximate the model posterior in the fully-observed case
of our example, it is no longer sufficiently flexible in the incomplete data
case. We provide a further explanation when this situation may occur in
appendix A. As a result of the mismatch between the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and the variational
distribution ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, the KL
divergence term in eq. 2 may not be minimised, subsequently introducing a bias
to the fit of the model.

In the two right-most columns of fig. 1 we show the variational distributions
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ under fully-observed data ${\bm{x}}$
and approximations of the incomplete-data posteriors
$\mathbb{E}_{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]$,
which are good approximations of ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, respectively. The
two plots show that if the variational family used in the fully-observed case
well-approximates the model posterior, i.e.
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})\approx{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$,
then the imputation-mixture
$\mathbb{E}_{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]$
will also be a good approximation of the incomplete-data posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$. This observation
suggests that we can work with the same variational family in both the fully-
observed and incomplete data scenarios if we adopt a mixture approach. In the
rest of this paper, we investigate opportunities to improve VAE estimation
from incomplete data by constructing variational mixture approximations of the
incomplete-data posterior.

##  4 Fitting VAEs from incomplete data using mixture variational families

We propose working with mixture variational families in order to mitigate the
increase in posterior complexity due to missing data and improve the
estimation accuracy of VAEs when the training data are incomplete. This allows
us to use families of distributions for the mixture components that are known
to work well when the data is fully-observed, and use the mixtures to handle
the increased posterior complexity due to data missingness.

We propose two approaches for constructing variational mixtures. In section
4.1 we specify ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ as a
finite-mixture distribution that can be learnt directly using the
reparametrisation trick. In section 4.2 we investigate an imputation-based
variational-mixture where we specify
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\approx\mathbb{E}_{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]$.
Detailed evaluation of the proposed methods is provided in section 6.

###  4.1 Using finite mixture variational distributions to fit VAEs from
incomplete data

In section 3 we saw that a good approximation of the incomplete data posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ would be the
imputation-mixture
$\mathbb{E}_{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]$.
However, estimation of
${p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$ is
generally intractable for VAEs (Rezende et al., 2014; Mattei & Frellsen,
2018a; Simkus & Gutmann, 2023). Hence, we here consider a more tractable
approach and specify the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ in terms of a finite-
mixture distribution:

$\displaystyle{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})=\sum_{k=1}^{K}{q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}}){q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}}),$
(3)

where ${q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$ is a categorical
distribution over the components $k\in\\{1,\ldots,K\\}$ and each component
distribution ${q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
belongs to any reparametrisable distribution family. Both
${q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$ and
${q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ are amortised using
an encoder network, similar to section 2.

The “reparametrisation trick” is typically used in VAEs to efficiently
optimise the parameters ${\bm{\phi}}$ of the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, which requires that
the random variable ${\bm{z}}$ can be parametrised as a learnable
differentiable transformation
$t(\bm{\epsilon};{\bm{x}}_{\mathrm{obs}},{\bm{\phi}})$ of another random
variable $\bm{\epsilon}$ that follows a distribution with no learnable
parameters. However, reparametrising mixture-families requires extra care:
sampling the mixture ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ in
eq. 3 is typically done via ancestral sampling by first drawing
$k\sim{q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$ and then
${\bm{z}}\sim{q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, but
the sampling of the categorical distribution
${q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$ is non-differentiable, making
the direct application of the “reparametrisation trick” generally infeasible.

As a result, we consider two objectives for fitting VAEs with mixture-
variational distributions based on the variational ELBO (Kingma & Welling,
2013; Rezende et al., 2014):

$\displaystyle\mathcal{L}_{\mathrm{ELBO}}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle=\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})}\left[\log
w({\bm{z}})\right],\quad\text{and}$ (4)
$\displaystyle\mathcal{L}_{\mathrm{SELBO}}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle=\sum_{k=1}^{K}{q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})\mathbb{E}_{{q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})}\left[\log
w({\bm{z}})\right],$ (5) $\displaystyle\text{where}\quad w({\bm{z}})$
$\displaystyle=\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})}.$
(6)

The first objective $\mathcal{L}_{\mathrm{ELBO}}$ corresponds to the standard
ELBO, while $\mathcal{L}_{\mathrm{SELBO}}$ is the stratified ELBO (Roeder et
al., 2017, Section 4; Morningstar et al., 2021). When working with
$\mathcal{L}_{\mathrm{ELBO}}$, due to the mixture variational family, we will
need to optimise ${\bm{\phi}}$ with _implicit_ reparameterisation (Figurnov et
al., 2019). Implicit reparametrisation of mixture distributions requires that
it is possible to factorise the component distributions
${q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ using the chain
rule, i.e.
${q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})=\prod_{d}{q_{\bm{\phi}}^{k}}(z_{d}\mid{\bm{z}}_{<d},{\bm{x}}_{\mathrm{obs}})$,
and have access to the CDF (or other standardisation function) of each factor
${q_{\bm{\phi}}^{k}}(z_{d}\mid{\bm{z}}_{<d},{\bm{x}}_{\mathrm{obs}})$.
However, the chain rule requirement can be difficult to satisfy for some
highly flexible variational distribution families, such as normalising flows
(e.g. Papamakarios et al., 2021), and finding the (conditional) CDF of the
factors can also be hard if not already known in closed form. Consequently,
$\mathcal{L}_{\mathrm{ELBO}}$ with implicit reparametrisation may not be
usable with all families of variational distributions as components of the
mixture.

The second objective $\mathcal{L}_{\mathrm{SELBO}}$, on the other hand,
samples the mixture distribution with stratified sampling,555Stratified
sampling of mixture distributions typically draws an equal number of samples
from each component and weighs the samples by the component probabilities
${q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$ when estimating expectations.
It is commonly used to reduce Monte Carlo variance (Robert & Casella, 2004).
which avoids the non-differentiability of sampling
${q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})$, and as a result allows us to
use any family of reparametrisable distributions as the mixture components.

The importance-weighted ELBO (IWELBO, Burda et al., 2015) is often used as an
alternative to the standard ELBO as it can be made tighter. We here also
consider an ordinary version, $\mathcal{L}_{\mathrm{IWELBO}}$, and a
stratified version, $\mathcal{L}_{\mathrm{SIWELBO}}$ (Shi et al., 2019,
Appendix A; Morningstar et al., 2021):

$\displaystyle\mathcal{L}_{\mathrm{IWELBO}}^{I}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle=\mathbb{E}_{\\{{\bm{z}}_{j}\\}_{j=1}^{I}\sim{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})}\left[\log\frac{1}{I}\sum_{j=1}^{I}w({\bm{z}}_{j})\right],\quad\text{and}$
(7) $\displaystyle\mathcal{L}_{\mathrm{SIWELBO}}^{I}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle=\mathbb{E}_{\\{\\{{\bm{z}}_{j}^{k}\\}_{j=1}^{I}\sim{q_{\bm{\phi}}^{k}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\\}_{k=1}^{K}}\left[\log\sum_{k=1}^{K}{q_{\bm{\phi}}}(k\mid{\bm{x}}_{\mathrm{obs}})\frac{1}{I}\sum_{j=1}^{I}w({\bm{z}}_{j}^{k})\right],$
(8)

where $I$ is the number of importance samples in
$\mathcal{L}_{\mathrm{IWELBO}}$ and the number of samples per-mixture-
component in $\mathcal{L}_{\mathrm{SIWELBO}}$.

When the number of mixture-components is $K=1$ the lower-bounds above
correspond to the MVAE and MIWAE bounds in Mattei & Frellsen (2019) which are
among the most popular bounds for fitting VAEs from incomplete data. However,
as $K>1$ the proposed bounds can be tighter due to an increased flexibility of
the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ (Morningstar et al.,
2021, Appendix A), which potentially mitigates the problems caused by the
missing data (see section 3). Finally, the importance-weighted bounds in eqs.
7 and 8 maintain the asymptotic consistency guarantees of Burda et al. (2015)
and approaches the true marginal log-likelihood
$\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$ as $K\cdot
I\rightarrow\infty$, allowing for more accurate estimation of the model with
increasing computational budget.

We denote the four methods based on eqs. 4, 5, 7 and 8 by MissVAE, MissSVAE,
MissIWAE, and MissSIWAE respectively.

###  4.2 Using imputation-mixture distributions to fit VAEs from incomplete
data

In section 4.1, we dealt with the inference of the latents ${\bm{z}}$ (section
2) and the pitfalls of missing data (section 3) jointly by learning a finite-
mixture variational distribution. Here, we propose a second “decomposed”
approach to deal with the pitfalls of missing data.

Intuitively, if we had an oracle that was able to generate imputations of the
missing data based on the ground truth conditional distribution
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$, then the VAE
estimation task would reduce to the case of complete-data, that is, the
challenges affecting the estimation of the variational distribution
${q_{\bm{\phi}}}$ from section 3 would be mitigated. This suggests that an
_effective strategy would be to decompose the task of model estimation from
incomplete data into two (iterative) tasks: data imputation and model
estimation_ , akin to the Monte Carlo EM algorithm (Wei & Tanner, 1990;
Dempster et al., 1977). However, access to the oracle
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$ is unrealistic
and the exact sampling of
${p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$, as
required in EM, is generally intractable. To address this, we resort to (i)
approximate but computationally cheap conditional sampling methods for VAEs to
generate imputations (Rezende et al., 2014; Mattei & Frellsen, 2018a; Simkus &
Gutmann, 2023) and (ii) learning objectives for the model ${p_{\bm{\theta}}}$
and the variational distribution ${q_{\bm{\phi}}}$ that mitigate the pitfalls
caused by the missing data. We call the proposed approach DeMissVAE
(decomposed approach for handling missing data in VAEs).

We construct the variational distribution
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ for an
incomplete data-point ${\bm{x}}_{\mathrm{obs}}$ using a completed-data
variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
and an (approximate) imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$:

$\displaystyle{q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})=\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\left[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})\right].$
(9)

Assuming that the completed-data variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
well-represents the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$,
and that the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ draws plausible
imputations of the missing variables, then
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ will reasonably
represent ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ (see the
two right-most columns of fig. 1). In contrast to section 4.1 we here use a
continuous-mixture variational distribution, which is more flexible than a
finite-mixture distribution, albeit at an extra computational cost due to
sampling the (approximate) imputations (see appendix D).

We now derive the DeMissVAE objectives for fitting the generative model
${p_{\bm{\theta}}}({\bm{x}})$ and the completed-data variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$,
see appendix C for a more in-depth treatment.

##### Objective for ${p_{\bm{\theta}}}({\bm{x}},{\bm{z}})$.

With the variational distribution in eq. 9, we derive an ELBO on the marginal
log-likelihood, similar to eq. 2, to learn the parameters ${\bm{\theta}}$ of
the generative model:

$\displaystyle\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle\geq\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{z}})}{\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\left[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})\right]}\right]$
$\displaystyle=\underbrace{\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{z}})\right]}_{\overset{!}{=}\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}({\bm{x}}_{\mathrm{obs}};{\bm{\phi}},{\bm{\theta}},f^{t})}+\underbrace{\mathcal{H}\left[{q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\right]}_{\text{Const.
w.r.t.\ ${\bm{\theta}}$}}.$ (10)

This lower-bound can be further decomposed into log-likelihood and KL
divergence terms

$\displaystyle\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}({\bm{x}}_{\mathrm{obs}};{\bm{\phi}},{\bm{\theta}},f^{t})+\mathcal{H}\left[{q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\right]$
$\displaystyle=\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})-D_{\text{KL}}({q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\;||\;{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})),$
(11)

which means that if
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\approx{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
then maximising eq. 10 w.r.t. ${\bm{\theta}}$ performs approximate maximum-
likelihood estimation. Importantly, the missing variables
${\bm{x}}_{\mathrm{mis}}$ are marginalised-out, which adds robustness to the
potential sampling errors in
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$.

##### Objective for ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$.

We obtain the objective for learning the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ by marginalising the missing variables
${\bm{x}}_{\mathrm{mis}}$ from the complete-data ELBO in eq. 2 and then lower-
bounding the integral using
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ (see appendix
B):

$\displaystyle\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle\geq\underbrace{\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\right]}_{\overset{!}{=}\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}({\bm{x}}_{\mathrm{obs}};{\bm{\phi}},{\bm{\theta}},f^{t})}+\underbrace{\mathcal{H}\left[{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\right]}_{\text{Const.
w.r.t.\ ${\bm{\phi}}$}}.$ (12)

This lower-bound can also be decomposed into the log-likelihood term and two
KL divergence terms

$\displaystyle\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}({\bm{x}}_{\mathrm{obs}};{\bm{\phi}},{\bm{\theta}},f^{t})+\mathcal{H}\left[{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\right]$
$\displaystyle=\begin{aligned}
&\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})-D_{\text{KL}}({f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\;||\;{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}))\\\
&-\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\left[D_{\text{KL}}({q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})\;||\;{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}}))\right],\end{aligned}$
(13)

which means that the bound is maximised w.r.t. ${\bm{\phi}}$ iff
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
for all ${\bm{x}}_{\mathrm{mis}}$. Therefore, using the above objective to fit
${q_{\bm{\phi}}}$ corresponds directly to the complete-data case, and hence
the issues due to missingness identified in section 3 are mitigated.

If
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
for all ${\bm{x}}_{\mathrm{mis}}$, then maximising either of the bounds in
eqs. 12 and 10 w.r.t. the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ would correspond
to setting
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$.
However, directly learning an imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
is challenging (Simkus et al., 2023, Section 2.2). This motivates using
sampling methods to approximate the optimal imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
with samples. We draw samples from
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ using (cheap)
approximate conditional sampling methods for VAEs to obtain $K$ imputations
$\\{{\bm{x}}_{\mathrm{mis}}^{k}\\}_{k}^{K}$ and then use them to approximate
the expectations w.r.t.
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ in the above
objectives. We discuss the implementation of the algorithm in detail in
appendix D.

Finally, it is worth noting that the
$\mathcal{L}_{\mathrm{CVI}}^{{\bm{\theta}}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ objectives in eqs. 10 and 12 are
based on the standard ELBO. Extensions to the importance-weighted ELBO might
improve the method further by increasing the flexibility of the variational
posterior. However, unlike the standard ELBO used in eq. 10 where the density
of the imputation-based variational-mixture
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ can be dropped,
IWELBO requires computing the density of the proposal distribution
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, which is
generally intractable. We hence leave this direction for future work.

##  5 Related work

##### Fitting VAEs from incomplete data.

Since the seminal works of Kingma & Welling (2013) and Rezende et al. (2014),
VAEs have been widely used for density estimation from incomplete data and
various downstream tasks, primarily due to the computationally efficient
marginalisation of the model in eq. 1. Vedantam et al. (2017) and Wu & Goodman
(2018) explored the use of product-of-experts variational distributions,
drawing inspiration from findings in the factor analysis case with incomplete
data (Williams et al., 2018). Mattei & Frellsen (2019) used the importance-
weighted ELBO (Burda et al., 2015) for training VAEs on incomplete training
data sets. Ma et al. (2019) proposed the use of permutation invariant neural
networks to parametrise the encoder network instead of relying on zero-
masking. Nazábal et al. (2020) introduced hierarchical priors to handle
incomplete heterogeneous training data. Simkus et al. (2023) proposed a
general-purpose approach that is applicable to VAEs, not requiring the decoder
distribution to be easily marginalisable. Here, we further develop the
understanding of VAEs in the presence missing values in the training data set,
and propose variational-mixtures as a natural approach to improve VAE
estimation from incomplete data, building upon the motivation from imputation-
mixtures discussed in section 3.

##### Variational mixture distributions.

Mixture distributions have found widespread application in variational
inference and VAE literature. Roeder et al. (2017) introduced the stratified
ELBO corresponding to eq. 5. In the context of VAEs in multimodal domains, Shi
et al. (2019, Appendix A) introduced the stratified IWELBO corresponding to
eq. 8, but opted to use a looser bound instead, as detailed footnote 6. These
bounds were subsequently rediscovered by Morningstar et al. (2021) and Kviman
et al. (2023), who investigated their use for VAE estimation in fully-observed
data scenarios. Furthermore, Figurnov et al. (2019) introduced implicit
reparametrisation, enabling gradient estimation for ancestrally-sampled
mixtures, allowing the estimation of variational mixtures using eqs. 4 and 7.
Here, we build on the prior work, asserting that variational-mixtures are
well-suited for handling the posterior complexity increase due to missing data
(see section 3). Moreover, the imputation-mixture distribution used in
DeMissVAE is a novel type of variational mixtures specifically designed for
incomplete data scenarios.

##### Posterior complexity increase due to missing data.

Concurrent to this study, Sudak & Tschiatschek (2023) have recently brought
attention to a phenomenon related to the increase in posterior complexity due
to incomplete data, as discussed in section 3. They noted that, for any
${\bm{x}}_{\mathrm{obs}}$ and ${\bm{x}}_{\mathrm{obs}\setminus u}$, where $u$
is a subset of the observed dimensions, the model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}\setminus u})$ should
exhibit a strong dependency. However, because of the approximations in the
variational posterior (see e.g. Cremer et al., 2018; Zhang et al., 2021), the
variational approximations
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}\setminus u})$ may not
consistently capture this dependency. They refer to the lack of dependency
between ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}\setminus u})$, compared to
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}\setminus u})$, as
posterior inconsistency. Focused on improving downstream task performance,
they introduce regularisation into the VAE training objective to address
posterior inconsistency. In contrast to their work, we compare the fully-
observed and incomplete-data posteriors,
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, respectively. And,
with the goal of improving model estimation performance, we propose the use of
variational-mixtures to mitigate the posterior complexity gap between the
fully-observed and incomplete-data posteriors.

##### Marginalised variational bound.

In the standard ELBO derivation for incomplete data in eq. 2 the missing
variables are first marginalised (collapsed) from the likelihood, and then a
variational ELBO is established. This approach is sometimes referred to as
collapsed variational inference (CVI). In contrast, in the derivation of the
DeMissVAE encoder objective in eq. 12 we swap the order of marginalisation and
variational inference. Specifically, we start with the variational ELBO on
completed-data, and then marginalise the missing variables (see appendix B).
This approach bears similarity to the marginalised variational bound (MVB, or
KL-corrected bound) in exponential-conjugate variational inference literature
(King & Lawrence, 2006; Lázaro-Gredilla & Titsias, 2011; Hensman et al.,
2012). In these works, MVB has been preferred over CVI due to improved
convergence and guarantees that for appropriately formulated conjugate models
MVB is analytically tractable in cases where CVI is not (Hensman et al., 2012,
Section 3.3). While MVB remains intractable in the VAE setting with incomplete
data, similar to how the standard ELBO is intractable in fully-observed case,
we find the motivation behind MVB and DeMissVAE to be similar.

##  6 Evaluation

We here evaluate the proposed methods, MissVAE, MissSVAE, MissIWAE, MissSIWAE
(section 4.1), and DeMissVAE (section 4.2), on synthetic and real-world data,
and compare them to the popular methods MVAE and MIWAE that do not use mixture
variational distributions (Mattei & Frellsen, 2019). The methods are
summarised in table 1.

Method |  ${p_{\bm{\theta}}}$ objective |  ${q_{\bm{\phi}}}$ objective | # of components | Mixture sampling
---|---|---|---|---
MVAE$\dagger$ | eq. 4 | eq. 4 | $K=1$ | —
MissVAE | eq. 4 | eq. 4 | $K>1$ | Ancestral
MissSVAE | eq. 5 | eq. 5 | $K>1$ | Stratified
MIWAE$\dagger$ | eq. 7 | eq. 7 | $K=1$ | —
MissIWAE | eq. 7 | eq. 7 | $K>1$ | Ancestral
MissSIWAE | eq. 8 | eq. 8 | $K>1$ | Stratified
DeMissVAE | eq. 10 | eq. 12 | $K>1$ | Conditional VAE
Table 1: _Summary of the proposed and baseline methods._ The non-mixture
baselines ($\dagger$) are based on Mattei & Frellsen (2019) and the other
methods are proposed in this paper. Moreover, the methods using ancestral
sampling require implicit reparametrisation (Figurnov et al., 2019), whereas
the other methods work with the standard reparametrisation trick.

###  6.1 Mixture-of-Gaussians data with a 2D latent VAE

Evaluating log-likelihood on held-out data is generally intractable for VAEs
due to the intractable integral in eq. 1. We hence here choose a VAE with 2D
latent space, where numerical integration can be used to estimate the log-
likelihood of the model accurately (see section E.1 for more details). We fit
the model on incomplete data drawn from a mixture-of-Gaussians distribution.
By introducing missingness in the mixture-of-Gaussians data we introduce
multi-modality in the latent space (see fig. 1), which allows us to verify the
efficacy of mixture-variational distributions when the posteriors are multi-
modal due to missing data.

Figure 2: _Log-likelihood on held out data evaluated by numerically
integrating the 2D latent variables._ VAEs were fitted on mixture-of-Gaussians
data with 50% missingness. Each model is fitted with a computational budget of
5/15/25 samples from the variational distribution. The box plots show 1st and
3rd quartiles, the black lines are the medians, the dashed lines are the
means, and the whiskers show the data range over 5 independent runs. MVAE and
MIWAE ($\dagger$) are baseline methods by Mattei & Frellsen (2019). The other
five methods are proposed in this paper.

Results are shown in fig. 2. We first note that the stratified MissSVAE
approach performed better than MissVAE that uses ancestral sampling. The
reason for this is likely that stratified sampling reduces Monte Carlo
variance of the gradients w.r.t. ${\bm{\phi}}$ and hence enables a better fit
of the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ (see a further
investigation in section F.1.1). In line with this intuition, the MissVAE
results exhibit significantly larger variance than MissSVAE. Similarly, we
observe that the stratified MissSIWAE approach performed better than MissIWAE.
Importantly, we see that the use of mixture variational distributions in
MissSVAE and MissSIWAE improve the model fit over the MVAE and MIWAE baselines
that do not use mixtures to deal with the increased posterior complexity due
to missingness. Finally, we observe that DeMissVAE is capable of achieving
comparable performance to MIWAE and MissSIWAE, despite using a looser ELBO
bound, which shows that the decomposed approach to handling data missingness
can be used to achieve an improved fit of the model.

In section F.1.2, we analyse the model and variational posteriors of the
learnt models. We observe that the mixture approaches better-approximate the
incomplete-data posteriors, compared to the approaches that do not use
variational-mixtures. Moreover, we also observe that the structure of the
latent space is better-behaved when fitted using the decomposed approach in
DeMissVAE.

###  6.2 Real-world UCI data sets

Figure 3: _Estimate of the test log-likelihood using the IWELBO with
$I=50000$, on four UCI data sets._ Each data set was rendered incomplete by
applying uniform missingness of 20/50/80%. The curves show average performance
over 5 independent runs of the algorithms and the intervals show the 90%
confidence.

We here evaluate the proposed methods on real-world data sets from the UCI
repository (Dua & Graff, 2017; Papamakarios et al., 2017). We train a VAE
model with ResNet architecture on incomplete data sets with 20/50/80% uniform
missingness (see section E.2 for more details). We then estimate the log-
likelihood on complete test data set using the IWELBO bound with
$I=50\text{K}$ importance samples.777As $I\rightarrow\infty$ IWELBO approaches
$\log{p_{\bm{\theta}}}({\bm{x}})$. Moreover, as suggested by Mattei & Frellsen
(2018b), to improve the estimate on held-out data we fine-tune the encoder on
complete test data before estimating the log-likelihood. For additional
metrics see section F.2.

The results are shown in fig. 3. We first note that, similar to before, the
stratified MissSVAE approach performed better than MissVAE which uses
ancestral sampling. Importantly, we observe that using mixture variational
distributions in MissSVAE improves the fit of the model over MVAE (with the
exception on the Miniboone data set) that uses non-mixture variational
distributions. Furthermore, the gains in model accuracy typically increase
with data missingness, which verifies that MissSVAE performs better because it
handles the increased posterior complexity due to missing data better (see
fig. 1). Next, we observe that the performance of MIWAE, MissIWAE, and
MissSIWAE is similar, although we can note a small improvement by using
MissIWAE and MissSIWAE in large missingness settings. We observe only a
relatively small difference between the IWAE methods because the use of
importance weighted bound already corresponds to using a more flexible semi-
implicitly defined variational distribution (Cremer et al., 2017), which here
seems to be sufficient to deal with the complexities arising due to
missingness. Finally, we note that DeMissVAE results are in-between MissSVAE
and MIWAE. This verifies that the decomposed approach can be used to deal with
data missingness and, as a result, can improve the fit of the model.
Nonetheless, DeMissVAE is surpassed by the IWAE methods, which is likely due
to using the ELBO in DeMissVAE versus IWELBO in IWAE methods that can tighten
the bound more effectively.

###  6.3 MNIST and Omniglot data sets

Figure 4: _Estimate of the test log-likelihood using the IWELBO with $I=1000$,
MNIST and Omniglot data sets._ Each image in the data set was missing 2 out of
4 random quadrants. The box plots show 1st and 3rd quartiles, the black lines
are the medians, the dashed lines are the means, and the whiskers show the
data range over 5 independent runs.

In this section we evaluate the proposed methods on binarised MNIST (Garris et
al., 1991) and Omniglot (Lake et al., 2015) data sets of handwritten
characters. We fit a VAE model with a convolutional ResNet encoder and decoder
networks (see section E.3 for more details). The data is made incomplete by
masking 2 out of 4 quadrants of an image at random. Similar to the previous
section, we estimate the log-likelihood on a complete test data set using the
IWELBO bound with $I=1000$ importance samples.

On the MNIST data set we see that MVAE $\leq$ MissVAE $<$ MissSVAE similar to
the previous results but MIWAE $<$ MissSIWAE $<$ MissIWAE. This suggests that
MissIWAE, which uses ancestral sampling, was able to tighten the bound more
effectively compared to stratified MissSIWAE, and was able to fit the
variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ well despite the
potentially larger variance w.r.t. ${\bm{\phi}}$. Moreover, we also see that
MVAE $<$ MIWAE $<$ DeMissVAE, which further verifies that the decomposed
approach is able to handle the data missingness well.

On the Omniglot data we observe that the mixture approaches perform similarly
to MVAE and MIWAE, which do not use mixture variational distributions. This
suggests that either the posterior multi-modality is less prominent in the
Omniglot data set or that due to the reverse KL optimisation of the
variational distribution all mixture components have degenerated to a single
mode. Finally, DeMissVAE slightly outperforms MVAE, MissVAE, and MissSVAE, but
is surpassed by the importance-weighted approaches.

Interestingly, in this evaluation the stratified approaches (MissSVAE and
MissSIWAE) were outperformed by the approaches using standard ELBO and
implicit reparametrisation (MissVAE and MissIWAE). This suggests that the
performance of each approach can be data- and model-dependent and hence both
should be evaluated when possible.

##  7 Discussion

Handling missing data is a crucial task in machine learning for the
application of modern statistical methods in many practical domains
characterised by incomplete data. In the context of variational autoencoders
we have shown that incomplete data introduces posterior complexity over the
latent variables, compared to the fully-observed data case. Consequently,
accurate model fitting from incomplete data requires the use of more flexible
variational families compared to the complete case. We have then stipulated
that variational-mixtures are a natural approach for handling data
missingness. This allows us to work with the same variational families that
are known to work well when the data is fully-observed, enabling the transfer
of useful known inductive biases (Miao et al., 2022) from the fully-observed
to the incomplete data scenario.

Subsequently, we have introduced two methodologies grounded in variational
mixtures. First, we proposed using finite variational mixtures with the
standard and importance-weighted ELBOs using ancestral and stratified sampling
of the mixtures. Additionally, we have proposed a novel “decomposed”
variational-mixture approach, that uses cost-effective yet often coarse
conditional sampling methods for VAEs to generate imputations and ELBO-based
objectives that are robust to the sampling errors.

Our evaluation shows that using variational mixtures can improve the fit of
VAEs when dealing with incomplete data, surpassing the performance of models
without variational mixtures. Moreover, our observations indicate that,
although stratified sampling of the finite mixtures often yields better
results compared to ancestral sampling, the effectiveness of these methods can
be data- and model-dependent and hence both approaches should be evaluated
when possible. In our findings, we further note that the decomposed approach
in DeMissVAE outperforms all ELBO-based methods but falls short of surpassing
IWELBO-based methods. These results point towards promising research avenues,
suggesting potential improvements in VAE model estimation from incomplete
data. Future directions include extending the DeMissVAE approach to
incorporate IWELBO-based objectives and developing improved cost-effective
conditional sampling methods for VAEs.

| Budget |  Variational families |  Latent structure* | Evaluation rank
---|---|---|---|---
Method |  | _MoG_ | _UCI_ |  _MNIST_ +_Omniglot_
MissVAE | _small_ | _limited_ | _well-behaved_ | 5 | 5 | 5
MissSVAE | _medium_ | _any_ | _well-behaved_ | 4 | 4 | 4
MissIWAE | _medium_ | _limited_ | _potentially irregular_ | 3 | 2 | 1
MissSIWAE | _medium_ | _any_ | _potentially irregular_ | 1 | 1 | 2
DeMissVAE | _medium/high_ | _any_ |  _well-behaved_ + | 2 | 3 | 3
Table 2: _A coarse summary of advantages and disadvantages of the proposed
methods._ _Budget_ : small/medium/high depending on the number of latent
samples required or whether conditional sampling of VAEs is needed.
_Variational families_ : which families of distributions can be used as
mixture components—any reparametrisable families, or limited families, as
discussed in section 4.1. _Latent structure_ : methods with potential to learn
irregular latent spaces may have decreased downstream performance in certain
tasks. We have found ($+$) that DeMissVAE is able to achieve the most well-
behaved latent structures on the MoG data in section F.1.2. Please note (*)
that the learnt latent structure will depend on the chosen model architecture.
_Evaluation rank_ : the rank of the proposed methods in the evaluations in
sections 6.1, 6.2 and 6.3.

The choice between the proposed methods for fitting VAEs from incomplete data
depends on various factors such as computational budget, variational families,
model accuracy goals, and the specific requirements of downstream tasks,
discussed as follows and summarised in table 2.

Computational and memory budget.

The standard ELBO with ancestral sampling is the most suitable method for
small computational and memory budgets, since the objective can be estimated
using a single latent sample for each data point. On the other hand, methods
using stratified sampling or the importance-weighted ELBO require the use of
multiple latent samples for each data-point and hence may only be used if the
memory and compute budget allows. Furthermore, akin to the standard ELBO, the
DeMissVAE objectives can be estimated using a single latent sample, but the
approach incurs extra cost in sampling the imputations.

Variational families.

While the stratified and DeMissVAE approaches can use any reparametrisable
distribution family for the mixture components, the ancestral sampling methods
require the use of _implicit_ reparametrisation (Figurnov et al., 2019) and as
a result may not work with all distribution families (see discussion in
section 4.1).

Model accuracy.

Stratified sampling of mixtures can improve the model accuracy, compared to
ancestral sampling, by reducing Monte Carlo gradient variance. Additionally,
methods using the importance-weighted ELBO, compared to the standard ELBO, are
often able to tighten the bound more effectively by using multiple importance
samples, leading to improved model accuracy. DeMissVAE performance lies in
between the standard ELBO and importance-weighted ELBO approaches. Although
the introduced DeMissVAE objectives exhibit robustness to some imputation
distribution error, improved model accuracy can often be achieved by improving
the accuracy of imputations by using a larger budget for the imputation step.

Latent structure.

Different downstream tasks may prefer distinct latent structures, for example,
conditional generation from unconditional VAEs is often easier if the latent
space is well-structured (Engel et al., 2017; Gómez-Bombarelli et al., 2018).
To this end, observations in section F.1.2 show that the latent space of
DeMissVAE behaves well, and is comparable to a model fitted with complete
data. This characteristic makes it preferable for downstream tasks requiring
well-structured latent spaces. On the other hand, as noted by Burda et al.
(2015, Appendix C) and Cremer et al. (2018, Section 5.4), the use of
importance-weighted ELBO to mitigate the increased posterior complexity due to
missing data, may make the latent space less regular, compared to a model
trained on fully-observed data set, which potentially decreases the model’s
performance on downstream tasks.

Finally, we step back to note that this paper is focused on the class of
variational autoencoder models, a subset of the broader family of deep latent
variable models (DLVMs). Much like VAEs, DLVMs usually aim to efficiently
represent the intricate nature of data through a well-structured latent space,
implicitly defined by a learnable generative process. Building on our findings
in VAEs, where incomplete data led to an increased complexity in the posterior
distribution compared to the fully-observed case, we conjecture that a similar
effect may occur within the wider family of DLVMs, affecting the fit of the
model. We therefore believe that there is substantial scope to explore the
implications of incomplete data in other DLVM classes, particularly focusing
on the effects of marginalisation on latent space representations and the
associated generative processes. Investigating decomposed approaches, similar
to DeMissVAE or Monte Carlo EM (Wei & Tanner, 1990), presents promising
avenues for further research in this direction.

## References

  * Burda et al. (2015) Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.  Importance Weighted Autoencoders.  In _International Conference on Learning Representations (ICLR)_ , San Juan, Puerto Rico, September 2015.
  * Cremer et al. (2017) Chris Cremer, Quaid Morris, and David Duvenaud.  Reinterpreting Importance-Weighted Autoencoders.  In _ICLR Workshop_ , February 2017.
  * Cremer et al. (2018) Chris Cremer, Xuechen Li, and David Duvenaud.  Inference Suboptimality in Variational Autoencoders.  In _International Conference on Machine Learning (ICML)_ , May 2018.
  * Dempster et al. (1977) Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin.  Maximum Likelihood from Incomplete Data Via the EM Algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.  doi: 10.1111/j.2517-6161.1977.tb01600.x.
  * Dinh et al. (2017) Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using Real NVP.  In _International Conference on Learning Representations (ICLR)_ , February 2017.
  * Dua & Graff (2017) Dheeru Dua and Casey Graff.  UCI Machine Learning Repository, 2017.
  * Engel et al. (2017) Jesse Engel, Matthew Hoffman, and Adam Roberts.  Latent Constraints: Learning to Generate Conditionally from Unconditional Generative Models.  In _International Conference on Learning Representations (ICLR)_ , December 2017.
  * Figurnov et al. (2019) Michael Figurnov, Shakir Mohamed, and Andriy Mnih.  Implicit Reparameterization Gradients.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , January 2019.
  * Garris et al. (1991) Michael D. Garris, R. A. Wilkinson, and Charles L. Wilson.  Methods for enhancing neural network handwritten character recognition.  In _International Joint Conference on Neural Networks (IJCNN)_ , volume 1, pp. 695–700, Seattle, WA, USA, 1991. IEEE.  ISBN 978-0-7803-0164-1.  doi: 10.1109/IJCNN.1991.155265.
  * Gershman & Goodman (2014) Samuel J. Gershman and Noah D. Goodman.  Amortized Inference in Probabilistic Reasoning.  In _Annual Meeting of the Cognitive Science Society_ , volume 36, 2014.
  * Gómez-Bombarelli et al. (2018) Rafael Gómez-Bombarelli, Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Sánchez-Lengeling, Dennis Sheberla, Jorge Aguilera-Iparraguirre, Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik.  Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules.  _ACS Central Science_ , 4(2):268–276, February 2018.  ISSN 2374-7943.  doi: 10.1021/acscentsci.7b00572.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative Adversarial Networks.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , June 2014.
  * Hensman et al. (2012) James Hensman, Magnus Rattray, and Neil D. Lawrence.  Fast Variational Inference in the Conjugate Exponential Family.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , December 2012.
  * Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2017.
  * Ipsen et al. (2020) Niels Bruun Ipsen, Pierre-Alexandre Mattei, and Jes Frellsen.  Not-MIWAE: Deep Generative Modelling with Missing not at Random Data.  In _International Conference on Learning Representations (ICLR)_ , June 2020.
  * King & Lawrence (2006) Nathaniel J. King and Neil D. Lawrence.  Fast Variational Inference for Gaussian Process Models Through KL-Correction.  In _European Conference on Machine Learning (ECML)_ , 2006.
  * Kingma & Welling (2013) Diederik P. Kingma and Max Welling.  Auto-Encoding Variational Bayes.  In _International Conference on Learning Representations (ICLR)_ , December 2013.
  * Krishnan et al. (2016) Rahul G. Krishnan, Uri Shalit, and David Sontag.  Structured Inference Networks for Nonlinear State Space Models.  In _AAAI Conference on Artificial Intelligence_ , December 2016.  doi: 10.48550/arXiv.1609.09869.
  * Kviman et al. (2023) Oskar Kviman, Ricky Molén, Alexandra Hotti, Semih Kurt, Víctor Elvira, and Jens Lagergren.  Cooperation in the Latent Space: The Benefits of Adding Mixture Components in Variational Autoencoders.  In _International Conference on Machine Learning (ICML)_ , July 2023.  doi: 10.48550/arXiv.2209.15514.
  * Lake et al. (2015) Brenden M. Lake, Ruslan Salakhutdinov, and Joshua B. Tenenbaum.  Human-level concept learning through probabilistic program induction.  _Science_ , 2015.  doi: 10.1126/science.aab3050.
  * Lázaro-Gredilla & Titsias (2011) Miguel Lázaro-Gredilla and Michalis K. Titsias.  Variational heteroscedastic Gaussian process regression.  In _International Conference on Machine Learning (ICML)_ , June 2011.
  * Little & Rubin (2002) Roderick J. A. Little and Donald B. Rubin.  _Statistical Analysis with Missing Data: Second Edition_.  Wiley-Interscience, 2002.  ISBN 0-471-18386-5.
  * Ma & Zhang (2021) Chao Ma and Cheng Zhang.  Identifiable Generative Models for Missing Not at Random Data Imputation.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , October 2021.
  * Ma et al. (2019) Chao Ma, Sebastian Tschiatschek, Konstantina Palla, José Miguel Hernández-Lobato, Sebastian Nowozin, and Cheng Zhang.  EDDI: Efficient dynamic discovery of high-value information with partial VAE.  In _International Conference on Machine Learning (ICML)_ , pp. 7483–7504, 2019.  ISBN 9781510886988.
  * Mattei & Frellsen (2018a) Pierre-Alexandre Mattei and Jes Frellsen.  Leveraging the Exact Likelihood of Deep Latent Variable Models.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , February 2018a.
  * Mattei & Frellsen (2018b) Pierre-Alexandre Mattei and Jes Frellsen.  Refit your Encoder when New Data Comes by.  In _Workshop on Bayesian Deep Learning at Neural Information Processing Systems (NeurIPS)_ , pp. 4, Montreal, Canada, 2018b.
  * Mattei & Frellsen (2019) Pierre-Alexandre Mattei and Jes Frellsen.  MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets.  In _International Conference on Machine Learning (ICML)_ , 2019.
  * Meng (1994) Xiao-Li Meng.  On the Rate of Convergence of the ECM Algorithm.  _The Annals of Statistics_ , 22(1):326–339, March 1994.  ISSN 0090-5364, 2168-8966.  doi: 10.1214/aos/1176325371.
  * Miao et al. (2022) Ning Miao, Emile Mathieu, N. Siddharth, Yee Whye Teh, and Tom Rainforth.  On Incorporating Inductive Biases into VAEs.  In _International Conference on Learning Representations (ICLR)_ , February 2022.  doi: 10.48550/arXiv.2106.13746.
  * Morningstar et al. (2021) Warren Morningstar, Sharad Vikram, Cusuh Ham, Andrew Gallagher, and Joshua Dillon.  Automatic Differentiation Variational Inference with Mixtures.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pp. 3250–3258. PMLR, March 2021.
  * Nazábal et al. (2020) Alfredo Nazábal, Pablo M. Olmos, Zoubin Ghahramani, and Isabel Valera.  Handling Incomplete Heterogeneous Data using VAEs.  _Pattern Recognition_ , 107, 2020.  ISSN 0031-3203.  doi: 10.1016/j.patcog.2020.107501.
  * Papamakarios et al. (2017) George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked Autoregressive Flow for Density Estimation.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 30, 2017.
  * Papamakarios et al. (2021) George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing Flows for Probabilistic Modeling and Inference.  _Journal of Machine Learning Research_ , 22(57):1–64, 2021.
  * Rainforth et al. (2019) Tom Rainforth, Adam R. Kosiorek, Tuan Anh Le, Chris J. Maddison, Maximilian Igl, Frank Wood, and Yee Whye Teh.  Tighter Variational Bounds are Not Necessarily Better.  In _International Conference on Machine Learning (ICML)_ , March 2019.
  * Reddi et al. (2018) Sashank J. Reddi, Satyen Kale, and Sanjiv Kumar.  On the convergence of Adam and beyond.  In _International Conference on Learning Representations (ICLR)_ , pp. 1–23, 2018.
  * Rezende et al. (2014) Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic Backpropagation and Approximate Inference.  In _International Conference on Machine Learning (ICML)_ , Beijing, China, 2014.
  * Robert & Casella (2004) Christian P. Robert and George Casella.  _Monte Carlo Statistical Methods_.  Springer, 2004.  ISBN 0-387-21239-6.
  * Roeder et al. (2017) Geoffrey Roeder, Yuhuai Wu, and David K. Duvenaud.  Sticking the Landing: Simple, Lower-Variance Gradient Estimators for Variational Inference.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 30, 2017.
  * Shi et al. (2019) Yuge Shi, N. Siddharth, Brooks Paige, and Philip Torr.  Variational Mixture-of-Experts Autoencoders for Multi-Modal Deep Generative Models.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Simkus & Gutmann (2023) Vaidotas Simkus and Michael U. Gutmann.  Conditional Sampling of Variational Autoencoders via Iterated Approximate Ancestral Sampling.  _Transactions on Machine Learning Research_ , August 2023.  ISSN 2835-8856.
  * Simkus et al. (2023) Vaidotas Simkus, Benjamin Rhodes, and Michael U. Gutmann.  Variational Gibbs Inference for Statistical Model Estimation from Incomplete Data.  _Journal of Machine Learning Research_ , 24(196):1–72, 2023.  ISSN 1533-7928.
  * Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep Unsupervised Learning using Nonequilibrium Thermodynamics.  In _International Conference on Machine Learning (ICML)_ , November 2015.  doi: 10.48550/arXiv.1503.03585.
  * Sudak & Tschiatschek (2023) Timur Sudak and Sebastian Tschiatschek.  Posterior Consistency for Missing Data in Variational Autoencoders.  In _European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases (ECML-PKDD)_ , October 2023.  doi: 10.48550/arXiv.2310.16648.
  * Tieleman (2008) Tijmen Tieleman.  Training restricted Boltzmann machines using approximations to the likelihood gradient.  In _International Conference on Machine Learning (ICML)_ , pp. 1064–1071, 2008.  ISBN 9781605582054.  doi: 10.1145/1390156.1390290.
  * Tucker et al. (2018) George Tucker, Dieterich Lawson, Shixiang Gu, and Chris J. Maddison.  Doubly Reparameterized Gradient Estimators for Monte Carlo Objectives.  In _International Conference on Learning Representations (ICLR)_ , November 2018.
  * Vedantam et al. (2017) Ramakrishna Vedantam, Ian Fischer, Jonathan Huang, and Kevin P. Murphy.  Generative Models of Visually Grounded Imagination.  _International Conference on Learning Representations (ICLR)_ , May 2017.
  * Wei & Tanner (1990) Greg C. G. Wei and Martin A. Tanner.  A Monte Carlo Implementation of the EM Algorithm and the Poor Man’s Data Augmentation Algorithms.  _Journal of the American Statistical Association_ , 85(411):699–704, September 1990.  doi: 10.1080/01621459.1990.10474930.
  * Williams et al. (2018) Christopher K. I. Williams, Charlie Nash, and Alfredo Nazábal.  Autoencoders and Probabilistic Inference with Missing Data: An Exact Solution for The Factor Analysis Case.  _arXiv preprint_ , 1801.03851, January 2018.
  * Wu & Goodman (2018) Mike Wu and Noah D. Goodman.  Multimodal Generative Models for Scalable Weakly-Supervised Learning.  In _NeurIPS 2018_ , February 2018.
  * Younes (1999) Laurent Younes.  On the convergence of markovian stochastic algorithms with rapidly decreasing ergodicity rates.  _Stochastics and Stochastic Reports_ , 65(3-4):177–228, February 1999.  ISSN 1045-1129.  doi: 10.1080/17442509908834179.
  * Zhang et al. (2021) Mingtian Zhang, Peter Hayes, and David Barber.  Generalization Gap in Amortized Inference.  In _Workshop on Bayesian Deep Learning at Neural Information Processing Systems (NeurIPS)_ , pp. 6, 2021.

##  Appendix A Posterior complexity due to missing information

The complexity increase of the model posterior due to missing data, shown in
fig. 1, explains why flexible variational distributions (Burda et al., 2015;
Cremer et al., 2017) have been preferred when fitting VAEs from incomplete
data (Mattei & Frellsen, 2019; Ipsen et al., 2020; Ma & Zhang, 2021). We here
define the increase of the posterior complexity via the expected
Kullback–Leibler (KL) divergence as follows

$\displaystyle\mathbb{E}_{{p^{*}}({\bm{x}})}\left[D_{\text{KL}}({p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})\;||\;{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}}))\right]$
$\displaystyle=\mathbb{E}_{{p^{*}}({\bm{x}})}{\mathbb{E}_{{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}{{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})}\right]}=\mathcal{I}({\bm{z}},{\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}).$

As shown above the expected KL divergence equals the (conditional) mutual
information (MI) between the latents ${\bm{z}}$ and the missing variables
${\bm{x}}_{\mathrm{mis}}$.

The mutual information interpretation allows us to reason when a more flexible
variational family may be necessary to accurately estimate VAEs from
incomplete data. Specifically, when the MI is small then the two posterior
distributions, ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ are similar, in which
case a simple variational distribution may work sufficiently well. This
situation might appear when the observed ${\bm{x}}_{\mathrm{obs}}$ and
unobserved ${\bm{x}}_{\mathrm{mis}}$ variables are highly related and
${\bm{x}}_{\mathrm{mis}}$ provides little additional information about
${\bm{z}}$ over just ${\bm{x}}_{\mathrm{obs}}$, for example, when random
pixels of an image are masked it is “easy” to infer the complete image due to
strong relationship between neighbouring pixels. On the other hand, when the
MI is high then ${\bm{x}}_{\mathrm{mis}}$ provides significant additional
information about ${\bm{z}}$ over just ${\bm{x}}_{\mathrm{obs}}$, in which
case a more flexible variational family may be needed, for example, when the
pixels of an image are masked in blocks such that it introduces significant
uncertainty about what is missing.

##  Appendix B DeMissVAE: Encoder objective derivation

We derive the objective for learning the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ in DeMissVAE by first marginalising
the missing variables ${\bm{x}}_{\mathrm{mis}}$ from the complete-data ELBO in
eq. 2:

$\displaystyle\log\int$
$\displaystyle\exp\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}\geq\log\int\exp\left\\{\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\right]\right\\}\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}$
The l.h.s. is the marginal log-likelihood
$\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$, but the integral in r.h.s.
is intractable. We then lower-bound the integral in the r.h.s. using the
imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ and Jensen’s
inequality $\displaystyle\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
$\displaystyle\geq\log\int\frac{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\exp\left\\{\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\right]\right\\}\mathop{}\\!\mathrm{d}{\bm{x}}_{\mathrm{mis}}$
$\displaystyle=\log\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\left[\exp\left(-\log{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\right)\exp\left\\{\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\right]\right\\}\right]$
$\displaystyle=\log\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\left[\exp\left\\{\mathbb{E}_{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}}){f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\right]\right\\}\right]$
$\displaystyle\geq\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}}){f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}\right]$
$\displaystyle=\underbrace{\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\left[\log\frac{{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}},{\bm{z}})}{{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})}\right]}_{\overset{!}{=}\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}({\bm{x}}_{\mathrm{obs}};{\bm{\phi}},{\bm{\theta}},f^{t})}+\underbrace{\mathcal{H}\left[{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\right]}_{\text{Const.
w.r.t.\ ${\bm{\phi}}$}}.$

##  Appendix C DeMissVAE: Motivating the separation of objectives

The two DeMissVAE objectives $\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eqs. 10 and 12 correspond to
valid lower-bounds on $\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
irrespective of ${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$.
Moreover, both of them are tight at
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
and
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$.
_So, a natural question is why do we prefer
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ to learn ${p_{\bm{\theta}}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ to learn ${q_{\bm{\phi}}}$?_

##### Why use $\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ in eq. 10 over
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eq. 12 to learn
${p_{\bm{\theta}}}({\bm{x}})$?

Maximisation of the objective $\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in
iteration $t$ w.r.t. ${\bm{\theta}}$ would have to compromise between
maximising the log-likelihood $\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
and keeping the other two KL divergence terms in eq. 13 low. Specifically, the
compromise between maximising $\log{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}})$
and keeping
$D_{\text{KL}}({f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\;||\;{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}))$
low is equivalent to the compromise in the EM algorithm, which is known to
affect the convergence of the model (Meng, 1994). Moreover, if
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\not={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
then minimising the
$D_{\text{KL}}({f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\;||\;{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}))$
will fit the model ${p_{\bm{\theta}}}({\bm{x}})$ to the biased samples from
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$. On the other
hand, in $\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ the missing variables
${\bm{x}}_{\mathrm{mis}}$ are marginalised from the model, therefore it avoids
the compromise with
$D_{\text{KL}}({f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\;||\;{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}))$
and the potential bias of the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ affects the
model _only_ via the latents
${\bm{z}}\sim{q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$,
increasing the robustness to sub-optimal imputations.

##### Why use $\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eq. 12 over
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ in eq. 10 to learn
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$?

In the case of $\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$, if
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
then the bound is tightened when
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
for all ${\bm{x}}_{\mathrm{mis}}$, which is the same optimal ${q_{\bm{\phi}}}$
if we used $\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$. But, there is also at
least one more possible optimal solution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$,
which ignores the imputations and corresponds to the optimal solution of the
standard approach in section 2, and thus it means that the optimum is
(partially) unidentifiable and can make optimisation of ${q_{\bm{\phi}}}$
using $\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ difficult. Moreover, if
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\not={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
then in order to minimise
$D_{\text{KL}}({q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})\;||\;{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}}))$
w.r.t. ${\bm{\phi}}$ the variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
would have to compensate for the inaccuracies of
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ by adjusting the
probability mass over the latents ${\bm{z}}$, such that
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
is correct on average, i.e.
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})=\mathbb{E}_{{f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}}[{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})]\approx{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$.
These two issues make optimising ${\bm{\phi}}$ via
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ such that
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})\approx{p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
difficult. On the other hand, in $\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ the
optimal ${q_{\bm{\phi}}}$ is always at
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$,
irrespective of the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$, hence the
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ objective in eq. 12 is well-defined
and more robust to inaccuracies of
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ for the
optimisation of
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$.

Figure 5: _A control study on a VAE model with 2D latent space (see additional
details in section E.1), examining the sensitivity of the proposed method
(DeMissVAE, green) and two control methods (blue and yellow) to the accuracy
of the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$._ Left:
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
represented using rejection sampling. Center: an oracle imputation function
that gets progressively “wider” from left-to-right of the figure. Right: an
oracle imputation distribution that towards the right of the figure more
significantly oversamples low-probability posterior modes. The log-likelihood
is computed on a held-out test data set by numerically integrating the 2D
latent space of the VAE. The horizontal axis on the two right-most figures
shows the Jensen–Shannon divergence between the imputation distribution and
the ground-truth conditional
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$.

In fig. 5 we verify the efficacy of DeMissVAE via a control study on a small
VAE model ${p_{\bm{\theta}}}({\bm{x}})$ with 2D latent space fitted on
incomplete samples from a ground truth mixture-of-Gaussians (MoG) distribution
${p^{*}}({\bm{x}})$. We evaluate fitting the VAE using only
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ in eq. 10 (CVI-VAE, blue), only
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eq. 12 (MVB-VAE, yellow), and
using the proposed two-objective approach (DeMissVAE, green). In the left-most
figure we evaluate the three methods where we represent the imputation
distribution
$f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$
using rejection sampling, which corresponds to the optimal imputation
distribution w.r.t.
$D_{\text{KL}}({f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\;||\;{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}}))=0$.
We see that the proposed approach (green) dominates over the other two control
methods (blue and yellow), and importantly that marginalisation of the missing
variables in DeMissVAE (green) improves the model accuracy compared to an EM-
type handling of the missing variables (yellow). Furthermore, in the remaining
two figures we investigate the sensitivity of the methods to the accuracy of
imputations in ${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$.
In Oracle 1 we start with the ground-truth conditional
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$ and, along the
x-axis of the figure, investigate how the methods perform when the imputation
distribution becomes “wider”: first interpolating from
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$ to an
independent unconditional distribution
$\prod_{d\in\text{idx}({\bm{m}})}{p^{*}}(x_{d})$ and then further towards and
independent Gaussian distribution. And in Oracle 2 we investigate what happens
when the sampler “oversamples” posterior modes: we interpolate the imputation
distribution from
${p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$ to
$\frac{1}{C}\sum_{c}^{C}{p^{*}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}},c)$,
where $c$ is the component of the mixture distribution with a total of $C$
components. As we see in the figure, the proposed DeMissVAE approach (green)
performs similar or better than the MVB-VAE (yellow) and CVI-VAE (blue)
control methods, with an exception when the
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ are extremely
inaccurate (last two points on the middle figure) which is expected since
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ in eq. 9 can be
arbitrarily far from ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
when
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})={p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}},{\bm{x}}_{\mathrm{mis}})$
but
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\not\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$.

Finally, in fig. 6 we investigate what happens if we used only
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ in eq. 10 or
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eq. 12 to fit the VAE model, in
contrast to the two separate objectives for encoder and decoder in DeMissVAE.
We use the LAIR sampling method (Simkus & Gutmann, 2023) as detailed in
appendix D to obtain approximate samples from
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$.
And, we observe that DeMissVAE achieves a better fit of the model, in line
with our motivation in this section.

Figure 6: _A control study on a VAE model with 2D latent space (see additional
details in section E.1), investigating the importance of the two-objective
approach in DeMissVAE (green) and two control methods (blue and yellow)._ In
CVI-VAE (blue) we fit both the encoder and decoder using eq. 10, and in MVB-
VAE (yellow) we fit both the encoder and decoder using eq. 12. The log-
likelihood is computed on a held-out test data set by numerically integrating
the 2D latent space of the VAE.

##  Appendix D DeMissVAE: Implementing the training procedure

DeMissVAE requires optimising two objectives
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eqs. 10 and 12 and drawing
(approximate) samples to represent
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}\approx{p_{\bm{\theta}}}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})$.
Our aim is to implement this efficiently to minimise redundant computation.

Algorithm 1 Shared computation of the DeMissVAE learning objectives

1:parameters ${\bm{\theta}}$ and ${\bm{\phi}}$, number of latent samples $L$,
completed data-point
$({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik})$

2:${\bm{\psi}}^{ik}\leftarrow\texttt{Encoder}({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik};{\bm{\phi}})$
$\triangleright$ Compute parameters of the variational distribution

3:${\bm{z}}_{1},\ldots,{\bm{z}}_{L}\sim
q({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik};{\bm{\psi}}^{ik})$
$\triangleright$ Sample latents ${\bm{z}}$

4:${\bm{\eta}}_{l}\leftarrow\texttt{Decoder}({\bm{z}}_{l};{\bm{\theta}})$ for
$\forall l\in[1,L]$ $\triangleright$ Compute parameters of the generative
distribution

5:def
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$(${\bm{z}}_{1},\ldots,{\bm{z}}_{L},{\bm{\eta}}_{1},\ldots,{\bm{\eta}}_{L}$):
$\triangleright$ Procedure for estimating eq. 10

6: return $\frac{1}{L}\sum_{l=1}^{L}\log
p({\bm{x}}_{\mathrm{obs}}^{i},{\bm{z}}_{l};{\bm{\eta}}_{l})$

7:$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}({\bm{z}}_{1},\ldots,{\bm{z}}_{L},{\bm{\eta}}_{1},\ldots,{\bm{\eta}}_{L}),\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}({\bm{\psi}}^{ik},{\bm{z}}_{1},\ldots,{\bm{z}}_{L},{\bm{\eta}}_{1},\ldots,{\bm{\eta}}_{L})$

8:def
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$(${\bm{\psi}}^{ik},{\bm{z}}_{1},\ldots,{\bm{z}}_{L},{\bm{\eta}}_{1},\ldots,{\bm{\eta}}_{L}$):
$\triangleright$ Procedure for estimating eq. 12

9: return $\frac{1}{L}\sum_{l=1}^{L}\log
p({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik},{\bm{z}}_{l};{\bm{\eta}}_{l})-\log
q({\bm{z}}_{l}\mid{\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik};{\bm{\psi}}^{ik})$

The algorithm starts with a randomly-initialised target VAE model
${p_{\bm{\theta}}}({\bm{x}},{\bm{z}})$, an amortised variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$, and an incomplete data set
$\mathcal{D}=\\{{\bm{x}}_{\mathrm{obs}}^{i}\\}_{i}$. And then, to represent
the imputation distribution
${f^{0}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$, $K$ imputations
$\\{{\bm{x}}_{\mathrm{mis}}^{ik}\\}_{k=1}^{K}$ are generated for each
${\bm{x}}_{\mathrm{obs}}^{i}\in\mathcal{D}$ using some simple imputation
function such as sampling the marginal empirical distributions of the missing
variables. The algorithm then iterates between the following two steps:

  1. 1.

Imputation. Update the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ using cheap
approximate sampling methods such as pseudo-Gibbs (Rezende et al., 2014),
Metropolis-within-Gibbs (MWG, Mattei & Frellsen, 2018a), or latent-adaptive
importance resampling (LAIR, Simkus & Gutmann, 2023). Moreover, since the
model and the variational distributions are initialised randomly, we skip the
imputation step during the first epoch over the data.

  2. 2.

Parameter update. Update the parameters using stochastic gradient ascent on
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eqs. 10 and 12 with the
imputations from
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$.

##### Efficient parameter update.

While the two objectives for ${p_{\bm{\theta}}}$ and ${q_{\bm{\phi}}}$ in eqs.
10 and 12 are different, a major part of the computation can be shared, as
shown in algorithm 1. As usual, the objectives are approximated using Monte
Carlo averaging and require only one evaluation of the generative model,
including the encoder, decoder, and prior, for each completed data-point
$({\bm{x}}_{\mathrm{obs}}^{i},{\bm{x}}_{\mathrm{mis}}^{ik})$. Therefore, only
backpropagation needs to be performed separately and the overall per-iteration
computational cost of optimising the two objectives is about 1.67 times the
cost of a fully-observed VAE optimisation (instead of 2 times if implemented
naïvely).999The cost of backpropagation is about 2 times the cost of a forward
pass (Burda et al., 2015).

##### Efficient imputation.

To make the imputation step efficient, the imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ is “persitent”
between iterations, that is, the imputation distribution from the previous
iteration ${f^{t-1}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ is
used to initialise the iterative approximate VAE sampler at iteration
$t$.101010Persistent samplers have been used in the past to increase
efficiency of maximum-likelihood estimation methods (Younes, 1999; Tieleman,
2008; Simkus et al., 2023). Moreover, an iteration of a pseudo-Gibbs, MWG, or
LAIR samplers uses the same quantities as the objectives
$\mathcal{L}_{\mathrm{CVI}}^{\bm{\theta}}$ and
$\mathcal{L}_{\mathrm{LMVB}}^{\bm{\phi}}$ in eqs. 10 and 12, and hence the
cost of one iteration of the sampler in the imputation step can be shared with
the cost of computation of the learning objectives. However, it is important
to note that the accuracy of imputations affects the accuracy of the estimated
model, and hence better estimation can be achieved by increasing the
computational budget for imputation.

##  Appendix E Experiment details

In this appendix we provide additional details on the experiments.

###  E.1 Mixture-of-Gaussians data with a 2D latent VAE

We generated a random 5D mixture-of-Gaussians model with 15 components by
sampling the mixture covariance matrices from the inverse Wishart distribution
$\mathcal{W}^{-1}(\nu=D,\Psi=\mathbf{I})$, means from the Gaussian
distribution $\mathcal{N}(\mu=\mathbf{0},\sigma=\mathbf{3})$ and the component
probabilities from Dirichlet distribution $\mathrm{Dir}(\alpha=\mathbf{1})$
(uniform). The model was then standardised to have a zero mean and a standard
deviation of one. The pairwise marginal densities of the generated
distribution is visualised in fig. 7 showing a highly-complex and multimodal
distribution, and the generated parameters and data used in this paper are
available in the shared code repository. We simulated a 20K sample data set
used to fit the VAEs.

Figure 7: The pairwise marginals of the ground-truth Mixture-of-Gaussians
distribution.

We then fitted a VAE model with 2-dimensional latent space using diagonal
Gaussian encoder and decoder distributions, and a fixed standard Normal prior.
For the decoder and encoder networks we used fully-connected residual neural
networks with 3 residual blocks, 200 hidden dimensions, and ReLU activations.
To optimise the model parameters we have used AMSGrad optimiser (Reddi et al.,
2018) with a learning rate of $10^{-3}$ for a total of 500 epochs.

The hyperameters are listed in table 3, note that the total number of samples
was the same for all methods (i.e. $5/15/25$). Moreover, we have used
“sticking-the-landing” (STL) gradients (Roeder et al., 2017) to reduce
gradient variance for all methods.111111We have also evaluated the doubly-
reparametrised gradients (DReG, Tucker et al., 2018) for IWAE methods but
found STL to perform similar or better.

Method | $Z$ | $K$ | $I$ | Mixture sampling
---|---|---|---|---
MVAE | $5/15/25$ | 1 | — | —
MissVAE | $5/15/25$ | $5/15/25$ | — | Ancestral
MissSVAE | 1 | $5/15/25$ | — | Stratified
MIWAE | 1 | 1 | $5/15/25$ | —
MissIWAE | 1 | $5/15/25$ | $5/15/25$ | Ancestral
MissSIWAE | 1 | $5/15/25$ | 1 | Stratified
DeMissVAE | 1 | $5/15/25$ | — |  LAIR (1 iteration, $R=0$) (Simkus & Gutmann, 2023)
Table 3: _Method hyperparameters on MoG data._

###  E.2 UCI data sets

We fit the VAEs on four data sets from the UCI repository (Dua & Graff, 2017)
with the preprocessing of (Papamakarios et al., 2017). The VAE model uses
diagonal Gaussian encoder and decoder distributions regularised such that the
standard deviation $\geq 10^{-5}$ (Mattei & Frellsen, 2018a), and a fixed
standard Normal prior. The latent space is 16-dimensional, except for the
MINIBOONE where 32 dimensions were used.

The encoder and decoder networks used fully-connected residual neural networks
with 2 residual blocks (except for on the MINIBOONE dataset where 5 blocks
were used in the encoder) with 256 hidden dimensionality, and ReLU
activations. A dropout of 0.5 was used on the MINIBOONE dataset. The
parameters were optimised using AMSGrad optimiser (Reddi et al., 2018) with a
learning rate of $10^{-3}$ and cosine learning rate schedule for a total of
200K iterations (or 22K iterations on MINIBOONE). As before, STL gradients
(Roeder et al., 2017) were used to reduce the variance for all methods.
DeMissVAE used the LAIR sampler (Simkus & Gutmann, 2023) with $K=5$ $R=1$ and
10 iterations. Moreover we have used gradient norm clipping to stabilise
DeMissVAE with the maximum norm set to $1$ (except for POWER dataset where we
set it to $0.5$).

###  E.3 MNIST and Omniglot data sets

We fit a VAE on statically binarised MNIST and Omniglot data sets (Lake et
al., 2015) downsampled to 28x28 pixels. The VAE uses diagonal Gaussian decoder
distributions regularised such that the standard deviation $\geq 10^{-5}$
(Mattei & Frellsen, 2018a), a fixed standard Normal prior, and a Bernoulli
decoder distribution. The latent space is 50-dimensional.

For both MNIST and Omniglot we have used convolutional ResNet neural networks
for the encoder and decoder with 4 residual blocks, ReLU activations, and
dropout probability of 0.3. For MNIST, the encoder the residual block hidden
dimensionalities were 32, 64, 128, 256, and for the decoder they were
128,64,32,32. For Omniglot, the encoder the residual block hidden
dimensionalities were 64, 128, 256, 512, and for the decoder they were
256,128,64,64. We used AMSGrad optimiser (Reddi et al., 2018) with $10^{-4}$
learning rate, cosine learning rate schedule, and STL gradients (Roeder et
al., 2017) for 500 epochs for MNIST and 200 epochs for Omniglot.

For MVAE, we use 5 latent samples and for MIWAE we use 5 importance samples.
For MissVAE we use $K=5$ mixture components and sample 5 latent samples. For
MissSVAE we use $K=5$ mixture components and sample 1 sample from each
component, for a total of 5 samples. For MissIWAE we use $K=5$ components and
sample 5 importance samples. for MissSIWAE we use $K=5$ components and sample
1 sample from each component. For DeMissVAE we use $K=5$ imputations and
update them using a single step of pseudo-Gibbs (Rezende et al., 2014).

##  Appendix F Additional figures

In this appendix we provide additional figures for the experiments in this
paper.

###  F.1 Mixture-of-Gaussians data with a 2D latent VAE

In this section we show additional analysis on the mixture-of-Gaussians data,
supplementing the results in section 6.1.

####  F.1.1 Analysis of gradient variance with ancestral and stratified
sampling

In section 6.1 we observed that the model estimation performance can depend on
whether ancestral sampling (with implicit reparametrisation) or stratified
sampling is used to approximate the expectations in eqs. 4, 5, 7 and 8,
corresponding to MissVAE/MissIWAE and MissSVAE/MissSIWAE, respectively.

We analyse the signal-to-noise ratio (SNR) of the gradients w.r.t.
${\bm{\phi}}$ and ${\bm{\theta}}$ for the two approaches, which is defined as
follows (Rainforth et al., 2019)

$\displaystyle\text{SNR}({\bm{\phi}})=\left|\frac{\mathbb{E}\left[\Delta({\bm{\phi}})\right]}{\sigma[\Delta({\bm{\phi}})]}\right|,\quad\text{and}\quad\text{SNR}({\bm{\theta}})=\left|\frac{\mathbb{E}\left[\Delta({\bm{\theta}})\right]}{\sigma[\Delta({\bm{\theta}})]}\right|,$

where $\Delta(\cdot)$ denotes the gradient estimate, and $\sigma[\cdot]$ is
the standard deviation of a random variable. We estimate the SNR by computing
the expectation and standard deviation over the entire training epoch.

The SNR for ${\bm{\phi}}$ and ${\bm{\theta}}$ is plotted in fig. 8. We observe
that the stratified approaches (MissSVAE and MissSIWAE) generally have higher
SNR. This is possibly the reason why MissSVAE and MissSIWAE have achieved
better model accuracy than the ancestral approaches (MissVAE and MissIWAE) in
section 6.1.

Figure 8: _Signal-to-noise ratio (SNR, higher is better) of the gradients
w.r.t. encoder parameters ${\bm{\phi}}$ (left) and decoder parameters
${\bm{\theta}}$ (right)._ For all methods we used a budget of 5 samples from
the variational distribution (see section E.1 for more details). We show the
median SNR over 5 independent runs and a 90% confidence interval.

####  F.1.2 Analysis of the model posteriors

In figs. 9, 10 and 11 we visualise the model posteriors with complete and
incomplete data, ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, respectively, and
the variational distribution ${q_{\bm{\phi}}}({\bm{z}}\mid\cdot)$ that was
used to fit the model via the variational ELBO. For each method we have used a
budget of 25 samples from the variational distribution during training
(additional details are in section E.1). Each figure shows the posteriors for
5 training data-points using distinct colours.

Figure 9 shows MVAE, MissVAE, and MissSVAE model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, as well as the
variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, which approximates the
incomplete-data posterior. As motivated in section 3 we observe that the
Gaussian posterior in MVAE (first row) is not sufficiently flexible to
approximate the complex incomplete-data posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$. On the other hand,
the mixture-variational approaches, MissVAE (second row) and MissSVAE (third
row), are able to well-approximate the incomplete-data posteriors.

Figure 10 shows MIWAE, MissIWAE, and MissSIWAE model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, as well as the
variational proposal ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
and the importance-weighted semi-implicit distribution
$q_{{\bm{\phi}},{\bm{\theta}},I=25}^{\text{IW}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
that arises from sampling the variational proposal
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ and re-sampling using
importance weights
$w({\bm{z}})={p_{\bm{\theta}}}({\bm{x}}_{\mathrm{obs}},{\bm{z}})/{q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$
(Cremer et al., 2017). Similar to the MVAE case above, the variational
proposal ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ in MIWAE
(first row) is quite far from the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, but the importance-
weighted bound in eq. 7 is able to re-weigh the samples to sufficiently-well
match the model posterior, as shown in the fourth column. However, as
efficiency of importance sampling depends on the discrepancy between the
proposal and the target distributions, we can expect that more flexible
variational distributions may improve the performance of the importance-
weighted ELBO methods. Importantly, we show that the variational-mixture
approaches, MissIWAE (second row) and MissSIWAE (third row), are able to adapt
the variational proposals to the incomplete-data posteriors well, and as a
result achieve better efficiency than MIWAE.

Figure 11 shows DeMissVAE model posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ and
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$, the variational
distribution ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$, which approximates the
_complete_ -data posterior, and the imputation-mixture
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ approximated
using the 25 imputations in
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ at the end of
training. We observe similar behaviour to fig. 1, where the complete data
posteriors ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ are close to Gaussian but
the incomplete-data posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ are irregular. As we
show in section 6.1, DeMissVAE is capable of fitting the model well by
learning the completed-data variational distribution
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ (third column) and using the
imputation-mixture in eq. 9 to approximate the incomplete data posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$. Moreover, we observe
that the imputation-mixture
${q_{{\bm{\phi}},f^{t}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ (fourth column)
captures only one of the modes of the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$. This is a result of
using a small imputation sampling budget, that is, using only a single
iteration of LAIR to update the imputations (see more details in appendix D),
and hence better accuracy can be achieved by trading-off some computational
cost to obtain better imputations that would ensure a better representation of
the imputation distribution. Nonetheless, as observed in fig. 2, DeMissVAE
achieves good model accuracy despite potentially sub-optimal imputations,
further signifying the importance of the two learning objectives for the
encoder and decoder distributions in sections 4.2 and C.

Interestingly, by comparing the complete-data posteriors
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ (first column) in figs. 9, 10 and
11, we observe that they are slightly more irregular than in the complete case
in fig. 1, except for DeMissVAE whose posteriors are nearly Gaussian. The
irregularity is stronger for the importance-weighted ELBO-based methods in
fig. 10. This is in line with the observation by Burda et al. (2015, Appendix
C) and Cremer et al. (2018, Section 5.4) that VAEs trained with more flexible
variational distributions tend to learn a more complex model posterior. This
means that using the importance-weighted bounds, and to a lesser extent the
finite variational-mixture approaches from section 4.1, to fit VAEs on
incomplete data may result in worse-structured latent spaces, compared to
models fitted on complete data. On the other hand, we observe that DeMissVAE
learns a better-structured latent space, with the posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ close to a Gaussian, that is
comparable to the complete case. This suggests that the decomposed approach in
DeMissVAE may be important in cases where the latent space needs to be
regular, at the additional cost of obtaining missing data imputations (see
appendix D).

Figure 9: _Posterior distributions of MVAE, MissVAE, and MissSVAE._ First
column: the model posterior ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ under
complete data ${\bm{x}}$. Second column: the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ under incomplete data
${\bm{x}}_{\mathrm{obs}}$. Third column: variational approximation
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ of the incomplete
posterior ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ obtained at
the end of training.  Figure 10: _Posterior distributions of MIWAE, MissIWAE,
and MissSIWAE._ First column: the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ under complete data ${\bm{x}}$.
Second column: the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ under incomplete data
${\bm{x}}_{\mathrm{obs}}$. Third column: variational proposal
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ for an incomplete data-
point ${\bm{x}}_{\mathrm{obs}}$ obtained at the end of training. Fourth
column: importance-weighted variational distribution
${q_{\bm{\phi}}}^{\text{IW}}_{I=25}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ for
an incomplete data-point ${\bm{x}}_{\mathrm{obs}}$ obtained after re-weighting
samples from ${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ (Cremer et
al., 2017). Figure 11: _Posterior distributions of DeMissVAE._ First: the
model posterior ${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ under complete data
${\bm{x}}$. Second: the model posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}}_{\mathrm{obs}})$ under incomplete data
${\bm{x}}_{\mathrm{obs}}$. Third: variational approximation
${q_{\bm{\phi}}}({\bm{z}}\mid{\bm{x}})$ of the complete-data posterior
${p_{\bm{\theta}}}({\bm{z}}\mid{\bm{x}})$ obtained at the end of training.
Fourth: the variational imputation-mixture distribution in eq. 9 using the
imputation distribution
${f^{t}({\bm{x}}_{\mathrm{mis}}\mid{\bm{x}}_{\mathrm{obs}})}$ obtained at the
end of training, approximated using a Monte Carlo average with 25 imputations.

###  F.2 UCI data sets

In fig. 12 we plot the Fréchet inception distance (FID, Heusel et al., 2017)
versus training iteration on the UCI datasets. The results closely mimic the
log-likelihood results in section 6.2. Importantly, we observe that using
mixture variational distributions becomes more important as the missingness
fraction increases, causing the posterior distributions to be more complex.

Figure 12: _FID (lower is better) between the model and the complete test data
versus training iterations._ The FID is computed using features of the last
encoder layer of an independent VAE model trained on complete data. Lines show
the median of 5 independent runs and the intervals show 90% confidence.