<figure><img src="image_1.png" /></figure>

HAL
open science

# From Unsupervised to Guided Clustering: A Variational Implementation

Violaine Courrier, Christophe Biernacki

► To cite this version:

Violaine Courrier, Christophe Biernacki. From Unsupervised to Guided Clustering: A Variational Implementation. 2026. hal-05489363

HAL Id: hal-05489363

https://inria.hal.science/hal-05489363v1

Preprint submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# From Unsupervised to Guided Clustering:
A Variational Implementation

Violaine Courrier*

Withings, Inria, Université de Lille, CNRS
and

Christophe Biernacki

Inria, Université de Lille, CNRS

## Abstract

Clustering is viewed as an unsupervised technique, but in practice it requires guidance to uncover meaningful structures. We formalize this with guided clustering, a paradigm that uses a guiding variable to steer the discovery process, and introduce the Guided Clustering Variational Autoencoder (GCVAE) as its deep generative realization. GCVAE learns a latent space structured as a Gaussian Mixture Model by optimizing a variational objective that forces the representation to be maximally informative about the guiding variable. This framework allows the resulting clustering to be reoriented by changing the guiding variable, yielding clusters that are meaningful for the specified context. Experiments on public (MNIST-SVHN) and proprietary connected health devices data demonstrate GCVAE's ability to discover coherent and task-relevant clusters in complex, high-dimensional settings.

*Keywords:* clustering, variational inference

<footnote>*Corresponding author: violaine.courrier@inria.fr</footnote>
# 1 Introduction

Clustering is never a purely unsupervised task. In practice, analysts implicitly guide the discovery process through choices of parameters, distance metrics, or feature selection to find “interesting” results. This hands-on necessity stems from a more fundamental challenge: the inherent ambiguity of the clustering task itself. A single dataset can contain multiple, equally valid partitions, and the most meaningful one depends entirely on the analytical goal. Without an explicit objective, the search for a relevant partition relies on an often informal trial-and-error guidance. We argue that this external knowledge should not be an afterthought, but a formal component of the model definition.

We propose to formalize this via the so-called *guided clustering*. In this new paradigm, the guiding variable is not the end goal but rather the lens through which we discover the most relevant and coherent partition of the data. We propose a specific implementation of this novel guided clustering approach using a deep generative model that learns to compress the input into a representation that is simultaneously organized into a discrete mixture of clusters and optimized to be maximally predictive of the guiding variable. The Guided Clustering Variational Autoencoder (GCVAE) model directly addresses the limitations of prior works by learning non-linear, structurally meaningful clusters within a compact space tailored by the guiding variable.

# 2 Related works

The inherent ambiguity of unsupervised clustering often necessitates guidance to discover partitions that are relevant to a specific analytical context. The guiding principle can be broadly categorized based on its source: either internal to the data itself or provided by an external, contextual variable.
Internal guidance seeks to find structure by modifying the input data **x**. The motivation is that not all features are equally relevant. This has led to variable selection techniques in model-based clustering that identify a feature subset to best reveal latent groups [Raftery and Dean, 2006, Maugis et al., 2009, Maugis et al., 2011]. More advanced frameworks, like Multi-Partition Clustering, concurrently discover multiple valid partitions, each defined by a different feature subset [Marbac and Vandewalle, 2019]. While powerful, these methods are fundamentally introspective, defining “interesting” clusters based solely on the statistical properties of the input data **x**.

In contrast, external guidance aligns the clustering process with an extrinsic goal. The necessity of utilizing auxiliary information to resolve the ambiguity of unsupervised representation learning has been formally established in frameworks such as Identifiable VAEs (iVAE) [Khemakhem et al., 2020]. However, while iVAE leverages auxiliary variables to recover continuous latent factors (disentanglement), our focus is on discovering discrete structures. In this context, the most direct form of guidance is semi-supervised clustering, where the discovery is steered by observing a subset of the actual cluster labels [Kingma et al., 2014, Maaløe et al., 2017] or pairwise constraints [Basu et al., 2008]. However, this approach inherently requires partial access to the ground truth partition, which is often unavailable in exploratory settings. When the cluster labels are unobserved, guidance must instead come from a contextual variable.

The most explicit form of this approach, and the focus of this work, is predictive clustering [Chickering et al., 2000]. This paradigm reframes clustering as a tool for supervised discovery, where the objective is to cluster **x** into groups that are predictive of a target variable **y**. The pursuit of this objective has led to the development of diverse methods. For instance, the Predictive Clustering Trees (PCTs), which adapt decision tree algorithms to recursively partition the joint data space [Blockeel et al., 2000, Stepišnik and Kocev, 2021]. Another
distinct, probabilistic approach is the family of Finite Mixture of Regression (FMR) models, which assume the data arises from subpopulations each governed by its own regression model linking **x** to **y** [Khalili and Chen, 2007, Marbac et al., 2022].

We distinguish our proposed guided clustering paradigm from standard predictive clustering, where the primary metric is prediction accuracy and low performance is typically viewed as a failure. In contrast, guided clustering (our proposal) views the guiding variable not as a target to be perfectly predicted, but as a lens to steer the discovery process. From this perspective, the goal is meaningful organization: even a weak predictor can successfully structure the latent space into clusters that are relevant to the context defined by **y**.

Despite their conceptual differences, these methods have a common limitation: they operate directly in the raw feature space **x**. This direct approach struggles with the curse of dimensionality and fails to capture complex non-linearities. Specifically, PCTs are constrained by greedy, axis-aligned splits, while FMRs are limited by strong parametric assumptions that are often too rigid for real-world data. This challenge motivates a paradigm shift towards deep representation learning, mapping data **x** to a compact latent space **z** to overcome the limitations of operating in the raw data space. More specifically, we adopt a deep generative modeling perspective. The central idea is to discover latent clusters in **x** that are generative of the target **y**. This requires a framework capable of learning a compact, non-linear representation of the data that is structured for this generative-predictive task.

The Variational Autoencoder (VAE) [Kingma and Welling, 2014] serves as a cornerstone for such representation learning, inspiring a diverse range of models. Models like VaDE [Jiang et al., 2016] and GMVAE [Dilokthanakul et al., 2016] learn latent representations for clustering. However, their objective is to reconstruct **x**, meaning the latent space is structured to preserve information about **x** alone. Architectures like Conditional VAEs (CVAEs) [Sohn et al., 2015], Multimodal VAEs (MVAEs) [Wu and Goodman, 2018], and
Characteristic Capturing VAEs (CCVAEs) [Joy et al., 2020] successfully incorporate external information; however, it is used to condition the generation of $x$ or to learn a joint space across modalities.

No existing paradigm fully synthesizes deep representation learning with the specific guided clustering goal. To bridge this gap, we draw inspiration from the Variational Information Bottleneck (VIB) [Alemi et al., 2017]. This framework provides a principled way to extract a representation that is maximally relevant to a target variable while filtering out unrelated noise. However, VIB views the latent space as a compressed encoding computed downstream from the input. Our approach adapts this information-theoretic balance into a deep generative framework. Instead of merely filtering the input, we assume the latent clusters are the pre-existing underlying structure that generates the guiding variable. This shift allows us to use the input not just to predict a target, but to uncover the latent groups that caused it. We realize this paradigm with the so-called Guided Clustering Variational AutoEncoder (GCVAE), which structures the latent space as a probabilistic mixture of components. By optimizing this generative structure to be maximally informative about the guiding variable, GCVAE forces the representation to be both predictively relevant and organized into meaningful clusters.

# 3 Contribution

Our primary contributions are:

* We propose and formalize guided clustering, a framework that reframes the role of contextual information in unsupervised learning. It shifts the objective from using clusters to predict a variable, to using a variable which actively guides the discovery of meaningful data partitions.
* We introduce the Guided Clustering Variational Autoencoder (GCVAE), a deep generative model that operationalizes our paradigm. Its core technical innovation is the integration of a probabilistic mixture structure directly into a guided variational framework, forcing the model to learn a latent space that is both clustered and maximally informative about the guiding variable.

* Empirical validation on open-source and private real-world datasets, demonstrating that our model discovers coherent clusters, even in high-dimensional and non-linear settings.

The remainder of the paper is organized as follows. Section 4 details the generative process, inference model, and objective function of GCVAE. In Section 5, we present our comprehensive experimental evaluation, and finally, we conclude in Section 6.

# 4 The variational implementation of the model

To operationalize the guided clustering paradigm, we propose the Guided Clustering Variational Autoencoder (GCVAE), a deep generative model designed to discover a latent cluster structure in data **x** that is maximally informative for a guiding variable **y**.

The model architecture is built on two core principles. First, it uses an inference model (encoder) to learn a compressed latent representation **z** of the input **x**. This representation is optimized to act as an information bottleneck, retaining only the information from **x** that is necessary for the second component: a generative model (decoder) that predicts the guiding variable **y** from **z**.

