Unsupervised learning with contrastive latent variable models

11/14/2018 ∙ by Kristen Severson, et al. ∙ 0

In unsupervised learning, dimensionality reduction is an important tool for data exploration and visualization. Because these aims are typically open-ended, it can be useful to frame the problem as looking for patterns that are enriched in one dataset relative to another. These pairs of datasets occur commonly, for instance a population of interest vs. control or signal vs. signal free recordings.However, there are few methods that work on sets of data as opposed to data points or sequences. Here, we present a probabilistic model for dimensionality reduction to discover signal that is enriched in the target dataset relative to the background dataset. The data in these sets do not need to be paired or grouped beyond set membership. By using a probabilistic model where some structure is shared amongst the two datasets and some is unique to the target dataset, we are able to recover interesting structure in the latent space of the target dataset. The method also has the advantages of a probabilistic model, namely that it allows for the incorporation of prior information, handles missing data, and can be generalized to different distributional assumptions. We describe several possible variations of the model and demonstrate the application of the technique to de-noising, feature selection, and subgroup discovery settings.



There are no comments yet.


page 8

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In unsupervised learning, the goal is often to learn what is unique or interesting about a dataset. Given the subjective nature of this question, it can be useful to frame the problem in the context of what signal is enriched in one dataset, referred to as the target, relative to a second dataset, referred to as the background. An example of this is an exploration of a heterogeneous disease population, such as patients with Parkinson’s disease. The interesting sources of variation are those that are unique to the disease population. However, it is likely that some sources of variation are unrelated to the disease state, for instance variation due to aging. This is difficult to assess without a baseline population, therefore, it is useful to contrast the disease population with a population of healthy controls. Such contrastive analysis can discover nuisance variation that is common amongst the two populations and is uninteresting for the problem while highlighting variation unique to the disease population enabling downstream applications such as subgroup discovery.

Despite this natural setting for unsupervised learning, most techniques address individual data points, sequences, or paired data points. Few techniques generalize to the contrastive scenario where we have sets of data but no obvious correspondence between their members. Yet, there are many cases where datasets that can be used in a comparative setting arise naturally: control vs. study populations, pre- and post-intervention groups, and signal vs. signal free groups [1]. Each of these settings has possible nuisance variation, for example, population level variation, effects unrelated to intervention, and sensor noise variation.

The recently published contrastive principal component approach (cPCA) [1]

is one example of a technique that can be used for sets of data. cPCA builds on principal component analysis (PCA)


, a dimensionality reduction technique which projects data into a lower dimensional space while minimizing the squared loss. PCA and other dimensionality reduction techniques are popular because they allow high-dimensional data to be visualized while removing noise. cPCA seeks to find a projection to a lower dimensional space that discovers variation that is enriched in one dataset as compared to another by applying PCA to the empirical covariance matrix


where are the observations of interest, are the comparison data, and is a tuning parameter. The choice of

is a trade-off between maximizing the retained variance of the target set and minimizing the retained variance of the background set.

In this work, we develop probabilistic latent variable models applicable to the setting where contrastive analysis is desired. These models are based on the insight that it is possible to emphasize latent structures of interest while suppressing spurious, uninteresting variance in the data through carefully designed statistical models. Such models have several key advantages over deterministic approaches: it is straight forward to incorporate prior domain knowledge, missing and noisy data can naturally be modeled through appropriate noise distributions, model and feature selection can be performed through sparsity promoting prior distributions, and the model can more easily be incorporated into larger probabilistic systems in a principled manner. Through this paper, we advance the state-of-the-art in several ways. First, we develop latent variable models capable of contrastive analysis. We then demonstrate the generality of our framework by demonstrating how robust and sparse contrastive variants can be developed, learned and how automatic model selection can be performed. We also develop contrastive variants of the variational autoencoder, a deep generative model, and demonstrate its utility in modeling the density of noisy data. Finally, we vet our proposed models through extensive experiments on real world scientific data to demonstrate the utility of the proposed framework.

2 Contrastive Latent Variable Models

To achieve the aim of discovering patterns that are enriched in one dataset relative to another, we propose a latent variable model where some structure is shared across the two datasets and some structure is unique to the target dataset. Given a target dataset and a background dataset , the model is specified


where are the observed data, and are the latent variables, and are the corresponding factor loadings, are the dataset-specific means and are the noise. In general, we do not expect the number of samples in the two datasets to be the same, i.e. . Furthermore, there is no special relationship between the samples and in equation 2. The primary variables of interest are , which are the lower dimensional representation that is unique to the target dataset.

