# Leveraging the Exact Likelihood of Deep Latent Variables Models

Deep latent variable models combine the approximation abilities of deep neural networks and the statistical foundations of generative models. The induced data distribution is an infinite mixture model whose density is extremely delicate to compute. Variational methods are consequently used for inference, following the seminal work of Rezende et al. (2014) and Kingma and Welling (2014). We study the well-posedness of the exact problem (maximum likelihood) these techniques approximatively solve. In particular, we show that most unconstrained models used for continuous data have an unbounded likelihood. This ill-posedness and the problems it causes are illustrated on real data. We also show how to insure the existence of maximum likelihood estimates, and draw useful connections with nonparametric mixture models. Furthermore, we describe an algorithm that allows to perform missing data imputation using the exact conditional likelihood of a deep latent variable model. On several real data sets, our algorithm consistently and significantly outperforms the usual imputation scheme used within deep latent variable models.

• 15 publications
• 20 publications
02/13/2018

### Leveraging the Exact Likelihood of Deep Latent Variable Models

Deep latent variable models combine the approximation abilities of deep ...
01/11/2018

### Autoencoders and Probabilistic Inference with Missing Data: An Exact Solution for The Factor Analysis Case

Latent variable models can be used to probabilistically "fill-in" missin...
07/20/2022

### Maximum Likelihood Imputation

Maximum likelihood (ML) estimation is widely used in statistics. The h-l...
10/04/2021

### A moment-matching metric for latent variable generative models

It can be difficult to assess the quality of a fitted model when facing ...
08/22/2021

### Efficient Gaussian Neural Processes for Regression

Conditional Neural Processes (CNP; Garnelo et al., 2018) are an attracti...
11/18/2020

### Detecting Hierarchical Changes in Latent Variable Models

This paper addresses the issue of detecting hierarchical changes in late...
12/22/2021

### Density Regression with Bayesian Additive Regression Trees

Flexibly modeling how an entire density changes with covariates is an im...

## 1 Introduction

Dimension reduction aims at summarizing multivariate data using a small number of features that constitute a code. Earliest attempts rested on linear projections, leading to Hotelling’s (hotelling1933) principal component analysis (PCA) that has been vastly explored and perfected over the last century (jolliffe2016principal). In recent years, the field has been vividly animated by the successes of latent variable models, that probabilistically use the low-dimensional features to define powerful generative models. Usually, these latent variable models transform the random code into parameters of a simple output distribution. Linear mappings and Gaussian outputs were initially considered, giving rise to factor analysis (bartholomew2011)

and probabilistic principal component analysis

(tipping1999)

. Beyond the Gaussian distribution, Bernoulli or more general exponential family outputs have been considered. In recent years, much work has been done regarding nonlinear mappings parametrised by deep neural networks, following the seminal papers of

rezende2014 and kingma2014. These models have led to impressive empirical performance in unsupervised and semi-supervised generative modelling of images (narayanaswamy2017), molecular structures (kusner2017; gomez2018), or arithmetic expressions (kusner2017). This paper is an investigation of the statistical properties of these models, which remain essentially unknown.

### 1.1 Deep latent variable models

In their most common form, deep latent variable models (DLVMs) assume that we are in presence of a data matrix that we wish to explain using some latent variables . We assume that