Crucially, the latent space is not unstructured. We model its distribution as a Gaussian Mixture Model (GMM). This imposes a distinct cluster structure, forcing the informative bottleneck **z** to be organized into a discrete mixture of components **c**.
Formally, for an input $\mathbf{x} = (\mathbf{x}_1, ..., \mathbf{x}_n)$, with $\mathbf{x}_i \in \mathbb{R}^{d_x}$, and a corresponding guiding variable $\mathbf{y} = (\mathbf{y}_1, ..., \mathbf{y}_n)$, with $\mathbf{y}_i \in \mathbb{R}^{d_y}$, our model learns a continuous latent variable $\mathbf{z} = (\mathbf{z}_1, ..., \mathbf{z}_n)$, with $\mathbf{z}_i \in \mathbb{R}^J$, and infers a discrete cluster assignment $\mathbf{c} = (c_1, ..., c_n)$ with $c_i \in \{1, ..., K\}$, with $n$ the number of observations. We assume that $(\mathbf{x}_i, \mathbf{y}_i, \mathbf{z}_i, c_i)$ are independent and identically distributed. For simplicity, we will omit the index $i$ in the rest of this paper when there is no ambiguity.

Figure 1 provides a graphical overview of this architecture. The following subsections detail the probabilistic formulation of the generative process, the inference model, and the final training objective.

<figure><img src="image_3.png" /><figcaption>Figure 1: Overview of the GCVAE model. Diamond-shaped units denote latent variables, while round ones denote observations.</figcaption></figure>

## 4.1 Generative model

The generative process of our model consists of the following steps:

1. Choose a cluster $c$ using a Multinomial distribution:

$$p_{\pi}(c) = \mathrm{Cat}(c; \pi). \qquad (1)$$

2. Generate a latent vector $\mathbf{z}$ conditioned on the cluster $c$ using a spherical Gaussian
distribution:

$$p_{\mu_c, \sigma_c}(z|c) = \mathcal{N}(z; \mu_c, \sigma_c^2 I). \qquad (2)$$

3. Generate the variable $\mathbf{y}$ from the latent vector $\mathbf{z}$:

$$p_{\theta}(\mathbf{y}|\mathbf{z}) = \mathcal{N}(\mathbf{y}; f_{\theta}(\mathbf{z}), \mathbf{I}). \qquad (3)$$

Here, $\mathbf{\pi} = (\pi_1, ..., \pi_K) \in [0, 1]^K$, with $\pi_c$ the prior probability for cluster $c$, $\sum_{c=1}^K \pi_c = 1$, $\mu_c = \{\mu_{c,j}\}_{j=1,...,J}$ and $\sigma_c = \{\sigma_{c,j}\}_{j=1,...,J}$ are respectively the mean and the diagonal covariance of the multivariate normal distribution corresponding to cluster $c$, $J$ is the dimension of the latent space, $\mathbf{I}$ is an identity matrix, $f_{\theta}(\mathbf{z})$ is a network with input $\mathbf{z}$ and parametrized by $\mathbf{\theta}$.

<figure><img src="image_4.png" /><figcaption>Figure 2: A graphical representation of the generative process, with the GMM in Equation 2 (left) and the decoder in Equation 3 (right).</figcaption></figure>

Following the structure of established deep generative clustering models [Jiang et al., 2016, Dilokthanakul et al., 2016], we assume: $p_{\theta}(\mathbf{y}|\mathbf{z}, c) = p_{\theta}(\mathbf{y}|\mathbf{z})$. According to the generative process above, illustrated in Figure 2, the joint probability of the model is given by:

$$p_{\Theta}(\mathbf{z}, \mathbf{c}, \mathbf{y}) = p_{\Theta}(\mathbf{y}|\mathbf{z})p_{\Theta}(\mathbf{z}|\mathbf{c})p_{\Theta}(\mathbf{c}) \qquad (4)$$

with $\Theta = \{\mu_1, ..., \mu_K, \sigma_1, ..., \sigma_K, \pi_1, ..., \pi_K, \theta\}$.
## 4.2 Inference model

To estimate the generative parameters $\Theta$, a standard VAE procedure would typically infers latent variables from the same data it intends to reconstruct, seeking the posterior $p_{\Theta}(\mathbf{z}, \mathbf{c}|\mathbf{y})$. In contrast, the guided clustering paradigm is built on an informational asymmetry: the encoder must learn to map the input features $\mathbf{x}$ to a latent space that mimics the structure of the ideal posterior conditioned on $\mathbf{y}$, as illustrated in Figure 1.

We therefore define the inference task as finding the variational approximation $q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})$ [Blei et al., 2017] that best captures this ideal structure. Ideally, we would minimize the Kullback-Leibler divergence between our variational approximation and this target posterior, which gives us the central hypothesis of the guided clustering:

$$ \min_{\phi} \mathrm{KL} [q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x}) || p_{\Theta}(\mathbf{z}, \mathbf{c}|\mathbf{y})]. \qquad (5) $$

For the implementation, we illustrate the encoder structure in Figure 3. We assume a classical mean-field approximation for the variational distribution:

$$ q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x}) = q_{\phi}(\mathbf{z}|\mathbf{x})q(\mathbf{c}|\mathbf{x}) \qquad (6) $$

with

$$ q_{\phi}(\mathbf{z}_i|\mathbf{x}_i) = \mathcal{N}(\mathbf{z}_i; \tilde{\mu}_i, \tilde{\sigma}_i^2 I), \qquad (7) $$

$$ [\tilde{\mu}_i, \log \tilde{\sigma}_i] = g_{\phi}(\mathbf{x}_i) \qquad (8) $$

where $\tilde{\mu}_i \in \mathbb{R}^J$, $\tilde{\sigma}_i \in \mathbb{R}^{J+}$ for $i \in \{1, ..., n\}$. Here, $g_{\phi}$ is a neural network parameterized by $\phi$, and the detail of the estimation of $q(\mathbf{c}|\mathbf{x})$ is described in Section 4.3.
<figure><img src="image_5.png" /><figcaption>Figure 3: A graphical representation of the inference model.</figcaption></figure>

## 4.3 Estimation loss

Minimizing the divergence in Equation 5 is equivalent to maximizing the Evidence Lower Bound (ELBO):

$$ \text{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x}) || p_{\Theta}(\mathbf{z}, \mathbf{c}|\mathbf{y})] = \text{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x}) || p_{\Theta}(\mathbf{z}, \mathbf{c})] - \mathbb{E}_{q_{\phi}(\mathbf{z},\mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] + \log p_{\Theta}(\mathbf{y}) \\ \Rightarrow \log p_{\Theta}(\mathbf{y}) \ge \text{ELBO}_{\Theta,\phi}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) = \mathbb{E}_{q_{\phi}(\mathbf{z},\mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] - \text{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x}) || p_{\Theta}(\mathbf{z}, \mathbf{c})]. $$

**Remark** Minimizing the divergence 5 thus compels the inference model to satisfy two competing goals: it must organize the data **x** into the structural constraints of the mixture prior $p_{\Theta}(\mathbf{z}, \mathbf{c})$, while simultaneously retaining sufficient information to reconstruct the guiding variable **y**.

While this standard formulation provides a rigorous lower bound on the log-likelihood, we want a more flexible control over the information flow. Specifically, we need to regulate how much information from **x** is compressed into the latent structure ($\mathbf{z}, \mathbf{c}$) versus how much is used to predict **y**.

To control this trade-off, we frame the learning problem as a constrained optimization task, following the $\beta$-VAE framework [Higgins et al., 2017]. We condition the inference model on **x** while targeting the reconstruction of **y**. The goal is to maximize the reconstruction quality of **y** while constraining the information capacity of the latent representation learned from **x**, using Kullback-Leibler divergence (KL-divergence). With $\epsilon_{\beta}$ specifying the strength
of the applied constraint, this can be formally stated as:

$$ \max_{\Theta, \phi} \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] \text{ subject to } \mathrm{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{z}, \mathbf{c})] < \epsilon_{\beta}. $$

Rewriting it as a Lagrangian under the KKT conditions [Kuhn and Tucker, 1951], we obtain:

$$ \mathcal{F}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) = \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] - \beta (\mathrm{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{z}, \mathbf{c})] - \epsilon_{\beta}) $$

where the KKT multiplier $\beta$ is the regularisation coefficient.

Still following the $\beta$-VAE article principle and formulation, we can derive the final ELBO:

$$ \mathcal{F}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) \geq \mathcal{L}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) = \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] - \beta \mathrm{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{z}, \mathbf{c})]. $$

By adjusting $\beta$, we can prioritize the discovery of a structured latent space that is not just a passthrough for the data, but an effective information bottleneck. Importantly, the introduction of this weighting factor does not impact the theoretical convergence properties discussed in Section 4.4.

The ELBO can then be developed as below (cf. Appendix A):

$$ \mathcal{L}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) \approx -\frac{1}{L} \sum_{l=1}^{L} \sum_{i=1}^{n} \|\mathbf{y}_i - f_{\Theta}(\mathbf{z}_i^{(l)})\|_2^2 \\ + \beta \sum_{i=1}^{n} \left( -\frac{1}{2} \sum_{c=1}^{K} q(c|\mathbf{x}_i) \sum_{j=1}^{J} \left( \log \sigma_{c,j}^2 + \frac{\tilde{\sigma}_{i,j}^2}{\sigma_{c,j}^2} + \frac{(\tilde{\mu}_{i,j} - \mu_{c,j})^2}{\sigma_{c,j}^2} \right) \right. \\ \left. + \sum_{c=1}^{K} q(c|\mathbf{x}_i) \log \pi_c + \frac{1}{2} \sum_{j=1}^{J} (1 + \log \tilde{\sigma}_{i,j}^2) - \sum_{c=1}^{K} q(c|\mathbf{x}_i) \log q(c|\mathbf{x}_i) \right) $$