2.1 Gaussian likelihood and priors

To provide intuition into why eqn. 2

meets our goal of capturing patterns enriched in the target with respect to the background, we consider the case where the noise follows isotropic Gaussian distributions,

and and the latent variables are modeled using standard Gaussian distributions



is a multivariate normal distribution parameterized by mean

and covariance and denotes a identity matrix. The resulting marginal distributions for the observed data are


The covariance structure for the target data is additive and contains a term ( that is shared with the background data and a term that is unique to the target data (). This constructions allows the factor loading W to model the structure unique to the target. The model closely mirrors probabilistic PCA (PPCA) [3, 4] and is exactly PPCA applied to the combined datasets when the target factor loading dimensionality is zero. Similarly, this model is exactly PPCA applied to only the target dataset when the shared factor loading dimensionality

is zero. Expectation-maximization (EM)

[5] can be used to solve for the model parameters. Because EM requires conjugacy, most model formulations will not be solved this way. However, we present a summary of the EM steps to provide an intuition about the model. To provide interpretable equations in the below description, we consider the case where the factor loading matrices W and S are orthogonal.

The model parameters are S, W, , , and the latent variables are , , . The lower bound of the likelihood is


The M-step maximizes the lower bound of the likelihood with respect to the parameters. The update step for the shared factor loading is


where B is the sample covariance of the background data, T is the sample covariance of the target data, , and . The update step for the target factor loading is


Details on the derivation can be found in the supplemental information. It is useful to recall that the orthogonal projection onto the range space of a matrix A is given by and the orthogonal projection onto the nullspace of A is given by . In eqn. 7, can be expanded using the definition of M to . Similarly, in eqn. 6, can be expanded using the definition of Q to . When is small, these equations are similar to the projection onto the nullspace of S and W, respectively. This matches our intuition as to how these factor loading matrices are updated: in a sense, the part of the target data captured by the target factor loading space is projected away before updating the shared factor loading space, and vice versa. This behavior is similar to cPCA. In eqn. 1, as goes to infinity, directions not in the null space of the background data covariance are given an infinite penalty. When this is the case, cPCA projects the target data onto the null space of the background data and then performs PCA [1].

The update steps can also be compared to the PPCA updates. For the factor loading matrix W, the update step is:


which is the same as eqn. 7, except for the projection term.

2.2 Beyond Gaussian Models

The assumptions of Gaussianity are not necessary for recovering latent structure that enriches desired patterns in the target dataset. We can more generally express the proposed model as:


where and . The primary modeling decisions are to choose the appropriate likelihoods and priors on the loading matrices. The particular choices are governed by the application and domain specific knowledge.

However, this flexibility comes at a price: the posterior distributions are no longer guaranteed to be tractable. Consequently, the EM algorithm sketched in the previous section is no longer available and instead, we use variational inference [6]. In summary, the intractable posteriors are approximated with tractable surrogates and divergence is minimized with respect to the variational parameters . This is equivalent to maximizing the lower bound of the marginal likelihood,


where implies the parameters in except the parameters denoted in the set. Depending on the choice of and the expectations required for computing may themselves be intractable. We use recently proposed black box techniques [7, 8, 9, 10] to sidestep this additional complication. In particular, we approximate the intractable expectations in

with unbiased Monte-Carlo estimates,

. Because the latent variables of interest are continuous, we are able to use reparameterization gradients [8, 9] to differentiate through the sampling process and obtain low variance estimates of , . Using the noisy but unbiased gradients, optimization can proceed using a stochastic gradient ascent variant, e.g. ADAM [11]. In our experiments we use Edward [12], a library for probabilistic modeling, to implement these inference strategies for the proposed models. We sketch the pseudocode for variational learning in Algorithm 1.

1:  Input Model , variational approximations
2:  Output: Optimized and variational parameters
3:  Initialize and .
4:  repeat

     Use reparameterization trick to compute unbiased estimates of the gradients of the objective in Eqn. 

6:     Update ,
7:  until convergence
Algorithm 1 Pseudocode

Finally, we note that the black box inference framework does not restrict us to point estimates of . As we will illustrate in the next section, it is possible to infer variational distributions over by specifying an appropriate approximation .

2.3 cLVM Variants

We refer to the base structure of the model as provided in eqn. 9 as a contrastive latent variable model, cLVM. As previously noted, different choices for the distributions in eqn. 9 can be made to address the specific challenges of the application. Several models are introduced here and are summarized in Table 1.

Model name Prior Likelihood Variational Approximation
cLVM Gaussian
Sparse cLVM
cLVM with
model selection
Robust cLVM Student’s t
Gaussian parameterized

by neural network

Table 1: Summary of the model variants. For all of the models in the table, the latent variables are modeled as standard Gaussians and the variational distributions are also Gaussian, unless otherwise noted. The model choice depends on the application. The various models are not mutually exclusive and may also be combined.

Sparse cLVM

One application-specific problem is feature selection. In unsupervised learning, there is often a secondary goal of learning a subset of measurements that are of interest which is motivated by improved interpretability. This is especially important when the observed data is very high-dimensional. For instance, many biological assays result in datasets that have tens of thousands of measurements such as SNP and RNA-Seq data. During data exploration, discovering a subset of these measurements that is important to the target population can help guide further analysis. To learn a latent representation that is only a function of a subset of the observed dimensions, certain rows of the target factor loading, W, must be zero. The observed data corresponding to the zero rows in W then have no contribution to the latent representation t. Because there is no restriction on S, a sparsity requirement for W does not imply that the corresponding observation is zero.

One way to achieve this behavior is by using a regularization penalty on the model parameters. The penalty is added to the objective function to incite certain behavior. Regularization penalties can be related to priors by noting that , where is the penalty function. For feature selection, a group sparsity penalty [13] could be used. The rows of are penalized:


where is the row of W. This functional form is known to lead to sparsity at the group level, i.e. all members of a group are zero or non-zero. For increasing values of , the target factor loading matrix has a larger number of zero-valued rows.

Sparsity inducing priors such as the automatic relevance determination (ARD) [14, 15, 16] or global-local shrinkage priors such as the horseshoe [17, 18] can also be easily incorporated into the frameworkU̇sing the horseshoe prior as an example, the row of W is modeled,



is the half-Cauchy distribution with density

for . The horseshoe prior is useful for subset selection because it has heavy tails and an infinite spike at zero. Further discussion can be found in the supplemental information. For both the prior and regularization formulations, groups of rows in W could also be used instead of single rows if such a grouping exists.

cLVM with Automatic Model Selection

The ARD prior is more typically applied to the columns of a factor loading matrix. This use allows for automatic selection of the dimension of the matrix. This could also be done in the cLVM model. Although both latent spaces can have any dimension less than , which must be selected, we generally recommend setting the target dimension to two for visualization purposes. To select the dimension of the shared space, the percent variance explained can be analyzed or a prior, such as the ARD prior can be used. The columns of S are modeled


The ARD prior has been shown to be effective at model selection for PPCA models [19].

Robust cLVM

Another application-specific goal may be to systematically handle outliers in the dataset. Similar to PPCA, the cLVM model is sensitive to outliers and can produce poor results if outliers are not addressed. It may be possible to remove outliers from the dataset, however this is typically a manual process that requires domain expertise and an understanding of the process that generated the data. A more general approach to handling outliers uses a heavy-tailed distribution to describe the data. One approach for constructing heavy tailed distributions is through scale mixtures of Gaussians 

[20]. Consider,


The resulting marginal distribution of the observed data is


where St indicates a Student’s t-distribution [21]

. The larger probability mass in the tails of the Student’s t-distribution, as compared to the normal distribution, allows the model to be more robust to outliers.

2.4 Beyond Linear Models

Contrastive Variational Autoencoders

Thus far we have only considered models that linearly map latent variables and to the observed space. The linearity constraint can be relaxed, and doing so leads to powerful generative models capable of accounting for nuisance variance.


where , , and ,

represent non-linear transformations parameterized by neural networks. The latent variables are modeled using standard Gaussian distributions, as before. Observe that similar to the linear case (eqn. 

2) the target and background data share the projection while the target retains a private projection . This construction forces to model commonalities between the target and background data while allowing to capture structure unique to the target.

This model can be learned by maximizing the lower bound to the marginal likelihood , , analogously to eqn. 10. However, a large amount of data is typically required to learn such a non-linear model well. Moreover, since the number of latent variables proliferate with increasing data, it is computationally more efficient to amortize the cost of inferring the latent variables through inference networks shared between the data instances. In particular, we parametrize the variational posteriors and , where and are inference network parameters. Unlike eqn. 10 where the variational parameters grow with the number of data instances, the variational parameters and do not. is shared amongst the target instances while is shared between the background examples. This is an example of amortized variational inference [22, 23]. Finally, learning proceeds by maximizing the evidence lower bound,


with respect to and , . The KL terms are available to us in closed form, however the expectation terms are intractable and we again resort to Monte Carlo approximations and re-parameterized gradients to enable stochastic gradient ascent. We refer to this combination of the non-linear model and the amortized variational inference scheme as the contrastive variational auto encoder (cVAE).

3 Related Work

There are many techniques for dimensionality reduction, e.g. [2, 24, 25]. This review focuses on dimensionality techniques that use sets of data and/or address issues related to nuisance variation. Canonical correlation analysis (CCA) [26] and its probabilistic variant (PCCA) [27] use two (or more) sets of data, however requires that samples are paired views (or higher dimensional sets of views) of the same sample. For instance perhaps several tests are run on a single patient and therefore the tests are linked via the patient identity. In CCA, the number of samples in the sets must be equal, , however the dimensionality of each sample does not need to be the same. Damianou, Lawrence and Ek [28] proposed a nonlinear extension of PCCA where the mappings are sampled from a Gaussian process. The resulting model is a multi-view extension of GP-LVM [29], but still requires linking the samples across datasets.

In this work, we propose addressing nuisance variation in the dataset by introducing a structure to the latent representation. Schulam and Saria [30] investigate a similar idea with respect to sharing representations across different parts of the full data. In their work, a hierarchical model for disease trajectory is proposed where some of the model coefficients are shared across subsets of the data, e.g. total population and individual. This idea has also been proposed for the unsupervised analysis of time series data [31, 32]. Data samples are assumed to have a latent representation that can be partitioned into static and dynamic contributions. None of these works have considered a contrastive setting. There has also been work in addressing explicit sources of nuisance variation. Louizos et al. [33] explores a setting where certain variables within the dataset are a priori identified as nuisance and the remaining variables contribute to the latent representation. The observed data is modeled where s are the observed nuisance variables.

4 Experiments

Contrastive latent variable models have applications in subgroup discovery, feature selection, and de-noising, each of which is demonstrated here leveraging different modeling choices. We use examples from [1] to highlight the similarities and differences between the two approaches. The results of cLVM as applied to synthetic datasets can be found in the supplemental information.

Figure 1: cLVM is robust to missing data. Density plots of the subgroups revealed in the target latent representation of the mice protein expression data. Red and blue points are the control and trisomic mice samples, respectively. The rows use cPCA and robust cLVM to learn the latent representation, respectively. Each column uses a different level of missing data, starting with the leftmost column containing the natural level of missing data. PCA is unable to perform subgroup discovery (see supplemental information) and robust cLVM is better able to perform subgroup discovery in the presence of missing data.
Figure 2: cLVM variants allow for model and feature selection. (a) Subgroups revealed in the target latent representation for the RNA-Seq dataset using the model selection cLVM variant. (b) The percent variance explained by the ordered columns of the shared factor loading for LVM with and without ARD (model selection). The ARD model has over 100 fewer non-zero columns in the shared factor loading. (c) Subgroups revealed in the target latent representation for the mHealth dataset using sparse cLVM. (d) The norms of the rows of the target factor loading for sparse LVM where the different colors correspond to different sensor types. The six dimensions with zero-valued norms correspond to magnetometer readings.
Figure 3: cVAE recovers meaningful structure from noisy data. (a) Samples of the target noisy images of digits on grass and background grass images. (b) Generative samples of the de-noised target (top row) and background (bottom row) which are enabled by the cVAE structure. Note there is no correspondence between the samples in (a) and (b). (c) The 2D cVAE projection and a 2D tSNE projection of a VAE with 10 dimensional space. The colors represent different digits.

4.1 Subgroup Discovery for Incomplete Data

To demonstrate the use of cLVM for subgroup discovery, we use a dataset of mice protein expression levels [34]. The target dataset has 270 samples of two unknown classes of mice: trisomic (Down Syndrome model) and control. The background dataset has 135 known control samples. Each sample has 77 measurements. The dataset contains missing values at a level of approximately 1.6% due to technical artifacts and sampling that cannot be repeated. One of the advantages of the probabilistic approach is that it naturally handles missing data. Depending why the data is missing, missing data can either be ignored, marginalized, or explicitly modeled. For the mice protein dataset, we marginalize over the missing values by treating the missing data as latent variables and adding a corresponding normal variational approximation. Increasing levels of missing data were tested by artificially removing data from the target dataset. The robust variation of the model is applied to account for other possible data issues. The target and shared dimensionalities are both set to two. Fig. 7

shows the latent representation using cPCA and robust cLVM for the naturally occurring missing level, 25%, 50%, and 75% missing data. cPCA does not have natural handling for missing data therefore mean imputation was used to first fill-in. PCA is unable to recover the structure in the dataset (see supplemental information for results). Both cPCA and robust cLVM find the subsets, however, the proposed method is better able to discover the subgroups as the amount of missing data increases.

4.2 Subgroup Discovery for High Dimensional Data

To highlight the use of cLVM for subgroup discovery in high-dimensional data, we use a dataset of single cell RNA-Seq measurements [35]. The target dataset consists of expression levels from 7,898 samples of bone marrow mononuclear cells before and after stem cell transplant from a leukemia patient. The background contains expression levels from 1,985 samples from a healthy individual. Pre-processing of the data reduces the dimensionality from 32,738 to 500 [35, 1]. Given the size of the data to explore, it is useful in this setting to use an ARD prior to automatically select the dimensionality of the shared latent space. The target latent space is set to two and an prior is used for the columns of the shared factor loading. Fig. 8a shows the resulting latent representation, which is able to discover the subgroups, whereas PCA is not (see supplemental information). Fig. 8b compares the percent of variance explained in the ranked columns as compared to the cLVM model without model selection. The model with ARD uses over 100 fewer columns in the shared factor loading matrix and avoids an analysis to manually select the dimension.

4.3 Automatic Feature Selection using Sparse cLVM

The third example uses a dataset, referred to as mHealth, that contains 23 measurements of body motion and vital signs from four types of signals [36, 37]. The participants in the study complete a variety of activities. The target data is composed of the unknown classes of cycling and squatting and the background data is composed of the subjects lying still. In this application, we demonstrate feature selection by learning a latent representation that both separates the two activities and uses only a subset of the signals. A group sparsity penalty is used, as described in the methodology, on the target factor loading. The target dimension is two, the shared dimension is twenty, and is 400. is selected by varying its value and inspecting the latent representation. The latent representation using regularization is shown in Fig. 8c. The two classes are clearly separated. Fig. 8d shows the row-wise norms of the target factor loading. The last six dimensions, corresponding to the magnetometer readings, are all zero which indicates that the magnetometer measurements are not important for differentiating the two classes and can be excluded from further analysis.

4.4 De-noised Generative Modeling using cVAE

Finally, to demonstrate the utility of cVAE, we consider a dataset of corrupted images (see Fig. 3a). This dataset was created by overlaying a randomly selected set of MNIST [38]

digits on randomly selected images of the grass category from Imagenet 

[39]. The background is grass images. We train a cVAE with a two-dimensional target latent space and an eight-dimensional shared space. We use fully connected encoder and decoder networks with two hidden layers with and hidden units employing rectified-linear non-linearities. For the cVAE, both the target and shared decoders and use identical architectures. We compare against a standard variational autoencoder with an identical architecture and employ a latent dimensionality of ten, to match the combined dimensionality of the shared and target spaces of the contrastive variant. Fig. 3c presents the results of this experiment. The latent projections for the cVAE cluster according to the digit labels. VAE on the other hand confounds the digits with the background and fails to recover meaningful latent projections. Moreover, cVAE allows us to selectively generate samples from the target or the background space, Fig. 3b. The samples from the target space capture the digits, while the background samples capture the coarse texture seen in the grass images. Additional comparisons with a VAE using a two dimensional latent space is available in the supplemental.

5 Conclusions

Dimensionality reduction methods are important tools for unsupervised data exploration and visualization. We propose a probabilistic model for improved visualization when the goal is to learn structure in one dataset that is enriched as compared to another. The latent variable model’s core characteristic is that it shares some structure across the two datasets and maintains unique structure for the dataset of interest. The resulting cLVM model is demonstrated using robust, sparse, and nonlinear variations. The method is well-suited to scenarios where there is a control dataset, which is common in scientific and industrial applications.


  • [1] A. Abid, M. J. Zhang, V. K. Bagaria, and J. Zou, “Exploring patterns enriched in a dataset with contrastive principal component analysis,” Nature Communications, vol. 9, p. 2134, 2018.
  • [2] H. Hotelling, “Analysis of a complex of statistical variables into principal components,” Journal of Educational Psychology, vol. 24, p. 417, 1933.
  • [3] M. E. Tipping and C. M. Bishop, “Probabilistic principal component analysis,” Journal of the Royal Statistical Society, Series B, vol. 61, pp. 611–622, 1999.
  • [4] S. T. Roweis, “EM algorithms for PCA and SPCA,” in NIPS, 1998, pp. 626–632.
  • [5] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” Journal of the Royal Statistics Society. Series B (Methodological), vol. 39, pp. 1–38, 1977.
  • [6] M. J. Wainwright and M. I. Jordan, “Graphical models, exponential families, and variational inference,”

    Foundations and Trends® in Machine Learning

    , vol. 1, no. 1–2, pp. 1–305, 2008.
  • [7] R. Ranganath, S. Gerrish, and D. M. Blei, “Black box variational inference,” in AISTATS, 2014, pp. 814–822.
  • [8] D. P. Kingma and M. Welling, “Stochastic gradient vb and the variational auto-encoder,” in ICLR, 2014.
  • [9] D. J. Rezende, S. Mohamed, and D. Wierstra, “Stochastic backpropogration and approximate inference in deep generative models,” in ICML.   PMLR, 2014, pp. 1278–1286.
  • [10] M. Titsias and M. Lázaro-Gredilla, “Doubly stochastic variational Bayes for non-conjugate inference,” in ICML.   PMLR, 2014, pp. 1971–1979.
  • [11] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  • [12] D. Tran, A. Kucukelbir, A. B. Dieng, M. Rudolph, D. Liang, and D. M. Blei, “Edward: A library for probabilistic modeling, inference, and criticism,” arXiv preprint arXiv:1610.09787, 2016.
  • [13] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” Journal of the Royal Statistical Society, Series B, vol. 68, pp. 49–67, 2007.
  • [14] C. M. Bishop, “Variational principal components,” in ICANN.   IEE, 1999, pp. 509–514.
  • [15] S. Virtanen, J. Jia, A. Klami, and T. Darrell, “Factorized multi-modal topic model,” in

    Proceedings of the Conference on Uncertainty in Artificial Intelligence

    .   ACM, 2011, pp. 843–851.
  • [16] A. Klami, S. Virtanen, and S. Kaski, “Bayesian canonical correlation analysis,” Journal of Machine Learning Research, vol. 14, pp. 965–1003, 2013.
  • [17] C. M. Carvalho, N. G. Polson, and J. G. Scott, “Handing sparsity via the horseshoe,” in AISTATS.   JMLR, 2009, pp. 73–80.
  • [18] ——, “The horseshoe estimator for sparse signals,” Biometrika, vol. 97, pp. 465–480, 2010.
  • [19] C. M. Bishop, “Bayesian PCA,” in NIPS, 1999, pp. 382–388.
  • [20] M. West, “On scale mixtures of normal distributions,” Biometrika, vol. 74, pp. 646–648, 1987.
  • [21] C. Archambeau, N. Delannay, and M. Verleysen, “Robust probabilistic projections,” in ICML.   ACM, 2006, pp. 33–40.
  • [22] P. Dayan, G. E. Hinton, R. M. Neal, and R. S. Zemel, “The helmholtz machine,” Neural computation, vol. 7, no. 5, pp. 889–904, 1995.
  • [23] S. Gershman and N. Goodman, “Amortized inference in probabilistic reasoning,” in Proceedings of the Annual Meeting of the Cognitive Science Society, vol. 36, no. 36, 2014.
  • [24] L. van der Maaten and G. Hinton, “Visualizing data using t-SNE,” Journal of Machine Learning Research, vol. 9, pp. 2579–2605, 2008.
  • [25] M. A. A. Cox and T. F. Cox, “Multidimensional scaling,” in

    Handbook of Data Visualizations

    .   Springer, 2008, pp. 315–347.
  • [26] H. Hotelling, “Relations between two sets of variables,” Biometrika, vol. 28, pp. 321–377, 1936.
  • [27] F. R. Bach and M. I. Jordan, “A probabilistic interpretation of canonical correlation analysis,” University of California, Berkeley, Tech. Rep., 2005.
  • [28] A. Damianou, N. D. Lawrence, and C. H. Ek, “Multi-view learning as a nonparametric nonlinear inter-battery factor analysis,” arXiv preprint arXiv:1604.04939, 2016.
  • [29] N. Lawrence, “Probabilistic non-linear principal component analysis with Gaussian process latent variable models,” Journal of Machine Learning Research, vol. 6, pp. 1783–1816, 2005.
  • [30] P. Schulam and S. Saria, “A framework for individualizing of disease trajectories by exploiting multi-resolution structure,” in NIPS, 2015, pp. 748–756.
  • [31] W.-N. Hsu, Y. Zhang, and J. Glass, “Unsupervised learning of disentangled and interpretable representations from sequential data,” in NIPS, 2017, pp. 1878–1889.
  • [32] Y. Li and S. Mandt, “Disentangled sequential autoencoder,” in ICML.   PMLR, 2018, pp. 5656–5665.
  • [33] C. Louizos, K. Swersky, Y. Li, M. Welling, and R. Zemel, “The variational fair autoencoder,” in ICLR, 2016.
  • [34] C. Higuera, K. J. Gardiner, and K. J. Cios, “Self-organizing feature maps identify proteins critical to learning in a mouse model of down syndrome,” PLOS ONE, vol. 10, p. e0129126, 2015.
  • [35] G. X. Y. Zheng, J. M. Terry, P. Belgrader, and P. e. Ryvkin, “Massively parallel digital transcriptional profiling of single cells,” Nature Communications, vol. 8, p. 14049, 2017.
  • [36] O. Banos, R. Garcia, J. A. Holgado, M. Damas, H. Pomares, I. Rojas, A. Saez, and C. Villalonga, “mhealthdroid: A novel framework for agile development of mobile health applications,” in 6th International Work-conference on Ambient Assisted Living and Daily Activities.   Springer, 2014, pp. 91–98.
  • [37] O. Banos, C. Villaonga, G. Rafael, A. Saez, M. Damas, J. A. Holgado-Terriza, S. Lee, H. Pomares, and I. Rojas, “Design, implementation, and validation of a novel open framework for agile development of mobile health applications,” BioMedical Engineering Online, vol. 14, pp. 1–20, 2015.
  • [38] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” in Proceedings of the IEEE, 1998, pp. 2278–2324.
  • [39] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and F.-F. Li, “Imagenet large scale visual recognition challenge,”

    International Journal of Computer Vision

    , vol. 115, pp. 211–252, 2015.
  • [40] M. P. Wand, J. T. Ormerad, S. A. Padoan, and Frühwirth, “Mean field variational bayes for elaborate distributions,” Bayesian Analysis, vol. 6, pp. 847–900, 2011.

Appendix A Expectation-maximization

For the EM algorithm, we consider the conditional joint distribution of



To provide intuition about the model, we consider the simplifying case where the factor loading matrices are orthogonal, therefore . The distributions can then be written:


and the corresponding expectations are:


where and


where . Because of the orthogonality constraint, the expectation of the outer product of and is the outer products of the expectations. The update equations are derived by taking the derivative of the lower bound of the likelihood with respect to each of the parameters. The update equation for the shared factor loading is


Following the same ideas for W


These equations can be simplified by plugging in the definitions of the expectations, as is shown in the main text.

Appendix B Model variants

Several different priors can be used to obtain a sparse target loading matrix. Here we present a discussion of the pros and cons of these modeling choices. As mentioned in the manuscript, the horseshoe prior is useful for subset selection because it has heavy tails and an infinite spike at zero. The infinite spike at zero provides strong shrinkage while the heavy tails allow coefficients to “escape” penalization as compared to the ARD prior. In the context of sparse cLVM, the horseshoe prior is used


In practice to use the horseshoe prior, we suggest reparameterizing the half-Cauchy distributions as inverse gamma distributions


where IG is the inverse gamma distribution [40]. We find that standard exponential family distributions are better able to approximate the inverse gamma distributions leading to improvements during the learning phase. The variational approximation used for the inverse gamma distribution is log-normal. In the experiments section below, the sparse cLVM using the horseshoe prior is applied to the mHealth dataset.

Finally, we note that the horseshoe prior will shrink the weight values but will not set them equal to zero, therefore a thresholding rule is required to prune the weights entirely. We suggest a pruning rule that considers the posterior distribution of the scales . For example, prune rows of W, when the probability that the scale is less than a sufficiently small value is greater than , i.e. .

Appendix C Experiments

In this section, we present additional results on synthetic datasets and further comparisons for the real datasets.

c.1 Synthetic data

Figure 4: LVM applied to the synthetic test data. The title of the plot is the shared dimension of the latent space for the target and background samples. The dimension of the target latent space is 2. When the shared dimension is zero, the method is equivalent to applying PPCA to the target dataset.
Figure 5: PCA, cPCA, and cLVM with automatic model selection applied to the generative synthetic test data. PCA is unable to discover the subtypes in the data whereas cPCA and cLVM
Figure 6: The target latent representation without outliers is plotted using empty circles. The left plot shows the target latent representation learned in the presence of outliers using cLVM (square points). The left plot shows the target latent representation learned in the presence of outliers using robust cLVM (square points). The robust cLVM model variant is better able to recover the target latent representation without outliers.
Subgroup A
Subgroup B
Subgroup C
Subgroup D
Features 1-10
Features 11-20
Features 21-30
Table 2: A summary of the distributional assumptions used to generate the synthetic data used in analysis.

cLVM is demonstrated on two synthetic datasets. The first is based on the synthetic dataset proposed in [1]. 30-dimensional data is simulated using the distributional assumptions summarized in Table 2. The target and background datasets each have 400 samples. The target latent representation for different choices of the shared latent dimensionality are shown in Fig. 4. The second synthetic dataset is generated using the distributional assumptions using eqn 3 where the target latent variables have cluster specific means. The target latent representation for different dimensionality reduction techniques are shown in Fig. 5

. cLVM with automatic model selection is able to correctly discover the dimensionality of the shared factor loading. A third example is shown by adding outliers to the second synthetic dataset. Outliers are drawn from a uniform distribution [-20, 20] and 20 outliers each are added to the target and background datasets. cLVM and robust cLVM are applied to the dataset. The target latent representations are shown in Fig. 


c.2 Mice Protein Expression Data

The latent representation of the mice protein expression data is also compared to PCA (see Fig. 7, first row). For all levels of missing data, PCA is unable to recover the latent subgroups.

Figure 7: Density plots of the subgroups revealed in the target latent representation of the mice protein expression data. Red and blue points are the control and trisomic mouse samples, respectively. The rows use PCA, cPCA, and robust cLVM to learn the latent representation, respectively. Each column uses a different level of missing data, starting with the leftmost column containing the natural level of missing data. PCA is unable to perform subgroup discovery and LVM is better able to perform subgroup discovery in the presence of missing data.

c.3 Single Cell RNA-Seq Data

The latent representation of the single cell RNA-Seq data is compared to PCA and cPCA in Fig. 8

. PCA is unable to recover the latent subgroups. cPCA and cLVM recover similar latent representation, however cLVM does not require a priori selection of tuning parameters whereas cPCA uses a heuristic to select the value of

, which controls the tradeoff between maximizing the variance in the target dataset and minimizing the variance in the background dataset.

Figure 8: The latent representation of the single cell RNA-Seq data. Red and blue points are the post- and pre-treatment samples, respectively. From left to right, the plots use PCA, cPCA, and cLVM with model selection to learn the latent representation, respectively. PCA is unable to perform subgroup discovery. cPCA is able to perform subgroup discovery and uses a heuristic to selection , the tuning parameter for the empirical covariance matrix (eqn. 1). cLVM with model selection finds a similar latent representation to cPCA without the need to select the shared dimensionality a priori by using automatic relevance detection.

c.4 mHealth Sensor Data

The target latent representation for the mHealth data with and without regularization is presented in Fig. 9. In both cases, there are two clear subgroups in the target latent representation, which correspond to squatting and cycling. For comparison to a probabilistic result, we also present the target latent representation using the horseshoe prior in Fig. 10. By using the horseshoe prior, we also have information about the distribution of the factor loading matrix. The distributions of the scales which govern the variance of the factor loading matrix are shown in Fig. 11

Figure 9: The target latent representations for the mHealth dataset with and without regularization. In both cases, there are clearly two subgroups, which correspond to the squatting and cycling activities. (c) The row norms of the target factor loading matrix. In the case with regularization, the measurements corresponding the to the magnetometer (last six measurements) are zero.
Figure 10: (a) The target latent representation for the mHealth dataset using a horseshoe prior for the target factor loading. (b) The row norms of the target factor loading. The sensor measurements are ordered the same as in Fig. 9.
Figure 11: The distribution of the scales () for each of the rows of the target factor loading matrix W. The title of the plot is the type of sensor: accelerometer (Acc), electrocardiogram (EGC), gyroscope (Gyro), and magnetometer (Mag). X, Y, and Z indicate the axis of the measurement. Note that the y-axis is not shared amongst the plots and the magnetometer measurements are highly peaked near zero. Increasing values of a pruning threshold will cause an increasing number of the magnetometer sensors coefficients to go to zero. A reasonably chosen threshold may also result in pruning the Acc Chest X-axis sensor.

c.5 MNIST Digits on Grass

Latent representations recovered by contrastive and regular VAEs on the MNIST on Grass dataset are shown in Fig. 12.

Figure 12: Latent projections from cVAE, VAE with 2D latent space, 2D tSNE projections of VAE with 10 dimensional latent space.