1 Introduction
Many core machine learning tasks are concerned with density estimation and manifold discovery. Probabilistic graphical models are a dominating approach for constructing sophisticated density estimates, but they often present computational difficulties in practice. For example, undirected models, such as the Boltzmann machine
(Smolensky, 1986; Hinton et al., 2006) are able to achieve compact and efficientlycomputed latent variable representations at the cost of only providing unnormalized density estimates. Directed belief networks (Pearl, 1988; Neal, 1992; Adams et al., 2010), on the other hand, enable one to specify a priori marginals of hidden variables and are easily normalized, but require costly inference procedures. Bayesian nonparametric density estimation (e.g., (Escobar and West, 1995; Rasmussen, 2000; Adams et al., 2009)) is another flexible approach, but it often requires costly inference procedures and does not typically scale well to highdimensional data.Manifold learning provides an alternative way to implicitly characterize the density of data via a lowdimensional embedding, e.g., locallylinear embeddding (Roweis and Saul, 2000), IsoMap (Tenenbaum et al., 2000), the Gaussian process latent variable model (Lawrence, 2005), kernel PCA (Schölkopf et al., 1998), and tSNE (Van der Maaten and Hinton, 2008)
. Typically, however, these methods have emphasized visualization as the primary motivation. A notable exception is the autoencoder neural network
(Cottrell et al., 1987; Hinton and Salakhutdinov, 2006), which seeks embeddings in representation spaces that themselves can be high dimensional. Unfortunately, the autoencoder does not have a clear probabilistic interpretation (although see (Rifai et al., 2012) for a discussion).Some approaches, such as manifold Parzen windows (Vincent and Bengio, 2002)
, have attempted to tackle the combined problem of density estimation and manifold learning directly, but have faced difficulties due to the curse of dimensionality. Other approaches, such as the Bayesian GPLVM
(Titsias and Lawrence, 2010), characterize the manifold implicitly in terms of a nonlinear mapping from a representation space to the observed space. To define a coherent probabilistic model, however, it is necessary to find an invertible map between these spaces so that the density of a datum can be evaluated in the latent space without integrating over the preimage. It has proven difficult to flexibly parameterize the space of such invertible maps, however, let alone find a transformation that results in a tractable density on the representation space. Independent components analysis
(Bell and Sejnowski, 1995) and the related idea of a density network (MacKay and Gibbs, 1997)are examples of bijective models that exploit invertible linear transformations; these, however, have rather limited expressiveness. DiffeoMap
(Walder and Schölkopf, 2008) establishes a bijection close to a lowerdimensional subspace, and then projects to it. Other approaches, such as the the backconstrained GPLVM (Lawrence and Candela, 2006) attempt to approximate this bijection.In this work, we introduce the deep density model (DDM), an approach that bridges manifold discovery and density estimation. We exploit ideas from deep learning to introduce a rich and flexible class of bijective transformations of the observed space. We optimize over these transformations to obtain a map under which the implied distribution on the representation space has an approximately factorized form with known marginals. The invertibility of the map ensures that measure is not collapsed across the transformation, and as such, the determinant of the Jacobian can be computed. This leads to fullynormalized probability densities without a partition function.
The combination of rich bijective transformations with density estimation enables us to explore a variety of modeling directions for highdimensional data. As the approach is generative, we can easily sample data from a trained model without Markov chain Monte Carlo (MCMC). We present a variety of applications to the CIFAR and MNIST datasets, for proofofconcept. The deep density model also provides new possibilities for supervised learning by building Bayesian classifiers that have wellcalibrated classconditional probabilities. This additionally permits to exploit densities of unlabeled data to perform unsupervised learning, by constructing mixtures of models and training them coherently with expectation maximization.
In developing the deep density model, we also provide insight into a variety of fundamental concepts for latent variable models. Using information theoretic tools, we identify important connections between sparsity and the independence of the latent dimensions. These connections allow finding a map leads to an approximately factorized latent distribution. By understanding the distribution of the data in representation space and the transformation that gives rise to it, we can characterize the entropy of the distribution over data in the observed space. This enables making informed choices in model selection.
2 Bijections and Normalized Densities
We are interested in learning a distribution over data in a highdimensional space . We denote this (unknown) distribution as . An axiomatic assumption in machine learning is that the data contain structure, and this corresponds to distribution having most of its probability mass on a lowerdimensional, but very complicated, manifold in . Tractably parametrizing the space of such manifolds and then fitting the resulting distributions to data is a significant challenge.
Similarly, studying this distribution directly in the observed space presents both theoretical and computational difficulties. Instead, we consider how the data might be the result of a transformation from an unobserved representation space . We denote the distribution on this space as , and we assume the observed data arise from a transformation . We assume that the latent distribution has a simple factorized form:
(1) 
where the marginal factors have a simple and known univariate form.
We further make the assumption that is bijective and that is available analytically. In this case, the probability density for a point can be computed:
(2) 
In this paper, we introduce the deep density model, which discovers rich bijective transformations to map from simple latent distributions into complex observed densities. By optimizing over a large and flexible class of such bijective transformations, it is possible to discover structure in highdimensional data sets while still having a manageable, normalized density estimator.
Bijectivity is critical for ensuring that the density in Eq. (2) is normalized. This bijectivity is in contrast to many neural network approaches to latent representation where the latent space is often smaller (for an information bottleneck in, e.g., an autoencoder) or larger (for an overcomplete representation) than the observed space. When the representation space is smaller, then cannot be surjective and so multiple points in may map to the same point in , leading to overestimates of the density. In the overcomplete case, we also sacrifice bijectivity since cannot be surjective; in other words, the image of under will not span . As such, will have support beyond . That is, the latent normalization includes mass that appears in . However, taking this mass into account is a very challenging problem, whose difficulty unfortunately increases with the richness of . Thus, by using instead of the normalized density
(3) 
we underestimate the true density and are in a situation similar to that in energybased models with partition functions.
Thus, when learning a representation of data, we must make a choice between overcompleteness and bijectivity. Our model will be able to use bijectivity to produce normalized density estimates, but at the cost of requiring the representation space to be of equal dimension to the observed space. Overcomplete representation spaces are thought to allow for greater parametric flexibility, but we view the ability to produce a normalized density estimate as a significant win.
3 The Deep Density Model
Our objective is to learn a density estimate on the space from a set of training examples . We seek to produce an invertible parametric map such that the observed data under are welldescribed by a simple, factorized distribution. We denote this function as the decoder, and the representationspace data as . We construct the map via a sequence of layers, each of which applies an invertible linear transformation followed by an elementwise nonlinearity. This nonlinearity must be bijective; we use the standard sigmoid (logistic) function . The parameters for layer are square matrix
and bias vector
. Layer performs transformation . Given layers, we write the entire map as the composition(4) 
where is the composition operator and .
We also need to construct a nonbijective encoder function , which maps the observed data into the representation space. This function will be optimized to imbue the latent distribution with the properties outlined in Section 2; the necessity of will be expanded upon in Subsection 3.4.2. is composed of layers, and the th layer of has hidden units (with ) and parameters and . We denote the parameters for as and use the notational shortcuts above to write
(5) 
3.1 Regularizing the Transformations
We will train the model on data by minimizing an objective composed of several parts:
Divergence Penalty
: This determines the fit of the current encoding transformation. It forces the marginal densities of the empirical distribution of the representationspace data to match a target distribution of our choice, by penalizing divergence from it.
Invertibility Measure
: This ensures the invertibility of by penalizing poorlyconditioned transformations.
Reconstruction Loss
: This jointly penalizes the encoder and decoder to ensure that on the data.
Each of these participates in the overall objective given by:
(6) 
where are the weights of each term. We will examine each of these terms in more detail in the proceeding sections.
3.2 Divergence Penalty: Sculpting the Latent Marginals
In order to fit the deep density model to data, it is necessary to specify a measure of distance between the empirical distribution and the model distribution. We achieve this via a divergence penalty, which forces the model to distribute the mass of the data in the representation space so as to be similar to a distribution chosen a priori. This construction has several advantages: it 1) results in a known, fixed distribution on the representation space that can be used to generate fantasy data, 2) enables sparsity to be enforced as a constraint rather than a penalty weighed against the reconstruction cost, and 3) combats overfitting by explicitly requiring that some of the data have low probability under the model. This third advantage is subtle, but critical: some data must live in the tail of the distribution, in contrast to the maximum of the posterior which is simply a weighing of the MLE against the mode of the prior. See Figure 1 for a comparison of distributions produced by a various regularization techniques.
Concretely, we assume that the representation space is a unit hypercube, i.e., . As before, we assume there are data , which (for a given ) are mapped into to give . Ideally, for representation dimension , the divergence penalty would measure the difference between the marginal empirical distribution
(7) 
and a target univariate distribution that we define. In practice, we approximate by finding the best fit of a tractable parametric family and then computing the symmetrized Kullback–Liebler divergence:
(8) 
where
(9) 
Since, in this case the representation space is the unit hypercube, we choose our objective distribution to be a member of the Beta family:
(10) 
Given the data representations
, we estimate the empirical distribution of each dimension with a Beta distribution, using momentmatching to approximate its parameters:
(11)  
(12) 
where and
are the sample mean and variance, respectively. We note that a Beta distribution with parameter
produces a very sharp peak at 0, and as such allows to pursue sparsity in distribution (under the assumption that elements of small magnitude cannot be distinguished from each other).With these in hand, we get a closedform expression for our divergence penalty:
(14)  
where is the digamma function.
Furthermore, for each example, we impose an example divergence penalty, denoted as , which penalizes the distance between our objective distribution and the empirical distribution over the elements of that particular example:
Finally, our total divergence penalty is
3.2.1 Sparsity in Distribution
The traditional pursuit of sparsity entails the application of an type regularization that directly penalizes the activations in the representation space. This has several undesirable properties. First, the penalty does not differentiate between the cases where the activated units are distributed evenly among examples, and where a fixed set of units is always activated at all examples while others never are. Furthermore, it is discomforting that, in the limit of small reconstruction cost in the objective function, the regularization term is optimized if and only if all examples are identically mapped to the same point, namely zero. This forces all the activations to be small in order to have some of them vanish; it artificially forces the activation distribution to be contained in a small region around zero. Another implication of direct activity penalization is that we must search the parameter space of the regularization coefficient in order to attain our desired sparsity structure.
Instead of inducing sparsity directly, we achieve it in distribution, across examples. In practice, this arises from penalizing the KLdivergence between the empirical distribution — the actual distribution in the representation space for the given set of observations — and an appropriatelychosen target distribution , which has a peak at 0. The difference between this and the traditional approach to sparsity can be understood by considering the optimization problem as the dualization of the sparsity constraints. In the case of an type regularization, these constraints directly bound the space under the prescribed distance metric. In the case, we thus have a situation with two different points on the same contour of these constraints, one of which a more desirable of a solution than the other. On the other hand, the constrains that emerge from the divergence penalization are imposed within the probability simplex. The advantage is that a contour of these constraints corresponds to a locus of distributions with similar sparsity structures and thus similar desirability.
3.3 Invertibility: condition number penalty
Invertibility of is critical to providing a normalized density. A standard autoencoder contracts volumes around observed examples only, due to the reconstruction penalty approximating an invertible map at the observations. A true bijection, however, will guarantee conservation of volume not just at the data, but also at points we have never seen before. This will allow computation of the determinant of the Jacobian of
, and precisely specify how probability mass is reshuffled by the transformation. To ensure invertibility of the transformation, we must ensure the invertibility of each layer. As the nonlinear activation functions are fixed, invertibility is determined by the condition numbers of the
matrices in . We therefore introduce a regularization term that ensures invertibility:(15) 
Here, and
are the maximum and minimum eigenvalues of
. In this case, the curse of dimensionality becomes a blessing of dimensionality: in high dimensions, not only orthogonality is easily attained, it is difficult to escape. We find that in practice the invertibility requirement is easily satisfied, and does not constrain the algorithm at all.3.4 Reconstruction and the independence of latent dimensions
Since we now have the ability to dictate the marginal distributions in the representation space, we can shed light on the connection between entropy, sparsity and independence as a function of the transformation to the representation space. In order to attain a tractable distribution in the representation space, we must eliminate dependencies between the latent dimensions. We refer to this process as diversification, as it reduces the overlap of information learnt by distinct dimensions; see the effects of this process in Figure 2.
Diversification arises naturally when considering the effect of simultaneously increasing sparsity in the latent representation, while decreasing the entropy of the target marginals. The reconstruction penalty demands that information be preserved in the representation space. However, as the entropy decreases, the information capacity of the marginals decreases; it is necessary for the dimensions of the latent distribution to become more independent in order to reduce redundancy and continue to reconstruct successfully. Hence, we increase sparsity by limiting the capacity of the encoder, until we get to the point of minimum marginal entropy under sufficient conservation of information (which we measure by reconstruction of the observations). At this point, we expect the marginals to be approximately independent.
In the second line below, notice that we may write the observed space entropy in terms of the representation entropy and a term that accounts for the contraction of volumes under the transformation. We then write the joint entropy in terms of the marginal entropies and the mutual information information between the joint distribution and the independent marginal factorization:
(16)  
The latent dimensions are independent if and only if the mutual information is zero, and as such, we seek to minimize it. The entropy in the observed space is fixed, but we have control of the marginal entropies, and we can decrease by placing an information bottleneck on . Thus, by minimizing both these terms, we can minimize latent dependencies.
To that end, we approximate
(17) 
Note that an inherent property of the sigmoid nonlinearity is its asymptotic flatness as its output approaches zero. As such, the decoder
becomes more unable to distinguish between points as their magnitude decreases. Thus, for each representation dimension, by increasing the parameter of our objective distribution , we can shift more probability mass towards zero: this not only decreases the marginal entropies , but also increases the information bottleneck. We add that the flatness of the sigmoid still does not imply true sparsity: a deep transformation can distinguish between small but nonzero values, and as such ”undo” the sparsity effect—information is not gone, only bitshifted. As such, in order to ensure true information loss in the neighbourhood of zero, we also introduce a threshold to the encoder, that kills any representation element less than some ; namely, we use the encoder .3.4.1 Entropy characterization
A direct consequence of the above is our ability to now place an upper bound on , and, once the mutual information is minimized, to characterize it. This now allows us to make informed choice of a representation space—for example, we understand the interplay between the latent dimensionality and sparsity. Minimizing the mutual information in the way demonstrated above corresponds to selecting this space to be maximally sparse under the constraint of retaining the information encapsulated in the input examples.
We furthermore note that the entropy of model can be thought of as the ”effective number of configurations” that data drawn from its distribution can take. As such, it governs how mass is allocated between the training examples and outofsample data: thus, the diversification procedure described dictates the generalizability of the model in a very transparent way.
3.4.2 On why we need
We can now understand the motivation for introducing even with a bijective map. We need to induce points in the observed space to be close together in , both to control the information capacity as well as to determine the marginals in . On the other hand, should not expand measure in mapping from back to ; this ensures that there is an information bottleneck even under bijectivity. However, asking and to serve this double purpose is contradictory, as one function exactly undoes the contraction of measure performed by the other function. Furthermore, is helpful from a computational point of view. In the process of shaping the latent distribution, we must define our requirements as penalties on , which we proceed to backpropagate. Since asymptotic flatness of means asymptotic steepness of , performing numerical optimization on this function—let alone demanding representation points to be clustered in this asymptotically steep regime, as the divergence penalty requires—is clearly numerically unstable.
4 Training
We train the model on a GPU cluster. In cases of continuous input dimensions, we preprocess the data by whitening and normalizing it. We additionally use PCA to reduce the dimensionality of CIFAR.
In the pretraining stage, we recursively train singlelayer deep density models with stochastic gradient descent, where we take the representation space examples of one iteration as the observed space examples of the next.
In the finetuning stage, we optimize the objective in Eq. (6) by performing block coordinate descent: we iterate through the layers, and for each layer we take a step to minimize the objective as a function only of that layer’s parameters. Since different gradient magnitudes of distinct layers are not mixed here, this sidesteps the problem of having the gradient exponentially decay during backpropagation.
The optimization procedure is designed to have few free parameters: most are actively adapted in the process of optimization. The step size is chosen adaptively via an inexact Armijo’s rule line search. Exact line search on the overall objective is computationally expensive and not possible using minibatches of data. Secondly, the coefficients are adapted to maintain a specified ratio of gradient magnitudes, as the step direction is a mixture of the various penalties and is thus dictated by the relative proportions of their gradients.
We sprinkle masking noise onto the inputs to attain robust solutions, as presented in Vincent et al. (2008). We also add momentum, but implicitly: we define a window size, and sweep it across the shuffled indices of the training examples, in increments that are a fraction of the window size. As such, each example will be presented in several minibatches in a row, but with each minibatch still introducing new training examples into the objective. We found this implicit momentum technique gives the algorithm a more stable convergence.
4.1 Distribution sequencing and initialization
Enforcing the divergence penalty immediately after initialization results in a highly nonconvex and thus a very challenging optimization problem. The diversification procedure often terminates with an objective distribution with . As such, once the algorithm settles near a solution whose empirical distribution in the representation space takes this shape, it will be very difficult for it to “reshuffle” the representation data to attain an alternative solution also with the same distribution. To that end, instead of penalizing directly with the final objective distribution, we penalize through a family of distributions over iterations , where we have as an “easy” initial problem, and then have the parameter sequences converge as . The transitions between consecutive elements in these sequences are also adaptive: the objective distribution parameters are updated once the current empirical distribution is sufficiently close to the objective distribution.
In addition, we choose the sequence of hyperparameters such that
, as maintaining a constant expectation improves stability. In accordance with this, we initialize the biases as , which similarly gives rise to a consistent initial expectation. In order to initialize the weight matrices to satisfy the invertibility constraint, we select them to be random orthonormal matrices, with a scaling such that the representation distribution approximately matches in variance the initial objective distribution .5 Empirical results
In this section, we examine the properties of the model and learning algorithm on benchmark data: the MNIST digits (60,000 binary handwritten numerals)^{3}^{3}3http://yann.lecun.com/exdb/mnist/ and the CIFAR10 image data^{4}^{4}4http://www.cs.toronto.edu/~kriz/cifar.html (50,000 color images). In particular, we examine the global and local properties of the learned densities, via generation and perturbation. We emphasize that the density estimates we report here arise are fully normalized due to the tractability of the Jacobian determinant of our transformation.
Dataset  MNIST  CIFAR10 