with $\mathbf{z}_i^{(l)} \sim \mathcal{N}(\tilde{\mu}_i, \tilde{\sigma}_i^2 I)$, and $[\tilde{\mu}_i, \log \tilde{\sigma}_i^2] = g_{\phi}(\mathbf{x}_i)$ with $\tilde{\mu}_i = \{\tilde{\mu}_{i,j}\}_{j=1,...,J}$ and $\tilde{\sigma}_i = \{\tilde{\sigma}_{i,j}\}_{j=1,...,J}$. Recall that $J$ is the dimension of $\mathbf{z}_i$.

We finally approximate $q(c|\mathbf{x})$ using the SGVB estimator (cf. Appendix B), where $L$ is the number of Monte Carlo samples:

$$ q(c|\mathbf{x}) = \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[p_{\Theta}(c|\mathbf{z})] \approx \frac{1}{L} \sum_{l=1}^{L} \frac{p_{\Theta}(\mathbf{z}^{(l)}|c)p_{\Theta}(c)}{\sum_{c'=1}^{K} p_{\Theta}(\mathbf{z}^{(l)}|c')p_{\Theta}(c')}. $$
To encourage sharp clusters assignment while maintaining the differentiability of the model, we employ the Gumbel-Softmax reparameterization trick. This provides a differentiable approximation to sampling from the categorical distribution $q(c|x)$, thereby ensuring that gradients can flow through the entire model.

## 4.4 Consistency

To establish the theoretical validity of our estimator, we analyze the convergence properties of the GCVAE. We utilize the framework of Generalized Variational Inference (GVI) [Knoblauch et al., 2019], which generalizes standard Bayesian inference to posterior beliefs derived from arbitrary loss functions. Our analysis relies on the frequentist consistency results established in Theorem 2 of [Knoblauch, 2019]. By applying this framework to the specific structure of the GCVAE, we can state the following:

**Proposition 1 (GCVAE consistency).** Let $\psi = \{\Theta, \phi\} \in \Psi$ denote the complete set of learnable global parameters in the GCVAE model. If the data are i.i.d. and $\Psi$ is compact, the estimator $\hat{\psi}_n$ converges to a point mass at the population-optimal parameter $\psi^*$ as number of observations $n \to \infty$.

The proof is available in [Knoblauch, 2019]. The direct application of this theorem is non-trivial due to the non-standard nature of our mixed variational family, which involves both continuous and discrete latent variables. As the rigorous verification of these conditions is rarely detailed in the literature, we provide a full verification in Appendix C.

# 5 Experiments

In this section, we evaluate the performance of the GCVAE on two datasets. The model has been implemented in Python using Pytorch [Paszke et al., 2019], and the code is available
at https://github.com/vcourrier/gcvae.

## 5.1 Practical implementation guidelines

To effectively apply the GCVAE in practice, specific attention must be paid to the architectural choices and the optimization of the ELBO. While the generative framework is agnostic to the specific family of neural networks used, the choice of hyperparameters impacts the model’s ability to uncover meaningful cluster structures.

### Initializing GMM
In this work, pre-training is used to initialize GMM parameters $(\pi_c, \mu_c, \sigma_c)$, a common practice in deep clustering [Jiang et al., 2016, Xie et al., 2016]. We pretrain the model without the clusters in the latent space, leading to the following ELBO:

$$ \mathcal{L}_{\Theta, \phi, \beta}^{\text{pretrain}}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) = \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] - \beta \text{KL}[q_{\phi}(\mathbf{z}|\mathbf{x})||\mathcal{N}(\mathbf{z}; 0_J, \mathbf{I})] $$

with $0_J$ a vector null of dimension $J$.

After pretraining the encoder and decoder for a few epochs, we can fit a GMM in the latent space to initialize its parameters.

### Network architectures and the information bottleneck
The encoder and decoder functions $(g_{\phi}$ and $f_{\theta})$ can be parameterized using various architectures suited to the data modality, such as Multilayer Perceptrons (MLPs), Convolutional Neural Networks (CNNs), Recurrent Neural Networks (RNNs)... However, from a statistical perspective, we caution against the use of excessively complex networks. If the encoder capacity is too high, the network may “absorb” the data’s structure into its parameters, rendering the latent space $\mathbf{z}$ a simple passthrough rather than a structured representation. To ensure the model learns robust clusters, we design the networks to act as an effective information bottleneck. By intentionally constraining the network capacity (e.g., reducing the number of layers or units), we force the model to compress the essential features of $\mathbf{x}$ into the latent clusters $\mathbf{c}$ and
variable **z** to maximize the predictive accuracy of **y**. This parsimonious approach ensures that the clustering structure captures the signal rather than the network weights.

**The regularization parameter β** The coefficient **β** controls the trade-off between the reconstruction fidelity of the guiding variable and the adherence of the latent space to the GMM prior. In our experiments, we observe that a small weighting (β < 1) typically yields a more tractable optimization problem [Burgess et al., 2018]. This prevents the complex KL-divergence term associated with the GMM from dominating the loss. Multiple techniques can be used to find a suitable **β**, via standard grid search or KL annealing [Fu et al., 2019], where **β** is gradually increased during training.

**Optimization and hyperparameters** The model parameters are estimated using the Adam optimizer [Kingma and Ba, 2014], a standard for stochastic gradient-based optimization. Beyond **β**, the key hyperparameters requiring selection include the learning rate, the dimensions of the neural networks, and the number of latent clusters *K*. The choice of *K* may be driven by prior domain knowledge (as seen in our clinical application) or determined via model selection criteria compatible with the VAE framework.

## 5.2 MNIST-SVHN dataset

<figure><img src="image_6.png" /><figcaption>Figure 4: Example of SVHN (top) and MNIST (bottom) data.</figcaption></figure>

**Dataset** We evaluate our proposed approach on a dataset composed of paired MNIST and black-and-white SVHN images, with the SVHN image as our input **x** and the MNIST
image as our guiding variable **y**, where each pair represents the same digit class, following
the methodology introduced in the MMVAE model [Shi et al., 2019]. Each instance of a
digit class (in either dataset) is randomly paired with 30 instances of the same digit class
from the other dataset. We use the standard training partitions of MNIST and SVHN
(see details in the respective datasets description). As illustrated in Figure 4, the dataset
presents a variety of styles, making the task of capturing the digit classes challenging.

**Implementation details** For this experiment, we use CNNs for the encoder of SVHN, and a MLP for the decoder of MNIST. The latent space dimension is set to 20. For learning, we use the Adam optimizer [Kingma and Ba, 2014]. We set $\beta = 0.1$. During the pre-training, we use a learning rate of 0.0001 and run it for 5 epochs. During the training, we use a learning rate of 0.0005 for the parameters of the encoder and decoder, and 0.0005 too for the parameters of the GMM, and set the number of clusters to 10, run it for 50 epochs.

**Impact of the guiding variable** To assess the influence of the predictive variable on our model's clustering performance, we conduct experiments on a multitude of clustering methods.

Consistent with prior research in image clustering, we evaluate performance using global classification accuracy (ACC), where a cluster-to-class mapping is determined via the Hungarian algorithm [Kuhn, 1955]. Without a guiding variable **y**, with the same model but reconstructing **x**, the model achieves an ACC of 11.6% (average on 10 runs). In contrast, using the MNIST data as the guiding variable significantly improves the ACC to 62.1% (average on 10 runs).

**Comparison on standard benchmarks**
<table><caption>Table 1: ACC on standard clustering benchmarks.</caption><thead><tr><th>MODEL</th><th>ACC</th></tr></thead><tbody><tr><td colspan="2"><em>CLUSTERING MODELS WITH IMAGE-SPECIFIC TRANSFORMATIONS</em></td></tr><tr><td>DTI K-MEANS [MONNIER ET AL., 2020]</td><td>44.5%</td></tr><tr><td>SCAE [KOSIOREK ET AL., 2019]</td><td>55.3%</td></tr><tr><td>DTI GMM [MONNIER ET AL., 2020]</td><td>57.4%</td></tr><tr><td>ACOL-GAR [KILINC AND UYSAL, 2018]</td><td>76.8%</td></tr><tr><td colspan="2"><em>CLUSTERING MODELS WITH DOMAIN-AGNOSTIC DESIGNS</em></td></tr><tr><td>GMM [DEMPSTER ET AL., 1977]</td><td>11.6%</td></tr><tr><td>DEC [XIE ET AL., 2016]</td><td>11.9%</td></tr><tr><td>K-MEANS [MACQUEEN, 1967]</td><td>12.2%</td></tr><tr><td>DEEPCLUSTER-V2 [CARON ET AL., 2020]</td><td>20.6%</td></tr><tr><td>VADE [JIANG ET AL., 2016]</td><td>30.8%</td></tr><tr><td>MFCVAE [FALCK ET AL., 2021]</td><td>56.3%</td></tr><tr><td>IMSAT [HU ET AL., 2017]</td><td>57.3%</td></tr><tr><td><strong>GCVAE</strong> (OUR PROPOSAL)</td><td>62.1%</td></tr></tbody></table>