are independent and identically distributed (i.i.d.) random variables driven by the following generative model:

 {z∼p(z)pθ(x|z)=Φ(x|fθ(z)). (1)

The unobserved random vector

is called the latent variable and usually follows marginally a simple distribution called the prior distribution. The dimension of the latent space is called the intrinsic dimension, and is usually smaller than the dimensionality of the data. The collection is a parametric family of output densities with respect to a dominating measure (usually the Lebesgue or the counting measure). The function is called a decoder or a generative network, and is parametrised by a (deep) neural network whose weights are stored in . This parametrisation allows to leverage recent advances in deep architectures, such as deep residual networks (kingma2016), recurrent networks bowman2016, or batch normalisation (sonderby2016). Several output families have been considered:

• [leftmargin=3ex,topsep=0ex]

• In case of discrete multivariate data, products of Multinomial distributions.

• In case of real multivariate data, multivariate Gaussian distributions.

• When dealing with images, several specific proposals have been made (e.g. the discretised logistic mixture approach of salimans2017).

• Dirac outputs correspond to deterministic decoders, that are used e.g. within generative adversarial networks (goodfellow2014), or non-volume preserving transformations (dinh2017).

The Gaussian and Bernoulli families are the most widely studied, and will be the focus of this article. The latent structure of these DLVMs leads to the following marginal distribution of the data:

 pθ(x)=∫Rdpθ(x|z)p(z)dz=∫RdΦ(x|fθ(z))p(z)dz. (2)

### 1.2 Scalable learning through amortised variational inference

The log-likelihood function of a DLVM is, for all ,

 ℓ(θ)=logpθ(X)=n∑i=1logpθ(xi), (3)

which is an extremely challenging quantity to compute that involves potentially high-dimensional integrals. Estimating by maximum likelihood appears therefore out of reach. Consequently, following rezende2014 and kingma2014, inference in DLVMs is usually performed using amortised variational inference. Variational inference approximatively maximises the log-likelihood by maximising a lower bound known as the evidence lower bound (ELBO, see e.g. blei2017):

 ELBO(θ,q)=EZ∼q[logp(X,Z)q(Z)]=ℓ(θ)−KL(q||p(⋅|X))≤ℓ(θ). (4)

where the variational distribution is a distribution over the space of codes . The variational distribution plays the role of a tractable approximation of the posterior distribution of the codes; when this approximation is perfectly accurate, the ELBO is equal to the log-likelihood. Amortised inference builds using a neural network called the inference network , whose weights are stored in :

 qγ,X(Z)=qγ,X(z1,...,zn)=n∏i=1Ψ(zi|gγ(xi)), (5)

where is a parametric family of distributions over – such as Gaussians with diagonal covariances (kingma2014). Other kinds of families – built using e.g. normalising flows (rezende2015; kingma2016), auxiliary variables (maaloe2016; ranganath2016), or importance weights (burda2016; cremer2017) – have been considered for amortised inference, but they will not be central focus of in this paper. Within this framework, variational inference for DLVMs solves the optimisation problem using variants of stochastic gradient ascent (see e.g. roeder2017, for strategies for computing gradients estimates of the ELBO).

As emphasised by kingma2014

, the ELBO resembles the objective function of a popular deep learning model called an

autoencoder (see e.g. goodfellow2016, Chapter 14). This motivates the popular denomination of encoder for the inference network and variational autoencoder (VAE) for the combination of a DLVM with amortised variational inference.

#### Contributions.

In this work, we revisit DLVMs by asking: Is it possible to leverage the properties of to understand and improve deep generative modelling? Our main contributions are:

• [leftmargin=3ex,topsep=0ex]

• We show that maximum likelihood is ill-posed for continuous DLVMs and well-posed for discrete ones. We link this undesirable property of continuous DLVMs to the mode collapse phenomenon, and illustrate it on a real data set.

• We draw a connection between DLVMs and nonparametric statistics, and show that DLVMs can be seen as parsimonious submodels of nonparametric mixture models.

• We leverage this connection to provide a way of finding an upper bound of the likelihood based on finite mixtures. Combined with the ELBO, this bound allows us to provide useful “sandwichings” of the exact likelihood. We also prove that this bound characterises the large capacity behaviour of DLVMs.

• When dealing with missing data, we show how a simple modification of an approximate scheme proposed by rezende2014 allows us to draw according to the exact conditional distribution of the missing data. On several data sets and missing data scenarios, our algorithm consistently outperforms the one of rezende2014, while having the same computational cost.

## 2 Is maximum likelihood well-defined for deep latent variable models?

In this section, we investigate the properties of maximum likelihood estimation for DLVMs with Gaussian and Bernoulli outputs.

### 2.1 On the boundedness of the likelihood of deep latent variable models

Deep generative models with Gaussian outputs assume that the data space is , and that the family of output distributions is the family of -variate full-rank Gaussian distributions. The conditional distribution of each data point is consequently

 pθ(x|z)=N(x|μθ(z),Σθ(z)), (6)

where and are two continuous functions parametrised by neural networks whose weights are stored in a parameter . These two functions constitute the decoder of the model. This leads to the log-likelihood

 ℓ(θ)=n∑i=1log(∫RdN(xi|μθ(z),Σθ(z))p(z)dz). (7)

This model can be seen as a special case of infinite mixture of Gaussian distributions. However, it is well-known that maximum likelihood is ill-posed for finite Gaussian mixtures (see e.g. lecam1990). Indeed, the likelihood function is unbounded above, and its infinite maxima happen to be very poor generative models. This problematic behaviour of a model quite similar to DLVMs motivates the question: is the likelihood function of DLVMs bounded above?

In this section, we will not make any particular parametric assumption about the prior distribution of the latent variable . Both simple isotropic Gaussian distributions and more complex learnable priors have been proposed, but all of them are affected by the negative results of this section. We simply make the natural assumptions that is continuous and has zero mean. Many different neural architectures have been explored regarding the parametrisation of the decoder. For example, kingma2014

consider multilayer perceptrons (MLPs) of the form

 μθ(z)=Vtanh(Wz+a)+b,Σθ(z)=Diag(exp(αtanh(Wz+a)+β)), (8)

where . The weights of the decoder are , and . The integer is the (common) number of hidden units of the MLPs. Much more complex parametrisations exist, but we will see that this one, arguably one of the most rigid, is already too flexible for maximum likelihood. Actually, we will show that an even much less flexible family of MLPs with a single hidden unit is problematic. Let and let be a sequence of nonnegative real numbers such that as . For all , with the parametrisation (8), we consider the sequences of parameters . This leads to the following simplified decoders:

 (9)

As shown by next theorem, these sequences of decoders lead to the divergence of the log-likelihood function.

###### Theorem 1.

For all and , we have

###### Proof.

A detailed proof is provided in Appendix A. Its cornertone is the fact that the sequence of functions converges to a function that outputs both singular and nonsingular covariances, leading to the explosion of while all other terms in the likelihood remain bounded below by a constant. ∎

Using simple MLP-based parametrisations such a the one of kingma2014 therefore brings about an unbounded log-likelihood function. A natural question that follows is: do these infinite suprema lead to useful generative models? The answer is no. Actually, none of the functions considered in Theorem 1 are particularly useful, because of the use of a constant mean function. This is formalised in the next proposition, that exhibits a strong link between likelihood blow-up and the mode collapse phenomenon.

###### Proposition 1.

For all , , and , the distribution is spherically symmetric and unimodal around .

###### Proof.

Because of the constant mean function and the isotropic covariance, the density of is a decreasing function , hence the spherical symmetry and the unimodality. Note that there are several different definitions for multivariate unimodality (see e.g. dharmadhikari1988). Here, we mean that the only local maximum of the density is at .∎

The spherical symmetry implies that the distribution of these “optimal” deep generative model will lead to uncorrelated variables, and the unimodality will lead to poor sample diversity. This behaviour is symptomatic of mode collapse, which remains one of the most challenging drawbacks of generative modelling (arora2018). Unregularised gradient-based optimisation of a tight lower bound of this unbounded likelihood is therefore likely to chase these (infinitely many) bad suprema. This gives a theoretical foundation to the necessary regularisation of VAEs that was already noted by rezende2014 and kingma2014. For example, using weight decay as in kingma2014

is likely to help avoiding these infinite maxima. This difficulty to learn the variance may also explain the choice made by several authors to use a constant variance function of the form

, where can be either fixed (zhao2017) or learned via approximate maximum likelihood (pu2016).

#### Tackling the unboundedness of the likelihood.

Let us go back to a parametrisation which is not necessarily MLP-based. Even in this general context, it is possible to tackle the unboundedness of the likelihood using additional constraints on . Specifically, for each , we will consider the set Note that . This simple spectral constraint allows to end up with a bounded likelihood.

###### Proposition 2.

Let . If the parametrisation of the decoder is such that the image of is included in for all , then the log-likelihood function is upper bounded by .

###### Proof.

For all , we have

 p(xi|μθ,Σθ)=∫RdN(x|μθ(z),Σθ(z))N(z|0,Id)dz≤1(2πξ)p/2,

using the fact that the determinant of is lower bounded by for all and that the exponential of a negative number is smaller than one. Therefore, the likelihood function is bounded above by . ∎

Similar constraints have been proposed to solve the ill-posedness of maximum likelihood for finite Gaussian mixtures (e.g. hathaway1985; biernacki2011). In practice, implementing such constraints can be easily done by adding a constant diagonal matrix to the output of the covariance decoder.

We chose a specific and natural parametrisation in order to obtain a constructive proof of the unboundedness of the likelihood. However, virtually any other deep parametrisation that does not include covariance constraints will be affected by our result, because of the universal approximation abilities of neural networks (see e.g. goodfellow2016, Section 6.4.1).

#### Bernoulli DLVMs do not suffer from unbounded likelihood.

When , Bernoulli DLVMs assume that is the family of

-variate multivariate Bernoulli distributions (i.e. the family of products of

univariate Bernoulli distributions). In this case, maximum likelihood is well-posed.

###### Proposition 3.

Given any possible parametrisation, the log-likelihood function of a deep latent model with Bernoulli outputs is everywhere negative.

###### Proof.

This is a direct consequence of the fact that the density of a Bernoulli distribution is always smaller than one. ∎

### 2.2 Towards data-dependent likelihood upper bounds

We have determined under which conditions maximum likelihood estimates exist, and have computed simple upper bounds on the likelihood functions. Since they do not depend on the data, these bounds are likely to be very loose. A natural follow-up issue is to seek tighter, data-dependent upper bounds that remain easily computable. Such bounds are desirable because, combined with ELBOs, they would allow sandwiching the likelihood between two bounds.

To study this problem, let us take a step backwards and consider a more general infinite mixture model. Precisely, given any distribution over the generic parameter space , we define the nonparametric mixture model (see e.g. lindsay1995, Chapter 1) as:

 pG(x)=∫HΦ(x|η)dG(η). (10)

Note that there are many ways for a mixture model to be nonparametric (e.g. having some nonparametric components, an infinite but countable number of components, or un uncountable number of components). In this case, this comes from the fact that the model parameter is the mixing distribution , which belongs to the set

of all probability measures over

. The likelihood of any is given by

 ℓ(G)=n∑i=1logpG(xi). (11)

When has a finite support of cardinal , is a finite mixture model with components. When the mixing distribution is generatively defined by the distribution of a random variable such that we exactly end up with a deep generative model with decoder . Therefore, the nonparametric mixture is a more general model that the DLVM. The fact that the mixing distribution of a DLVM is intrinsically low-dimensional leads us to interpret the DLVM as a parsimonious submodel of the nonparametric mixture model. This also gives us an immediate upper bound on the likelihood of any decoder :

Of course, in many cases, this upper bound will be infinite (for example in the case of unconstrained Gaussian outputs). However, under the conditions of boundedness of the likelihood of deep Gaussian models, the bound is finite and attained for a finite mixture model with at most components.

###### Theorem 2.

Assume that is the family of multivariate Bernoulli distributions or the family of Gaussian distributions with the spectral constraint of Proposition 2. The likelihood of the corresponding nonparametric mixture model is maximised for a finite mixture model of distributions from the family .

###### Proof.

A detailed proof is provided in Appendix B. The main tool of the proof of this rather surprising result is Lindsay’s (lindsay1983) geometric analysis of the likelihood of nonparametric mixtures, based on Minkovski’s theorem. Specifically, Lindsay’s (lindsay1983) Theorem 3.1 ensures that, when the trace of the likelihood curve is compact, the likelihood function is maximised for a finite mixture. For the Bernoulli case, compactness of the curve is immediate; for the Gaussian case, we use a compactification argument inspired by van1992. ∎

Assume now that the conditions of Theorem 2 are satisfied. Let us denote a maximum likelihood estimate of as . For all , we therefore have

 ℓ(θ)≤ℓ(^G), (12)

which gives an upper bound on the likelihood. We call the difference the parsimony gap (see Fig. 1). By sandwiching the exact likelihood between this bound and an ELBO, we can also have guarantees on how far a posterior approximation is from the true posterior:

 KL(q||p(⋅|X))≤ℓ(^G)−ELBO(θ,q). (13)

From a computational perspective, the estimate can be found using the expectation-maximisation algorithm for finite mixtures (dempster1977) – although it only ensures to find a local optimum. Some strategies guaranteed to find a global optimum have also been developed (e.g. lindsay1995, Chapter 6, or wang2007).

Now that computationally approachable upper bounds have been derived, the question remains whether or not these bounds can be tight. Actually, as shown by next theorem, tightness of the parsimony gap occurs when the decoder has universal approximation abilities. In other words, the nonparametric upper bound characterises the large capacity limit of the decoder.

###### Theorem 3 (Tightness of the parsimony gap).

Assume that

1. is the family of multivariate Bernoulli distributions or the family of Gaussian distributions with the spectral constraint of Proposition 2,

2. The decoder has universal approximation abilities: for any compact and continuous function , for all , there exists such that .

Then, for all , there exists such that .

###### Proof.

A detailed proof is provided in Appendix C. The main idea is to split the code space into a compact set made of several parts that will represent the mixture components, and an unbounded set of very small prior mass. The universal approximation property is finally used for this compact set. ∎

The universal approximation condition is satisfied for example by MLPs with nonpolynomial activations (leshno1993). Combined with the work of cremer2018, who studied the large capacity limit of the encoder, this result describes the general behaviour of a VAE in the large capacity limit (see Fig. 1).

## 3 Missing data imputation using the exact conditional likelihood

In this section, we assume that a variational autoencoder has been trained, and that some data is missing at test time. The couple decoder/encoder obtained after training is denoted by and . Let be a new data point that consists in some observed features and missing data . Since we have a probabilistic model of the data, an ideal way of imputing would be to generate some data according to the conditional likelihood

 pθ(xmiss|xobs)=∫Rdpθ(xmiss|xobs,z)p(z|xobs)dz. (14)

Again, this distribution appears out of reach because of the integration of the latent variable . However, it is reasonable to assume that, for all , is easily computable (this is the case for Bernoulli or Gaussian outputs). Under this assumption, we will see that generating data according to the conditional distribution is actually possible.

### 3.1 Pseudo-Gibbs sampling

rezende2014 proposed a simple way of imputing

by following a Markov chain

(initialised by randomly imputing the missing data with ). For all , the chain alternatively generates and until convergence. This scheme closely resembles Gibbs sampling (Geman1984), and actually exactly coincides with Gibbs sampling when the amortised variational distribution is equal to the true posterior distribution for all possible . Following the terminology of heckerman2000, we will call this algorithm pseudo-Gibbs sampling. Very similar schemes have been proposed for more general autoencoder settings (goodfellow2016, Section 20.11). Because of its flexibility, this pseudo-Gibbs approach is routinely used for missing data imputation using DLVMs (see e.g. li2016; rezende2016; du2018). rezende2014 proved that, when these two distributions are close in some sense, pseudo-Gibbs sampling generates points that approximatively follow the conditional distribution . Actually, we will see that a simple modification of this scheme allows to generate exactly according to the conditional distribution.

### 3.2 Metropolis-within-Gibbs sampling

At each step of the chain, rather than generating codes according to the approximate posterior distribution, we may use this approximation as a proposal within a Metropolis-Hastings algorithm (Metropolis1953; Hastings1970), using the fact that we have access to the unnormalised posterior density of the latent codes.

Specifically, at each step, we will generate a new code as a proposal using the approximate posterior . This proposal is kept as a valid code with acceptance probability

 ρt=Φ(xobs,^xmisst−1|fθ(~zt))p(~zt)Φ(xobs,^xmisst−1|fθ(zt−1))p(zt−1)Ψ(zt−1|gγ(xobs,^xmisst−1))Ψ(~zt|gγ(xobs,^xmisst−1)). (15)

This probability corresponds to a ratio of importance ratios, and is equal to one when the posterior approximation is perfect. This code-generating scheme exactly corresponds to performing a single iteration of an independent Metropolis-Hastings algorithm. With the obtained code , we can now generate a new imputation using the exact conditional . The obtained algorithm, detailed in Algorithm 1, is a particular instance of a Metropolis-within-Gibbs algorithm. Actually, it exactly corresponds to the algorithm described by gelman1993, and is ensured to asymptotically produce samples from the true conditional ditribution , even if the variational approximation is imperfect. Note that when the variational approximation is perfect, all proposals are accepted and the algorithm exactly reduces to Gibbs sampling.

The theoretical superiority of the Metropolis-within-Gibbs scheme compared to the pseudo-Gibbs sampler comes with almost no additional computational cost

. Indeed, all the quantities that need to be computed in order to compute the acceptance probability need also to be computed within the pseudo-Gibbs scheme – except for prior evaluations, which are assumed to be computationally negligible. However, a poor initialisation of the missing values might lead to a lot of rejections at the beginning of the chain, and to slow convergence. A good initialisation heuristic is to perform a few pseudo-Gibbs iterations at first in order to begin with a sensible imputation. Note eventually that, similarly to the pseudo-Gibbs sampler, our Metropolis-within-Gibbs scheme can be extended to many other variational approximations – like normalising flows

(rezende2015; kingma2016) – in a straightforward manner.

## 4 Empirical results

In this section, we investigate the empirical realisations of our theoretical findings on DLVMs.

#### Implementation.

Some computational choices are common for all experiments: the prior distribution is a standard Gaussian distribution; the chosen variational family is the family of Gaussian distributions with “diagonal + rank-1” covariance (as in rezende2014, Section 4.3); stochastic gradients of the ELBO are computed via the path derivative estimate of roeder2017; the Adam optimiser (kingma2014adam) is used with a learning rate of and mini-batches of size 10; neural network are initialised following the heuristics of glorot2010

; sampling for variational autoencoders is performed via the Distributions module of TensorFlow

(dillon2017). For the Frey faces data set, we used a Gaussian DLVM together with the parametrisation presented in Section 2.1 (with and ). The architecture of the inference network follows the same parametrisation. Constrained finite Gaussian mixtures were fit using the scikit-learn package (pedregosa2011). Regarding the data sets considered in the missing data experiments, we used the architecture of rezende2014, with hidden units and an intrinsic dimension of .

### 4.1 Witnessing likelihood blow-up

To investigate if the unboundedness of the likelihood of a DLVM with Gaussian outputs has actual concrete consequences for variational inference, we train two DLVMs on the Frey faces data set: one with no constraints, and one with the constraint of Proposition 2 (with ). The results are presented in Fig. 2. One can notice that the unconstrained DLVM finds models with very high likelihood but very poor generalisation performance. This confirms that the unboundedness of the likelihood is not a merely theoretical concern. We also display the two upper bounds of the likelihood. The nonparametric bound offers a slight but significant improvement over the naive upper bound. On this example, using the nonparametric upper bound as an early stopping criterion for the unconstrained ELBO appears to provide a good regularisation scheme – that perform better than the covariance constraints on this data set. This illustrates the potential practical usefulness of the connection that we drew between DLVMs and nonparametric mixtures.

### 4.2 Comparing the pseudo-Gibbs and Metropolis-within-Gibbs samplers

We compare the two samplers for single imputation of the test sets of three data sets: Caltech 101 Silhouettes and statically binarised versions of MNIST and OMNIGLOT. We consider two missing data scenarios: a first one with pixels missing uniformly at random (the fractions of missing data considered are , and ) and one where the top or bottom half of the pixels was removed. Both samplers use the same trained VAE and perform the same number of iterations ( for the first scenario, and for the second). Note that the convergence can be monitored using a validation set of complete data. The first iterations of the Metropolis-within-Gibbs sampler are actually pseudo-Gibbs iterations, in order to avoid having too many early rejections. The chains converge much faster for the missing at random (MAR) situation than for the top/bottom missing scenario. This is probably due to the fact that the conditional distribution of the missing half of an image is highly multimodal. The results are displayed on Fig. 4 and in the supplementary material. The Metropolis-within-Gibbs sampler consistently outperforms the pseudo-Gibbs scheme, especially for the most challenging scenarios where the top/bottom of the image is missing. Moreover, according to Wilcoxon signed-rank tests (see e.g. dalgaard2008, Chapter 5), the Metropolis-within-Gibbs scheme is significantly more accurate than the pseudo-Gibbs sampler in all scenarios (all -values are provided as supplementary material). One can see that the pseudo-Gibbs sampler appears to converge quickly to a stationary distribution that gives suboptimal results. The Metropolis-within-Gibbs algorithm converges slower, but to a much more useful distribution.

## 5 Conclusion

Although extremely difficult to compute in practice, the exact likelihood of DLVMs offers several important insights on deep generative modelling. An important research direction for future work is the design of principled regularisation schemes for maximum likelihood estimation.

The objective evaluation of deep generative models remains a widely open question. Missing data imputation is often used as a performance metric for DLVMs (e.g. li2016; du2018). Since both algorithms have essentially the same computational cost, this motivates to replace pseudo-Gibbs sampling by Metropolis-within-Gibbs when evaluating these models. Upon convergence, the samples generated by Metropolis-within-Gibbs do not depend on the inference network, and explicitly depend on the prior, which allows us to evaluate mainly the generative performance of the models.

We interpreted DLVMs as parsimonious submodels of nonparametric mixture models. While we used this connection to provide upper bounds of the likelihood, many other applications could be derived. In particular, the important body of work regarding consistency of maximum likelihood estimates for nonparametric mixtures (e.g. kiefer1956; van2003; chen2017) could be leveraged to study the asymptotics of DLVMs.

## Appendix A. Proof of Theorem 1

###### Theorem 1.

For all and , we have

###### Proof.

Let and . To avoid notational overload, we denote in the remainder of this proof. We will show that the contribution of the -th observation explodes while all other contributions remain bounded below.

A first useful remark is the fact that, since is continuous and has zero mean and , the univariate random variable is continuous and has zero mean. Therefore and .

Regarding the -th observation, we have for all ,

 pθk(xi) =∫RdN(xi|xi,Σθk(z))p(z)dz (16) ≥∫wTz≤0N(xi|xi,Σθk(z))p(z)dz (17) ≥∫wTz≤0|2πΣθk(z)|−1/2p(z)dz. (18)

Let such that . The function

 φ:α↦αtanh(αwTz)−α,

is strictly decreasing on . Indeed, its derivative is equal to

 φ′(α)=tanh(αwTz)−1+(1−tanh2(αwTz))αwT, (19)

which is strictly negative because the image of the hyperbolic tangent function is . Moreover,

 limα→∞αtanh(αwTz)−α=−∞.

Therefore, the sequence is strictly increasing and diverges to for all such that . Therefore , the monotone convergence theorem combined with the fact that insure that the right hand side of (16) diverges to , leading to .

Regarding the other contributions, let . Since , there exists such that . We have

 pθk(xj) ≥∫wTz≥εN(xj|xi,Σθk(z))p(z)dz =∫wTz≥εexp(−||xj−xi||222exp(αktanh(αkwTz)−αk))(2πexp(αktanh(αkwTz)−αk))p/2p(z)dz =1(2π)p/2∫wTz≥εexp(−||xj−xi||222exp(αktanh(αkwTz)−αk))p(z)dz ≥1(2π)p/2exp(−||xj−xi||222exp(αktanh(αkε)−αk))P(wTz≥ε),

and, since , we will have

 limsupk→∞pθk(xj)≥P(wTz)(2π)p/2exp(−||xj−xi||222)≥ε,

therefore

 limsupk→∞pθk(xj)>0.

By combining all contributions, we end up with . ∎

## Appendix B. Proof of Theorem 2

###### Theorem 2.

Assume that is the family of multivariate Bernoulli distributions or the family of Gaussian distributions with the spectral constraint of Proposition 2. The likelihood of the corresponding nonparametric mixture model is maximised for a finite mixture model of distributions from the family .

###### Proof.

Let us assume that there are distinct observations . Following lindsay1983, let us consider the trace of the likelihood curve

 Γ={(Φ(~xl|η))l≤L|η∈H}.

We will use the following theorem, which describes maximum likelihood estimators of nonparametric mixtures under some topological conditions.

###### Theorem (lindsay1983, Theorem 3.1).

If is compact, then the likelihood of the corresponding nonparametric mixture model is maximised for a finite mixture model of distributions from the family .

This theorem allows us to treat both cases:

#### Bernoulli outputs.

Since is compact and the function is continuous, is compact and Lindsay’s theorem can be directly applied.

#### Gaussian outputs.

The parameter space is not compact in this case, but we can get around this problem using a compactification argument similar to the one of van1992. Consider the Alexandroff compactification of the parameter space (kelley1955, p. 150). Because of the definition of , we have for all . Therefore, we can continuously extend the function from to using the conventions for all . The space is therefore compact, and we can use Lindsay’s theorem to deduce that the nonparametric maximum likelihood estimator for the compactified parameter space is a finite mixture of distributions from . However, this finite mixture can only contain elements from . Indeed, if any mixture component were associated with , the likelihood could be improved by emptying said component. Therefore, the maximum likelihood estimator found using the compactified space is also the maximum likelihood estimator of the original space, which allows us to conclude. ∎

## Appendix C. Proof of Theorem 3

###### Theorem 3.

Assume that

1. is the family of multivariate Bernoulli distributions or the family of Gaussian distributions with the spectral constraint of Proposition 2.

2. The decoder has universal approximation abilities : for any compact and continuous function , for all , there exists such that .

Then, for all , there exists such that

 ℓ(^G)≥ℓ(θ)≥ℓ(^G)−ε.
###### Proof.

The left hand side of the inequality directly comes from the first assumption and Theorem 2. We will now prove the right hand side.

Let . Since the logarithm is continuous, there exists such that

 0≤u≤δ⟹n∑i=1log(p^G(xi)−u)≥ℓ(^G)−ε. (20)

We will show that and can be made closer than for all .

Theorem 2 insures that is a finite mixture, so there exist some and such that

 ∀x∈X,p^G(x)=K∑k=1πkΦ(x|ηk).

We will partition the domain of integration of the latent variable into several parts:

• an infinite part of very small prior mass,

• a compact set of high prior mass over which we can apply the universal approximation property.

The compact set will be divided itself into parts: one part for each mixture component and small parts to insure the continuity of the decoder. More specifically, let , let

be the (prior) cumulative distribution function of

, and let:

Let . Consider a continuous function such that

 ∀z∈K,f(z)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩η1ifα1≤||z||2≤β1η2ifα2≤||z||2≤β2...ηKifαK≤||z||2≤βK.

For example, can be built using the Tietze extension theorem (see e.g. kelley1955, p. 242) together with the conditions of the above formula. According to the universal approximation property, we can find decoders arbitrarily close to . How close do they need to be? Invoking the continuity of the functions in , there exists such that

 ||ηk−η||∞≤Δ⟹|Φ(xi|η)−Φ(xi|ηk)|≤e.

We will consider now a decoder such that .

We can write, for all ,

 pθ(xi) =K∑l=1(∫αl≤||z||2≤βlΦ(xi|fθ(z))p(z)dz+∫βl≤||z||2≤αl+1Φ(xi|fθ(z))p(z)dz) ≥K∑l=1∫αl≤||z||2≤βlΦ(xi|fθ(z))p(z)dz.

Using the fact that , we have therefore, for all ,

 pθ(xi) ≥K∑l=1∫αl≤||z||2≤βl(Φ(xi|ηk)−e)p(z)dz =K∑l=1(Φ(xi|ηk)−e)(πk−e) =K∑l=1(πkΦ(xi|ηk)+e2−e(πk+Φ(xi|ηk))) ≥K∑l=1πkΦ(xi|ηk)−e(K∑l=1Φ(xi|ηk)+1).

Now, if we take such that

 e(K∑l=1Φ(xi|ηk)+1)≤δ,

for all , we end up with

 pθ(xi)≥p^G(xi)−δ.

Using (20), this eventually leads to . ∎

## Appendix G. Additional imputation experiments

We provide below all -values obtained for Wilcoxon signed-rank tests (see e.g. dalgaard2008, Chapter 5) based on the accuracy of the two imputation schemes over all images in the test sets. The Metropolis-within-Gibbs sampler always significantly outperforms the pseudo-Gibbs scheme.