Training examples  41.8  
Test examples (TE)  45.5  
Test examples  
rotated by  
TE with 10% of  
elements corrupted 
5.1 Density evaluation
We start by conducting a variety of tests to examine the quality of the density estimates produced by the DDM.
First, for a probability model to be useful, it should not overfit by distributing most of the mass to training examples, but also assign high probability to unseen outofsample data. Similarly, it should assign low density to points in data space that resemble real observations, which in fact are not.
In Table 1, we compute the probabilities assigned by models to their training examples, test examples, and distortions of the test examples.
As an interpretable example, we train a model on examples from MNIST’s digit 9 class, and consider the logprobability it assigns the digit 6 under rotation. Intuitively, we expect the highest density to be assigned to the upsidedown 6; see Figure 3. We also add a testset 9 to demonstrate the density calibration.
We also investigate the marginal entropy of MNIST. We train the model to have a final diversification marginals . This distribution is extremely peaked at 0 and 1, which justifies approximating the MNIST representation elements as Bernoullis by rounding off to 0 and 1 as , which corresponds to Bernoulli parameter
. We now note that, by the law of large numbers, as
, we expect the mean logprobability of the Bernoulli test examples to approach . Indeed, we compute that . Thus, we see that the marginal entropy of the model is close to our expectation of it.5.2 Generation
The approximate independence due to the process of diversification described in Subsection 3.4 combined with the invertibility property of the decoder allows us to produce real samples from the density estimation on the observed space, by sampling within the representation space. Such instantiations of the distribution provide us with visualizations of regions of high density under it; see Figure 4.
6 Discussion
In this paper, we have proposed a new way to construct normalized probability density estimates from highdimensional data. Our approach borrows ideas from deep learning, differential geometry, and information theory to ensure that the learned distributions are rich, but tractable and normalizable.
There are several interesting directions that we believe such a density model opens up. One advantage of a fullynormalized density model is that it enables Bayesian probabilistic classifiers to be constructed to provide calibrated conditional distributions over class membership. If the data are drawn from a mixture of classes, we may compute , the likelihood that an element is drawn from class by training a model with only the examples under that class. It is then possible to, e.g., estimate the most likely class:
Moreover, since we now have classconditional densities, we may reject classification of examples if no class model is confident enough to own them. Namely, we can set a threshold parameter , and reject classification of example if .
Lastly, we can not only control the density at the observations, but also the distribution of mass in the rest of the latent space. Specifically, instead of letting each class model see only its own examples, we can expose it to training examples from other classes, and demand it assigns them as little probability as possible. This not only teaches the model to recognize examples that lie on the manifold of its class, but also identify differences from other examples by mapping them away from this manifold.
We tested the above ideas on MNIST, and found that a raw mixture of Bayesian classifiers gives us an error rate of 9.5%. However, penalizing density assigned to foreign examples results in a model that achieves a much lower error rate of 1.614%. As the model provides calibrated probabilistic predictions, it is also able to assess its confidence when making classifications. Among examples in which the model is confident (approximately 95% of the test data), the Bayesian DDM classifier achieves 0.45% error. Although these are not stateoftheart rates, they show the flexibility of classifier construction in the DDM setting and how the normalized density can be leveraged across separately trained models, something not typically possible for energybased probabilistic approaches.
We can extend these ideas further to use the deep density model even if only a small fraction of the observations are labelled. That is, we can use density estimates to extract useful information from unlabeled data by leveraging our knowledge of the empirical density in the representation space. One possible approach is to run the expectationmaximization algorithm and train the DDM on weighted data as part of a mixture model.
References
 Smolensky (1986) P. Smolensky. Information processing in dynamical systems: Foundations of harmony theory. In Parallel Distributed Processing: Explorations in the Microstructure of Cognition. 1986.
 Hinton et al. (2006) G.E. Hinton, S. Osindero, and Y.W. Teh. A fast learning algorithm for deep belief nets. Neural computation, 18(7):1527–1554, 2006.
 Pearl (1988) J. Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference. Morgan Kaufmann, 1988.
 Neal (1992) Radford M. Neal. Connectionist learning in belief networks. Artificial Intelligence, 56:71–113, July 1992.
 Adams et al. (2010) Ryan P. Adams, Hanna M. Wallach, and Zoubin Ghahramani. Learning the structure of deep sparse graphical models. In Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, 2010.
 Escobar and West (1995) Michael D. Escobar and Mike West. Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association, 90(430):577–588, June 1995.