To quantify the value added by the guiding variable y, we contrast the performance of GCVAE against established unsupervised clustering methods. We categorize these baselines into two groups based on their reliance on domain-specific knowledge (Table 1).

In the first category, models rely on image-specific transformations or architectures, effectively acting as implicit guidance. DTI [Monnier et al., 2020] incorporates spatial or morphological transformations, SCAE [Kosiorek et al., 2019] segments images into part templates before reasoning about “object capsules”, and ACOL-GAR [Kilinc and Uysal, 2018] applies domain-specific transformations to generate pseudo parent classes and achieves high performance. Notably, ACOL-GAR achieves the highest performance on this benchmark
(76.8%), surpassing GCVAE (62.1%). However, this superiority stems from “hard-coded”
domain knowledge: the method relies on invariances specific to visual data to generate super-
vision. While highly effective for images, this reliance makes such methods less transferable
to non-visual tasks (e.g., tabular or sensor data) where such domain-specific invariants are
unknown or undefined.

In the second category, domain-agnostic methods seek broader applicability without depend-
ing on extensive image transformations. This includes classical approaches like GMM
[Dempster et al., 1977] and K-means [MacQueen, 1967], as well as deep methods like
DEC [Xie et al., 2016], VaDE [Jiang et al., 2016], MFCVAE [Falck et al., 2021], IMSAT
[Hu et al., 2017], [Caron et al., 2018] and DeepCluster-v2 [Caron et al., 2020]. As shown in
Table 1, these methods struggle to recover the digit classes solely from raw pixel statistics,
with accuracies ranging from 11.6% to 57.3%. GCVAE achieves 62.1% ACC, surpassing
these domain-agnostic models by leveraging the guiding variable exclusively in the generative
process.

This comparison highlights the specific niche of GCVAE: while it may not outperform
specialized models equipped with extensive domain-specific augmentations (like ACOL-
GAR), it significantly outperforms generic unsupervised methods by effectively leveraging
the guiding variable *y*. This offers a flexible solution that provides structure to complex
data without requiring the manual design of domain-specific transformations.

**Visualization of the clusters during the training** Figure 5 shows the t-SNE [van der Maaten and Hinton, 2008] visualizations of our model’s latent space at different training epochs (1, 5, 10, 25, 50) for 5,000 SVHN training examples. Each point in the plot corresponds to the latent vector of a single image, colored according to its ground-truth label. Initially, after the pretraining step, the latent-space representation is relatively unstructured. By epoch 5, clusters begin to appear, indicating that our model is starting to learn features
<figure><img src="image_7.png" /><figcaption>Figure 5: t-SNE visualisations at different epoch during the training of 5,000 training images.</figcaption></figure>

that help separate different digit classes. As training continues, these clusters become increasingly well-defined, and data points sharing the same label gather into tighter, more distinct regions. Notably, the boundaries between clusters also grow clearer, suggesting that the learned representations reflect class-specific properties more effectively over time. Overall, this progressive separation of clusters highlights how the guided training procedure refines the latent-space vectors over time to achieve better discrimination among different labels in an unsupervised setting.
## 5.3 Sleep dataset

**Dataset** To demonstrate the operational utility of our guided clustering paradigm, we apply GCVAE to a real-life dataset derived from Withings' proprietary data. Withings is a French company in the field of digital health, known for designing and manufacturing a wide range of connected health devices<fnref n="1" />. The experiment is designed to demonstrate how a guiding variable can help recover meaningful and coherent subgroups from complex data. While this dataset allows us to showcase our method's properties, the discovered subgroups are intended as a methodological proof-of-concept rather than a definitive clinical finding. Furthermore, the dataset has not undergone rigorous debiasing procedures.

Our analyses are based on sleep data collected from 50,000 individuals, with each contributing one night of data. All personally identifiable information has been removed in compliance with GDPR guidelines. We use seven features for the input vector **x**: sleep_duration, bmi, age, light_sleep_duration, deep_sleep_duration, nb_sleep_interruptions, and avg_night_hr. The guiding variable **y** explored is the apnea_hypopnea_index, a standard clinical metric that quantifies sleep apnea severity by measuring the number of breathing interruptions per hour of sleep. A detailed overview of the dataset, including feature descriptions and Apnea-Hypopnea Index (AHI) categorization, is provided in Appendix D. We split the dataset into train, test, and validation sets of repartition 70%, 20% and 10% respectively.

**Implementation details** The objective is to demonstrate that using the Apnea-Hypopnea Index (AHI) as a guiding variable **y** allows the model to discover more clinically coherent partitions of the user data **x** than an unguided approach. We set $K = 3$ clusters to search for distinct user phenotypes.

For a principled comparison, we contrast our GCVAE (which learns a mapping $\mathbf{x} \mapsto \mathbf{y}$)

<footnote n="1">¹https://www.withings.com/us/en/</footnote>
with a unguided baseline. This baseline uses the same architecture but is adapted for a different task. It receives the concatenated input $(x, y)$ and is trained to find a clustered representation that reconstructs both variables. This comparison allows us to differentiate between finding structure in the joint $(x, y)$ space versus finding structure within $x$ that is relevant to $y$.

For this experiment, we use MLPS in the encoder and the decoder. For learning, we use the Adam optimizer, as before. During the pre-training (5 epochs), we use a learning rate of 0.0005 and a $\beta$ of 0.001. During the training (50 epochs), we use a learning rate of 0.0001 for the parameters of the encoder and decoder, and 0.00001 for the parameters of the GMM, and a $\beta$ of 0.01.

<figure><img src="image_8.png" /><figcaption>Figure 6: Distribution of the AHI within each cluster discovered by the unguided model. Note the significant overlap between the distributions for Cluster 1 and Cluster 2.</figcaption></figure>

**Results and analysis** The unguided baseline, with direct access to AHI as an input, partitions the data based on its most prominent features. An analysis of the cluster profiles (see Appendix E for the full table) shows that the model excels at isolating a high-level of AHI cohort in Cluster 3, a task simplified by the direct visibility of $y$. For the remaining population, however, the model does not find further AHI-related structure. It defaults to
partitioning users based on the next largest source of variance: sleep duration. Consequently, Clusters 1 and 2 have nearly identical AHI profiles but represent different sleep behaviors. The AHI distributions in Figure 6 visually confirm this: while Cluster 3 is distinct, the distributions for Clusters 1 and 2 are almost perfectly overlapping.

In contrast, GCVAE uses $\boldsymbol{y}$ not as an input to be partitioned, but as a lens to find the most meaningful structure within $\boldsymbol{x}$. It is tasked with discovering groups of patients whose features in $\boldsymbol{x}$ are collectively indicative of different AHI levels. The resulting cluster profiles (see Appendix E for detailed profiles) reveal three coherent subgroups that align with a known clinical gradient of risk:

* Cluster 1 (Low-Risk): A group whose profile in $\boldsymbol{x}$ (lowest age and BMI, most deep sleep) corresponds to a low AHI.

* Cluster 2 (Intermediate-Risk): A transitional group whose profile in $\boldsymbol{x}$ indicates a moderately increased clinical risk.

* Cluster 3 (High-Risk): A cohort whose profile in $\boldsymbol{x}$ (highest age and BMI, least deep sleep) corresponds to a high AHI.

<figure><img src="image_9.png" /><figcaption>Figure 7: Distribution of the AHI within each cluster discovered by the GCVAE model.</figcaption></figure>
Figure 7 visualizes this result. The AHI distributions are clearly ordered, though they exhibit overlap. This overlap is a direct reflection of the challenging inference task. Unlike the baseline, which directly observes **y**, GCVAE must infer AHI severity from complex patterns in **x**. The overlap thus represents the inherent uncertainty of this relationship in real-world data. Crucially, this clustering is achieved solely from the sleep biometrics **x**; unlike the baseline, the GCVAE does not have access to the guiding variable **y** (AHI) during inference. The key achievement is that the model successfully learned to use the guidance from **y** to organize **x** into a coherent and clinically-aligned structure.

## 5.4 Data availability statement

The MNIST and SVHN datasets analyzed during the current study are available in the public domain via the torchvision library. The sleep dataset used in this study is proprietary to Withings and is not publicly available due to participant privacy and GDPR restrictions; however, details are provided in the article in Appendix D.

# 6 Conclusion

We introduced **guided clustering**, a paradigm that formalizes the implicit guidance inherent in any practical clustering analysis. We argue that for discovery to be meaningful, the analytical context must be explicitly integrated into the model's optimization objective, rather than being relegated to an informal, trial-and-error process. The Guided Clustering Variational Autoencoder (GCVAE) is presented as one effective deep generative realization of this principle, though the paradigm itself is model-agnostic.

As an implementation of this paradigm, we proposed the GCVAE. We demonstrated that by optimizing a latent representation to be maximally informative about a guiding variable, GCVAE discovers clusters that are not only coherent but also meaningful for a given
analytical task. Our experiments on both public and proprietary datasets confirmed its
ability to uncover relevant structures, outperforming domain-agnostic methods and revealing
clinically coherent user subgroups.

This paradigm opens several avenues for future work. The GCVAE itself can be extended
by replacing its encoder and decoder with more sophisticated models, such as Graph Neural
Networks or Transformers, to apply guided clustering to graph-structured or sequential
data.

By formally incorporating a guiding signal into the optimization objective, our approach
provides a principled method for navigating the inherent ambiguity of what constitutes a
“good” cluster, ensuring the discovered partitions are aligned with a specific analytical goal.
We believe this formalization of guidance is a promising step toward developing models that
can discover relevant, human-interpretable structures within vast, unannotated datasets.

References

[Alemi et al., 2017] Alemi, A. A., Fischer, I., Dillon, J. V., and Murphy, K. (2017). Deep variational information bottleneck. *International Conference on Learning Representations*.

[Basu et al., 2008] Basu, S., Davidson, I., and Wagstaff, K., editors (2008). *Constrained Clustering: Advances in Algorithms, Theory, and Applications*. Chapman and Hall/CRC, New York.

[Blei et al., 2017] Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017). Variational inference: A review for statisticians. *Journal of the American Statistical Association*, 112(518):859–877.

[Blockeel et al., 2000] Blockeel, H., De Raedt, L., and Ramon, J. (2000). Top-down induction of clustering trees. *Proc. 15th Intl. Conf. on Machine Learning*.
[Burgess et al., 2018] Burgess, C. P., Higgins, I., Pal, A., Matthey, L., Watters, N., Desjardins, G., and Lerchner, A. (2018). Understanding disentangling in β-vae. *arXiv preprint arXiv:1804.03599*.

[Caron et al., 2018] Caron, M., Bojanowski, P., Joulin, A., and Douze, M. (2018). Deep clustering for unsupervised learning of visual features. In *Proceedings of the European conference on computer vision (ECCV)*, pages 132–149.

[Caron et al., 2020] Caron, M., Misra, I., Mairal, J., Goyal, P., Bojanowski, P., and Joulin, A. (2020). Unsupervised learning of visual features by contrasting cluster assignments. *Advances in neural information processing systems*, 33:9912–9924.

[Chickering et al., 2000] Chickering, D. M., Heckerman, D., Meek, C., Platt, J. C., and Thiesson, B. (2000). Goal-oriented clustering. *Technical Report, MSR-TR-200-82*.

[Dempster et al., 1977] Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society: Series B*, 39:1–38.

[Dilokthanakul et al., 2016] Dilokthanakul, N., Mediano, P. A. M., Garnelo, M., Lee, M. C. H., Salimbeni, H., Arulkumaran, K., and Shanahan, M. (2016). Deep unsupervised clustering with gaussian mixture variational autoencoders. *CoRR*, abs/1611.02648.

[Falck et al., 2021] Falck, F., Zhang, H., Willetts, M., Nicholson, G., Yau, C., and Holmes, C. C. (2021). Multi-facet clustering variational autoencoders. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W., editors, *Advances in Neural Information Processing Systems*, volume 34, pages 8676–8690. Curran Associates, Inc.

[Fu et al., 2019] Fu, H., Li, C., Liu, X., Gao, J., Celikyilmaz, A., and Carin, L. (2019). Cyclical annealing schedule: A simple approach to mitigating kl vanishing. *arXiv preprint arXiv:1903.10145*.
[Higgins et al., 2017] Higgins, I., Matthey, L., Pal, A., Burgess, C., Glorot, X., Botvinick, M., Mohamed, S., and Lerchner, A. (2017). beta-VAE: Learning basic visual concepts with a constrained variational framework. In *International Conference on Learning Representations*.

[Hu et al., 2017] Hu, W., Miyato, T., Tokui, S., Matsumoto, E., and Sugiyama, M. (2017). Learning discrete representations via information maximizing self-augmented training. In *International conference on machine learning*, pages 1558–1567. PMLR.

[Jiang et al., 2016] Jiang, Z., Zheng, Y., Tan, H., Tang, B., and Zhou, H. (2016). Variational deep embedding: A generative approach to clustering. *CoRR*, abs/1611.05148.

[Joy et al., 2020] Joy, T., Schmon, S. M., Torr, P. H. S., Siddharth, N., and Rainforth, T. (2020). Capturing label characteristics in vaes. In *International Conference on Learning Representations*.

[Khalili and Chen, 2007] Khalili, A. and Chen, J. (2007). Variable selection in finite mixture of regression models. *Journal of the American Statistical Association*, 102(479):1025–1038.

[Khemakhem et al., 2020] Khemakhem, I., Kingma, D., Monti, R., and Hyvarinen, A. (2020). Variational autoencoders and nonlinear ica: A unifying framework. In *International conference on artificial intelligence and statistics*, pages 2207–2217. PMLR.

[Kilinc and Uysal, 2018] Kilinc, O. and Uysal, I. (2018). Learning latent representations in neural networks for clustering through pseudo supervision and graph-based activity regularization. In *International Conference on Learning Representations*.

[Kingma and Ba, 2014] Kingma, D. and Ba, J. (2014). Adam: A method for stochastic optimization. *International Conference on Learning Representations*.

[Kingma et al., 2014] Kingma, D. P., Rezende, D. J., Mohamed, S., and Welling, M. (2014). Semi-supervised learning with deep generative models. *Advances in neural information*
processing systems, 27.

[Kingma and Welling, 2014] Kingma, D. P. and Welling, M. (2014). Auto-encoding variational bayes. In Bengio, Y. and LeCun, Y., editors, *ICLR*.

[Kirk et al., 2017] Kirk, V., Baughn, J., D'Andrea, L., Friedman, N., Galion, A., Garetz, S., Hassan, F., Wrede, J., Harrod, C. G., and Malhotra, R. K. (2017). American academy of sleep medicine position paper for the use of a home sleep apnea test for the diagnosis of osa in children. *Journal of Clinical Sleep Medicine*, 13(10):1199-1203.

[Knoblauch, 2019] Knoblauch, J. (2019). Frequentist consistency of generalized variational inference. *arXiv preprint arXiv:1912.04946*.

[Knoblauch et al., 2019] Knoblauch, J., Jewson, J., and Damoulas, T. (2019). Generalized variational inference: Three arguments for deriving new posteriors. *arXiv preprint arXiv:1904.02063*.

[Kosiorek et al., 2019] Kosiorek, A., Sabour, S., Teh, Y. W., and Hinton, G. E. (2019). Stacked capsule autoencoders. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc.

[Kuhn, 1955] Kuhn, H. W. (1955). The Hungarian Method for the Assignment Problem. *Naval Research Logistics Quarterly*, 2(1-2):83-97.

[Kuhn and Tucker, 1951] Kuhn, H. W. and Tucker, A. W. (1951). Nonlinear Programming. In *Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability*, volume 2, pages 481-493. University of California Press.

[Maaløe et al., 2017] Maaløe, L., Fraccaro, M., and Winther, O. (2017). Semi-supervised generation with cluster-aware generative models. *arXiv preprint arXiv:1704.00637*.
[MacQueen, 1967] MacQueen, J. B. (1967). Some methods for classification and analysis of multivariate observations. In Cam, L. M. L. and Neyman, J., editors, *Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability*, volume 1, pages 281–297. University of California Press.

[Marbac et al., 2022] Marbac, M., Sedki, M., Biernacki, C., and Vandewalle, V. (2022). Simultaneous Semiparametric Estimation of Clustering and Regression. *Journal of Computational and Graphical Statistics*, 31(2):477–485. Publisher: Informa UK Limited.

[Marbac and Vandewalle, 2019] Marbac, M. and Vandewalle, V. (2019). A tractable multi-partitions clustering. *Computational Statistics & Data Analysis*, 132:167–179.

[Maugis et al., 2009] Maugis, C., Celeux, G., and Martin-Magniette, M.-L. (2009). Variable selection for clustering with Gaussian mixture models. *Biometrics*, 65(3):701–709.

[Maugis et al., 2011] Maugis, C., Celeux, G., and Martin-Magniette, M.-L. (2011). Variable selection in model-based discriminant analysis. *Journal of Multivariate Analysis*, 102(10):1374–1387.

[Monnier et al., 2020] Monnier, T., Groueix, T., and Aubry, M. (2020). Deep transformation-invariant clustering. *Advances in neural information processing systems*, 33:7945–7955.

[Paszke et al., 2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019). Pytorch: An imperative style, high-performance deep learning library. *Advances in neural information processing systems*, 32.

[Raftery and Dean, 2006] Raftery, A. E. and Dean, N. (2006). Variable selection for model-based clustering. *Journal of the American Statistical Association*, 101(473):168–178.

[Shi et al., 2019] Shi, Y., Paige, B., Torr, P., et al. (2019). Variational mixture-of-experts autoencoders for multi-modal deep generative models. *Advances in neural information processing systems*, 31(2):477–485.
processing systems, 32.

[Sohn et al., 2015] Sohn, K., Lee, H., and Yan, X. (2015). Learning structured output representation using deep conditional generative models. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R., editors, *Advances in Neural Information Processing Systems*, volume 28. Curran Associates, Inc.

[Stepišnik and Kocev, 2021] Stepišnik, T. and Kocev, D. (2021). Oblique predictive clustering trees. *Knowledge-Based Systems*, 227:107228.

[van der Maaten and Hinton, 2008] van der Maaten, L. and Hinton, G. (2008). Visualizing data using t-SNE. *Journal of Machine Learning Research*, 9:2579–2605.

[Wu and Goodman, 2018] Wu, M. and Goodman, N. (2018). Multimodal generative models for scalable weakly-supervised learning. *Advances in neural information processing systems*, 31.