Rasmussen (2000)
Carl Edward Rasmussen.
The infinite Gaussian mixture model.
In Advances in Neural Information Processing Systems 12, pages 554–560, 2000.  Adams et al. (2009) Ryan P. Adams, Iain Murray, and David J. C. MacKay. The Gaussian process density sampler. In Advances in Neural Information Processing Systems 21, 2009.
 Roweis and Saul (2000) S.T. Roweis and L.K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323–2326, 2000.
 Tenenbaum et al. (2000) J.B. Tenenbaum, V. De Silva, and J.C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319–2323, 2000.

Lawrence (2005)
N. Lawrence.
Probabilistic nonlinear principal component analysis with Gaussian process latent variable models.
Journal of Machine Learning Research, 6:1783–1816, 2005.  Schölkopf et al. (1998) B. Schölkopf, A. Smola, and K.R. Müller. Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5):1299–1319, 1998.
 Van der Maaten and Hinton (2008) L. Van der Maaten and G. Hinton. Visualizing data using tSNE. Journal of Machine Learning Research, 9(25792605):85, 2008.
 Cottrell et al. (1987) Garrison W. Cottrell, Paul Munro, and David Zipser. Learning internal representations from grayscale images: An example of extensional programming. In Conference of the Cognitive Science Society, 1987.
 Hinton and Salakhutdinov (2006) Geoffrey Hinton and Ruslan Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504–507, 2006.
 Rifai et al. (2012) Salah Rifai, Yoshua Bengio, Yann Dauphin, and Pascal Vincent. A generative process for sampling contractive autoencoders. In International Conference on Machine Learning, 2012.
 Vincent and Bengio (2002) Pascal Vincent and Yoshua Bengio. Manifold Parzen windows. In Advances in Neural Information Processing Systems 15, pages 825–832. MIT Press, 2002.
 Titsias and Lawrence (2010) M. Titsias and N. Lawrence. Bayesian Gaussian process latent variable model. In International Conference on Artificial Intelligence and Statistics, 2010.
 Bell and Sejnowski (1995) A.J. Bell and T.J. Sejnowski. An informationmaximization approach to blind separation and blind deconvolution. Neural computation, 7(6):1129–1159, 1995.
 MacKay and Gibbs (1997) David J.C. MacKay and Mark N. Gibbs. Density networks. In Proceedings of Society for General Microbiology Edinburgh meeting, 1997.
 Walder and Schölkopf (2008) Christian Walder and Bernhard Schölkopf. Diffeomorphic dimensionality reduction. In Advances in Neural Information Processing Systems 22, 2008.
 Lawrence and Candela (2006) Neil D. Lawrence and Joaquin Quiñonero Candela. Local distance preservation in the GPLVM through back constraints. In International Conference on Machine Learning, pages 513–520, 2006.

Vincent et al. (2008)
Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and PierreAntoine Manzagol.
Extracting and composing robust features with denoising autoencoders.
In Proceedings of the 25th International Conference on Machine Learning, pages 1096–1103, 2008.