[Xie et al., 2016] Xie, J., Girshick, R., and Farhadi, A. (2016). Unsupervised deep embedding for clustering analysis. In *International conference on machine learning*, pages 478–487. PMLR.
# SUPPLEMENTARY MATERIAL

## A Development of the ELBO

### A.1 Terms of the ELBO

Using Equation 6, we can decompose the ELBO in five terms:

$$ 
\begin{aligned}
\mathcal{L}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) &= \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] - \beta \text{KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{z}|\mathbf{c})p_{\Theta}(\mathbf{c})] \\
&= \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] + \beta (\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{z}|\mathbf{c})] \\
&\quad + \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{c})] - \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log q_{\phi}(\mathbf{z}|\mathbf{x})] - \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log q(\mathbf{c}|\mathbf{x})]).
\end{aligned}
$$

1st term

$$ 
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] &= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log p_{\Theta}(\mathbf{y}|\mathbf{z}) d\mathbf{z} \\
&= \int q_{\phi}(\mathbf{z}|\mathbf{x}) \log p_{\Theta}(\mathbf{y}|\mathbf{z}) d\mathbf{z} \\
&= \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] \\
&= \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})} \left[ \log \left( \frac{1}{\sqrt{(2\pi)^{d_y}}} \exp(-\frac{1}{2}||\mathbf{y} - f_{\Theta}(\mathbf{z})||^2) \right) \right] \\
&= -\frac{d_y}{2} \log(2\pi) - \frac{1}{2} \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})} [||\mathbf{y} - f_{\Theta}(\mathbf{z})||^2].
\end{aligned}
$$

Using the SGVB estimator, we can approximate it as:

$$ \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[\log p_{\Theta}(\mathbf{y}|\mathbf{z})] \propto -\frac{1}{L} \sum_{l=1}^{L} \sum_{i=1}^{n} ||\mathbf{y}_i - f_{\Theta}(\mathbf{z}_i^{(l)})||_2^2 $$

with $\mathbf{z}_i^{(l)} \sim \mathcal{N}(\mathbf{z}; \tilde{\mu}_i, \tilde{\sigma}_i^2 I)$, and $[\tilde{\mu}_i, \log \tilde{\sigma}_i^2] = g_{\phi}(\mathbf{x}_i)$. L is the number of Monte Carlo samples in the SGVB estimator.
2nd term

$$
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{z}|\mathbf{c})] &= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log p_{\Theta}(\mathbf{z}|\mathbf{c}) d\mathbf{z} \\
&= \sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \int \mathcal{N}(\mathbf{z}; \tilde{\mu}, \tilde{\sigma}^2 I) \log \mathcal{N}(\mathbf{z}; \mu_c, \sigma_c^2 I) d\mathbf{z}.
\end{aligned}
$$

Using Lemma A.1 in Appendix A.2, we have:

$$
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{z}|\mathbf{c})] &= -\sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \left[ \frac{J}{2} \log(2\pi) + \frac{1}{2} \sum_{j=1}^{J} \left( \log \sigma_{cj}^2 - \frac{\tilde{\sigma}_j^2}{\sigma_{cj}^2} - \frac{(\tilde{\mu}_j - \mu_{cj})^2}{\sigma_{cj}^2} \right) \right] \\
&= -\frac{J}{2} \log(2\pi) - \frac{1}{2} \sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \sum_{j=1}^{J} \left( \log \sigma_{cj}^2 + \frac{\tilde{\sigma}_j^2}{\sigma_{cj}^2} + \frac{(\tilde{\mu}_j - \mu_{cj})^2}{\sigma_{cj}^2} \right).
\end{aligned}
$$

3rd term

$$
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\Theta}(\mathbf{c})] &= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log p_{\Theta}(\mathbf{c}) d\mathbf{z} \\
&= \sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \log \pi_{\mathbf{c}}.
\end{aligned}
$$

4th term

$$
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log q_{\phi}(\mathbf{z}|\mathbf{x})] &= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log q_{\phi}(\mathbf{z}|\mathbf{x}) d\mathbf{z} \\
&= \int \mathcal{N}(\mathbf{z}; \tilde{\mu}, \tilde{\sigma}^2 I) \log \mathcal{N}(\mathbf{z}; \tilde{\mu}, \tilde{\sigma}^2 I) d\mathbf{z}.
\end{aligned}
$$

Using Lemma A.1 in Appendix A.2, we have:

$$
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log q_{\phi}(\mathbf{z}|\mathbf{x})] = -\frac{J}{2} \log(2\pi) - \frac{1}{2} \sum_{j=1}^{J} (1 + \log \tilde{\sigma}_j^2).
$$

5th term

$$
\begin{aligned}
\mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log q(\mathbf{c}|\mathbf{x})] &= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log q(\mathbf{c}|\mathbf{x}) d\mathbf{z} \\
&= \sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \log q(\mathbf{c}|\mathbf{x}).
\end{aligned}
$$
## A.2 Lemma

As presented in [Jiang et al., 2016], we have Lemma A.1.

**Lemma A.1.** Given two multivariate Gaussian distributions $q(z) = N(z; \tilde{\mu}, \tilde{\sigma}^2 I)$ and $p(z) = N(z; \mu, \sigma^2 I)$, we have:

$$ \int q(z) \log p(z) dz = -\frac{J}{2} \log(2\pi) - \frac{1}{2} \sum_{j=1}^{J} \left( \log \sigma_j^2 + \frac{\tilde{\sigma}_j^2}{\sigma_j^2} + \frac{(\tilde{\mu}_j - \mu_j)^2}{\sigma_j^2} \right) $$

where $\mu_j, \sigma_j, \tilde{\mu}_j$ and $\tilde{\sigma}_j$ simply denote the $j^{th}$ element of $\mu, \sigma, \tilde{\mu}$ and $\tilde{\sigma}$, respectively, and $J$ is the dimensionality of $z$.

*Proof.*

$$
\begin{aligned}
\int q(z) \log p(z) dz &= \int N(z; \tilde{\mu}, \tilde{\sigma}^2 I) \log N(z; \mu, \sigma^2 I) dz \\
&= \int \prod_{j=1}^{J} \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \log\left[\prod_{j=1}^{J} \frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left(-\frac{(z_j - \mu_j)^2}{2\sigma_j^2}\right)\right] dz \\
&= \sum_{j=1}^{J} \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \log\left[\frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left(-\frac{(z_j - \mu_j)^2}{2\sigma_j^2}\right)\right] dz_j \\
&= \sum_{j=1}^{J} \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \left[-\frac{1}{2}\log(2\pi\sigma_j^2)\right] dz_j \\
&\quad - \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \frac{(z_j - \mu_j)^2}{2\sigma_j^2} dz_j \\
&= \sum_{j=1}^{J} -\frac{1}{2}\log(2\pi\sigma_j^2) \\
&\quad - \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \frac{(z_j - \tilde{\mu}_j)^2 + 2(z_j - \tilde{\mu}_j)(\tilde{\mu}_j - \mu_j) + (\tilde{\mu}_j - \mu_j)^2 \tilde{\sigma}_j^2}{2\tilde{\sigma}_j^2} dz_j \\
&= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2} dz_j \\
&\quad - \underbrace{\frac{(\tilde{\mu}_j - \mu_j)}{2\sigma_j^2} \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) (z_j - \tilde{\mu}_j) dz_j}_{=0} \\
&\quad - \underbrace{\frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2} \int \frac{1}{\sqrt{2\pi\tilde{\sigma}_j^2}} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) dz_j}_{=1}
\end{aligned}
$$
$$= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \int \frac{1}{\sqrt{2\pi}\tilde{\sigma}_j^2} \exp\left(-\frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2}\right) \frac{(z_j - \tilde{\mu}_j)^2}{2\tilde{\sigma}_j^2} dz_j - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2}$$

$$= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \int \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x_j^2}{2}\right) \frac{x_j^2}{2} dx_j - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2} \quad \text{By change of variables}$$

$$= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \int \frac{1}{\sqrt{2\pi}} \left(-\frac{x_j}{2}\right) d\left(\exp\left(-\frac{x_j^2}{2}\right)\right) - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2} \quad \text{By change of variables}$$

$$= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \left[ \left[ \frac{1}{\sqrt{2\pi}} \left(-\frac{x_j}{2}\right) \exp\left(-\frac{x_j^2}{2}\right) \right]_{-\infty}^{+\infty} - \int \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x_j^2}{2}\right) d\left(-\frac{x_j}{2}\right) \right] - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2}$$

By integration by parts

$$= C - \frac{\tilde{\sigma}_j^2}{\sigma_j^2} \left[ (0 - 0) - \left(-\frac{1}{2}\right) \underbrace{\int \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x_j^2}{2}\right) dx_j}_{=1} \right] - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2}$$

$$= \sum_{j=1}^{J} -\frac{1}{2} \log (2\pi\sigma_j^2) - \frac{\tilde{\sigma}_j^2}{2\sigma_j^2} - \frac{(\tilde{\mu}_j - \mu_j)^2}{2\sigma_j^2}$$

$$= -\frac{J}{2} \log(2\pi) - \frac{1}{2} \sum_{j=1}^{J} \left( \log \sigma_j^2 + \frac{\tilde{\sigma}_j^2}{\sigma_j^2} + \frac{(\tilde{\mu}_j - \mu_j)^2}{\sigma_j^2} \right)$$

where $C$ denotes $\sum_{j=1}^{J} -\frac{1}{2} \log (2\pi\sigma_j^2)$ for simplicity. $\square$

<figure><img src="image_10.png" /></figure>

## B Approximate $q(c|x)$

We describe how to formulate $q(c|x)$ to maximize the ELBO. Specifically, our objective can be rewritten as:

$$\mathcal{L}_{\Theta, \phi, \beta}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c}) = \mathbb{E}_{q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})}[\log p_{\theta}(\mathbf{y}|\mathbf{z}, \mathbf{c})] - \beta \text{ KL}[q_{\phi}(\mathbf{z}, \mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{z}, \mathbf{c})]$$

$$= \int \sum_{\mathbf{c}} q_{\phi}(\mathbf{z}|\mathbf{x}) q(\mathbf{c}|\mathbf{x}) \log \frac{p_{\Theta}(\mathbf{y}|\mathbf{z}) p_{\Theta}(\mathbf{c}|\mathbf{z})^{\beta} p_{\Theta}(\mathbf{z})^{\beta}}{q_{\phi}(\mathbf{z}|\mathbf{x})^{\beta} q(\mathbf{c}|\mathbf{x})^{\beta}} d\mathbf{z}$$

$$= \int q_{\phi}(\mathbf{z}|\mathbf{x}) \log \frac{p_{\Theta}(\mathbf{y}|\mathbf{z}) p_{\Theta}(\mathbf{z})^{\beta}}{q_{\phi}(\mathbf{z}|\mathbf{x})^{\beta}} d\mathbf{z} + \int q_{\phi}(\mathbf{z}|\mathbf{x}) \sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) \log \frac{p_{\Theta}(\mathbf{c}|\mathbf{z})^{\beta}}{q(\mathbf{c}|\mathbf{x})^{\beta}} d\mathbf{z}$$

$$= \int q_{\phi}(\mathbf{z}|\mathbf{x}) \log \frac{p_{\Theta}(\mathbf{y}|\mathbf{z}) p_{\Theta}(\mathbf{z})^{\beta}}{q_{\phi}(\mathbf{z}|\mathbf{x})^{\beta}} d\mathbf{z} - \int q_{\phi}(\mathbf{z}|\mathbf{x}) \beta \text{ KL}[q(\mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{c}|\mathbf{z})] d\mathbf{z}.$$

As in [Jiang et al., 2016], the first term does not depend on $\mathbf{c}$ and the second term is
non-negative. Thus, maximizing the lower bound ELBO with respect to $q(\mathbf{c}|\mathbf{x})$ requires that $\text{KL}[q(\mathbf{c}|\mathbf{x})||p_{\Theta}(\mathbf{c}|\mathbf{z})] = 0$. Thus, with $\nu$ a constant, we have:

$$ \frac{q(\mathbf{c}|\mathbf{x})}{p_{\Theta}(\mathbf{c}|\mathbf{z})} = \nu. $$

Since $\sum_{\mathbf{c}} q(\mathbf{c}|\mathbf{x}) = 1$ and $\sum_{\mathbf{c}} p_{\Theta}(\mathbf{c}|\mathbf{z}) = 1$, we have:

$$ \frac{q(\mathbf{c}|\mathbf{x})}{p_{\Theta}(\mathbf{c}|\mathbf{z})} = 1. $$

Taking the expectation on both sides, we can obtain:

$$ q(\mathbf{c}|\mathbf{x}) = \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[p_{\Theta}(\mathbf{c}|\mathbf{z})]. $$

# C Verification of GVI assumptions

In this section, we provide the detailed verification that the GCVAE model satisfies the assumptions required for Theorem 2 (Consistency under Independence) in [Knoblauch, 2019].

We denote the global parameters by $\psi = \{\Theta, \phi\} \in \Psi$, where $\Theta$ contains the generative parameters (GMM and decoder weights) and $\phi$ contains the variational parameters (encoder weights). We denote by $\mathcal{X} = \mathcal{X}^o \times \mathcal{Z}$, where $\mathcal{X}^o$ and $\mathcal{Z}$ denote the spaces of the observables $(\mathbf{x}_i, \mathbf{y}_i)$ and the latent components $(\mathbf{z}_i, \mathbf{c}_i)$.

**Assumption 1.** The GVI problem is well-defined.

1. The loss function $\mathcal{L}: \Psi \times \mathcal{X} \to \mathbb{R}$ is discontinuous at most at finitely many points. By using the Gumbel-Softmax relaxation, the loss function is a composition of continuous functions and is therefore continuous everywhere (0 points of discontinuity).

2. For any $(\mathbf{x}_i, \mathbf{y}_i) \in \mathcal{X}^o$ and any $n$, $\mathcal{L}(\psi, \mathbf{x}_i, \mathbf{y}_i, \mathbf{z}_i, \mathbf{c}_i) < \infty$ for all $(\mathbf{z}_i, \mathbf{c}_i) \in \mathcal{Z}$. Since the log-densities of the Gaussian and Concrete distributions are finite everywhere
on their respective open supports (assuming $\sigma > 0$, $\pi > 0$ and $c_k > 0$), the pointwise
loss value is finite for any generated sample.

3. The minimizers $\hat{\psi}_n = \arg \min_{\psi} \left\{ \frac{1}{n} \sum_{i=1}^n \mathcal{L}(\psi, \mathbf{x}_i, \mathbf{y}_i, \mathbf{z}_i, \mathbf{c}_i) \right\} \in \Psi \text{ exist for all } n.$

Since the loss function is continuous and defined on a compact set $\Psi$ (by hypothesis),
the Weierstrass extreme value theorem guarantees that a global minimizer exists.

For a probability measure $\mu$ on $\mathcal{X}$:

4. The loss satisfies a law of large numbers, i.e. $\frac{1}{n} \sum_{i=1}^n \mathcal{L}(\psi, \mathbf{x}_i, \mathbf{y}_i, \mathbf{z}_i, \mathbf{c}_i) \xrightarrow{\mu \text{ a.s.}} \mathbb{E}_\mu[\mathcal{L}(\psi, \mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c})]$.

We make the standard statistical learning assumption that the training observations
$\{(\mathbf{x}_i, \mathbf{y}_i)\}_{i=1}^n$ are i.i.d. Under the i.i.d. assumption, the law of large numbers applies.

5. The $\mu$-population-minimizer $\psi^* = \arg \min_{\psi} \mathbb{E}_{\mu}[\mathcal{L}(\psi, \mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c})] \in \Psi \text{ exists and is unique.}$

* Existence: As established in Assumption 1.1, the pointwise loss function is continuous. Under the integrability condition (Assumption 1.2), this implies the population risk is continuous with respect to $\psi$. Since the parameter space $\Psi$ is compact, the Weierstrass extreme value theorem guarantees that it attains its minimum. Thus, a minimizer exists.

* Uniqueness: Neural networks are inherently non-identifiable due to permutation symmetries. However, we make the standard identifiability assumption that the global minimum is unique up to these permutation symmetries. That is, we assume that any two distinct parameter vectors achieving the minimum are functionally equivalent and differ only by a permutation of their indices.

6. The loss is finite in $\mu$-expectation, i.e. $\mathbb{E}_{\mu}[\mathcal{L}(\psi, \mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c})] < \infty$ for all $\psi \in \Psi$.

We previously established that the loss is finite pointwise. The integration over the
latent space $\mathcal{Z}$ is well-defined because the Gaussian variational distributions have
finite moments. Furthermore, since the neural networks are Lipschitz continuous on
the compact parameter space $\Psi$, the loss grows at most quadratically with respect
to the observations $(\mathbf{x}, \mathbf{y})$. Therefore, under the standard assumption that the data
distribution has finite second moments, the total expectation is finite.

7. One of the following holds true: $\mathbb{E}_{\mu}[\mathcal{L}(\psi, \mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{c})]$ is coercive in $\psi$ or $\Psi$ is compact.
We define $\Psi$ as a compact set.

**Assumption 2** The variational family $Q = Q^{\psi} \times Q^{z}$ with $Q^{\psi} = \{q(\psi | \kappa) : \kappa \in \mathcal{K}\}$ and $Q^{z} = \{q^{z}(\mathbf{z}, \mathbf{c} | \eta) : \eta \in \mathcal{H}\}$ consists of absolutely continuous densities with respect to the Lebesgue measure. Moreover, for all $n \in \mathbb{N}$ and any $(\psi^{*}, \mathbf{z}, \mathbf{c}) \in \Psi \times \mathcal{Z}$, there exist sequences $\{\kappa_k\}_{k=1}^{\infty}$ and $\{\eta_k\}_{k=1}^{\infty}$ of variational parameters so that $q(\psi | \kappa_k) \xrightarrow{D} \delta_{\psi^{*}}(\psi)$ and $q^{z}(\mathbf{z}, \mathbf{c} | \eta_k) \xrightarrow{D} \delta_z(\mathbf{z}, \mathbf{c})$ as $k \to \infty$.

The proposed variational family satisfies these regularity conditions by design. First,
regarding absolute continuity, the employment of the Gumbel-Softmax relaxation endows
the latent distribution with a valid, absolutely continuous density with respect to the
Lebesgue measure on the simplex. The weak convergence to the Dirac measure is achieved
in the distributional limit as the temperature $\tau \to 0$ and variance $\sigma^2 \to 0$.

**Assumption 3** The GVI uncertainty quantifier $D: \mathcal{P}(\Psi)^2 \to \mathbb{R}_+$ is a statistical divergence. Further, it is lower semi-continuous in its first argument with respect to the weak topology of $\mathcal{P}(\Psi)$.

We employ the Kullback-Leibler (KL) divergence as the uncertainty quantifier. The KL divergence is a well-defined statistical divergence and is rigorously known to be lower semi-continuous. Thus, the assumption is satisfied.
**Assumption 4** The prior $\pi$ and the GVI uncertainty quantifier $D$ are suitable for the variational family $Q^\psi$: For all $q \in Q^\psi$, $D(q\|\pi) < \infty$.

The prior is modeled as a Gaussian Mixture Model, and the variational family $Q$ is a mean-field product of Gaussian and Gumbel-Softmax distributions. Since the variational components share the same support as the prior and possess finite moments, the KL divergence is well-defined and finite for all valid variational parameters.

**Assumption 5** The prior belief $\pi_{prior}$ about $\psi$ is not infinitely bad: $E_{\pi_{prior}}[\mathbb{E}_{\mu}[\mathcal{L}(\psi, x, y, z, c)]] = C_{\pi_{prior}} < \infty$. Moreover, $Q^\psi$ contains the singleton $\pi_{prior}(\psi)$. In other words, $Q^\psi = \{q(\psi | \kappa) : \kappa \in \mathbf{K}\} \cup \{\pi_{prior}(\psi)\}$.

The prior $\pi_{prior}(\psi)$ is modeled as a standard Gaussian. Since the loss function is continuous and grows polynomially with respect to the weights (Lipschitz neural networks), and the Gaussian prior has finite moments of all orders, the expected loss under the prior is strictly finite. The second condition is satisfied by construction.

**Assumption 7** There exists a compact subset $A \subset \Psi$ so that (i) $\psi^* \in A$ and (ii) $\pi_{prior} \ge \bar{q}_n$ on $\Psi \setminus A$, for all $n \ge N$ for some $N < \infty$.

Since we explicitly defined the parameter space $\Psi$ to be compact (Assumption 1.7), we can choose the subset $A$ to be the entire space, i.e., $A = \Psi$.

(i) The minimizer $\psi^*$ must reside within the valid parameter space $\Psi$.

(ii) The complement set $\Psi \setminus A$ is the empty set.

Therefore our model verifies the assumption of the theorem:

**Theorem C.1 (GVI consistency under independence).** If Assumptions 1, 2, 3, 4, 5 and 7 hold and $(x_i, y_i) \xrightarrow{iid} (x_1, y_1)$, then the GVI posteriors are consistent. i.e., $q_n \xrightarrow{D} \delta_{\psi^*}$ $\mu$-almost surely, where $\mu$ is the probability measure on $(x_1, y_1)$.
# D Withings dataset

The dataset is composed of 50,000 individuals, each contributing one night of sleep data to avoid bias from repeated measurements for the same user. The data was recorded by the Withings Sleep Analyzer<fnref n="2" />, capturing comprehensive sleep and biometric information. The dataset includes an equal number of users across the three categories based on the Apnea-Hypopnea Index (AHI) as categorized by the American Academy of Sleep Medicine (AASM) [Kirk et al., 2017]: no or mild (<15 events/hour), moderate (15–30 events/hour), and severe (>30 events/hour).

## D.1 Variables overview

The seven variables in the input vector $x$ are:

* sleep_duration: Duration of the user's sleep during the night.

* light_sleep_duration: Duration of the user's light sleep during the night.

* deep_sleep_duration: Duration of the user's deep sleep during the night.

* nb_sleep_interruptions: Count of awakenings throughout the night.

* avg_night_hr: Mean heart rate during the night.

* bmi: User's average Body Mass Index measured over a year.

* age: User's age.

The guiding variable $y$ is the *apnea_hypopnea_index* (AHI), categorizing sleep apnea severity based on the number of apnea-hypopnea events per hour.

<footnote n="2">²https://www.withings.com/us/en/sleep</footnote>
## D.2 Dataset Summary

In summary, Table 2 shows a concise view of each variable's range, mean, and standard deviation.

<table><caption>Table 2: Descriptive statistics of the variables</caption><thead><tr><td>Variable (unit)</td><td>Range</td><td>Mean</td><td>Std. Dev.</td></tr></thead><tbody><tr><td>sleep_duration (seconds)</td><td>14880 – 36000</td><td>26224</td><td>4224</td></tr><tr><td>light_sleep_duration (seconds)</td><td>3600 – 31860</td><td>15747</td><td>4651</td></tr><tr><td>deep_sleep_duration (seconds)</td><td>3600 – 32220</td><td>10472</td><td>4022</td></tr><tr><td>nb_sleep_interruptions</td><td>0 – 20</td><td>2.74</td><td>2.34</td></tr><tr><td>avg_night_hr (bpm)</td><td>40 – 111</td><td>62.49</td><td>8.57</td></tr><tr><td>bmi (kg/m²)</td><td>16 – 50</td><td>27.53</td><td>5.14</td></tr><tr><td>age (years)</td><td>18 – 80</td><td>50</td><td>12.67</td></tr><tr><td>apnea_hypopnea_index</td><td>0 – 40</td><td>18.14</td><td>13.47</td></tr></tbody></table>

We applied Min-Max normalization to all the variables to ensure they fall on a comparable scale.

## E Clusters profile tables for the Withings case study

This section provides the detailed cluster profiles for the user subgroups discovered by the unguided baseline and our proposed Guided Model (GCVAE), as discussed in Section 5.3. All values for user features are denormalized and presented as Mean ± Standard Deviation.
<table><caption>Table 3: Cluster profiles for the unguided model.</caption><thead><tr><th></th><th>Cluster 1</th><th>Cluster 2</th><th>Cluster 3</th></tr></thead><tbody><tr><td>AHI</td><td>9.43 ± 6.92</td><td>28.55 ± 7.50</td><td>8.55 ± 6.69</td></tr><tr><td>Age (years)</td><td>49.61 ± 11.52</td><td>55.91 ± 11.78</td><td>49.36 ± 12.48</td></tr><tr><td>BMI (kg/m²)</td><td>27.87 ± 5.08</td><td>29.25 ± 5.16</td><td>26.93 ± 4.84</td></tr><tr><td>Sleep Duration (hrs)</td><td>6.63 ± 0.81</td><td>7.30 ± 1.18</td><td>8.28 ± 0.83</td></tr><tr><td>Deep Sleep (hrs)</td><td>2.94 ± 1.01</td><td>2.91 ± 1.16</td><td>3.32 ± 1.15</td></tr><tr><td>Light Sleep (hrs)</td><td>3.69 ± 1.00</td><td>4.39 ± 1.27</td><td>4.96 ± 1.19</td></tr><tr><td>Sleep Interruptions</td><td>2.28 ± 2.09</td><td>3.04 ± 2.50</td><td>2.75 ± 2.25</td></tr><tr><td>Avg. Night HR (bpm)</td><td>63.66 ± 8.59</td><td>63.99 ± 8.50</td><td>61.98 ± 8.38</td></tr><tr><td>Avg. Resp. Rate (rpm)</td><td>15.11 ± 2.14</td><td>15.11 ± 2.26</td><td>15.06 ± 2.10</td></tr></tbody></table>

<table><caption>Table 4: Clinical profiles of user subgroups discovered by GCVAE.</caption><thead><tr><th></th><th>Cluster 1</th><th>Cluster 2</th><th>Cluster 3</th></tr></thead><tbody><tr><td>AHI</td><td>10.4 ± 9.6</td><td>16.4 ± 10.8</td><td>21.6 ± 11.2</td></tr><tr><td>Age (years)</td><td>42.58 ± 9.88</td><td>53.35 ± 10.39</td><td>61.90 ± 9.73</td></tr><tr><td>BMI (kg/m²)</td><td>25.35 ± 3.67</td><td>28.25 ± 4.57</td><td>31.79 ± 5.78</td></tr><tr><td>Sleep Duration (hrs)</td><td>7.52 ± 1.12</td><td>7.36 ± 1.18</td><td>7.19 ± 1.19</td></tr><tr><td>Deep Sleep (hrs)</td><td>3.41 ± 1.11</td><td>3.04 ± 1.09</td><td>2.50 ± 0.98</td></tr><tr><td>Light Sleep (hrs)</td><td>4.10 ± 1.23</td><td>4.32 ± 1.27</td><td>4.68 ± 1.25</td></tr><tr><td>Sleep Interruptions</td><td>2.03 ± 1.89</td><td>2.64 ± 2.17</td><td>3.83 ± 2.74</td></tr><tr><td>Avg. Night HR (bpm)</td><td>61.49 ± 8.08</td><td>63.66 ± 8.54</td><td>65.01 ± 8.77</td></tr><tr><td>Avg. Resp. Rate (rpm)</td><td>15.08 ± 2.05</td><td>15.09 ± 2.19</td><td>15.15 ± 2.29</td></tr></tbody></table>