Bayesian Multi-View Graph Convolutional Network (BMGCN) for Integrative Multi-Omics Analysis with Survival Outcomes and Zero-Inflated Data

1. Introduction

The ability to measure multiple molecular layers—such as gene expression (RNA-seq), DNA methylation, microRNA, and proteomics—from the same subjects has revolutionized systems biology and precision medicine. Each of these “omics” views provides a unique window into the biological state of a system. The central hypothesis of multi-omics integration is that a unified analysis of these views can reveal emergent biological insights and yield more accurate predictive models than any single view alone [1].

However, the integrative analysis of such data is fraught with challenges. First, the data are heterogeneous, with different views having unique statistical properties, dimensionalities, and levels of noise. Second, many omics datasets, particularly those from single-cell and sequencing-based assays, are characterized by a high degree of sparsity and zero-inflation, where an excess of zero values complicates modeling [2]. These zeros can represent true biological absence or technical artifacts (e.g., “dropouts”), and distinguishing between them is crucial. Third, predicting clinical endpoints, such as patient survival, requires specialized statistical models that can properly handle right-censored data, where the event of interest is not observed for all subjects [3]. Finally, ignoring the rich web of prior biological knowledge, such as protein-protein interaction networks or gene regulatory pathways, is a missed opportunity to guide the model toward more robust and interpretable solutions.

To address these multifaceted challenges, we introduce the Bayesian Multi-view Graph Convolutional Network (BMGCN). BMGCN is a deep generative framework that provides a principled solution for integrating multi-omics data with censored survival outcomes. The design of BMGCN is guided by the following core objectives:

• Latent Abstraction: To learn a low-dimensional representation of the data by factorizing the signal into latent variables that are shared across all omics views (${\mathbf{Z}}_i$) and variables that are specific to each view (${\mathbf{Z}}_i^{\left(k\right)}$). This captures both common and view-specific biological processes.
• View Selection: To automatically determine the relevance of each omics view for predicting the survival outcome. We achieve this through a Bayesian spike-and-slab prior on the view-specific latents (${\gamma }_k$), which can effectively “turn off” uninformative views, enhancing interpretability [4].
• Principled Sparsity Modeling: To explicitly model the excess zeros common in omics data using a zero-inflated mixture likelihood [2]. This distinguishes between true biological absence (a point mass at zero) and low-level biological signal (a continuous component).
• Censored Outcome Modeling: To directly model time-to-event data through an integrated Cox proportional hazards module that operates on the learned latent representations [3].
• Integration of Biological Priors: To leverage prior knowledge in the form of biological networks (${\mathbf{A}}^{\left(k\right)}$) by using Graph Convolutional Network (GCN) encoders, which regularize the model and encourage biologically meaningful latent representations [5].

This paper details the statistical foundation of BMGCN. We first present the notation and model overview, followed by a detailed description of the full generative process. We then derive the variational inference framework and the Evidence Lower Bound (ELBO) objective function. Finally, we describe the architecture of the graph-based encoders and outline a rigorous strategy for model validation and diagnostics.

2. Methods

2.1. Notation and Core Entities

Let our dataset consist of N subjects and K omics views. The core entities of BMGCN are defined as follows:

Observed Data:

• ${\mathbf{X}}_i^{\left(k\right)}\in {\mathbb{R}}^{p_k}$: The feature vector for subject i in view k, where $p_k$ is the number of features in that view.
• ${\mathbf{A}}^{\left(k\right)}\in {\{0,1\}}^{p_k\times p_k}$: An optional adjacency matrix for view k, representing prior biological knowledge (e.g., a protein-protein interaction network).
• $(T_i,{\Delta }_i)$: The observed time for subject i, where $T_i\in {\mathbb{R}}^{+}$ is the survival or censoring time and ${\Delta }_i\in \{0,1\}$ is the event indicator (${\Delta }_i=1$ if the event was observed, ${\Delta }_i=0$ if censored).

Latent Variables:

• ${\mathbf{Z}}_i\in {\mathbb{R}}^L$: A shared latent vector for subject i, capturing variation common to all views.
• ${\mathbf{Z}}_i^{\left(k\right)}\in {\mathbb{R}}^{L_k}$: A view-specific latent vector for subject i in view k.
• ${\gamma }_k\in \{0,1\}$: A binary variable indicating whether view k is active in predicting survival.
• ${\mathit{\pi }}_k\in (0,1)$: The prior inclusion probability for view k.

2.2. The BMGCN Generative Model

BMGCN defines a full generative process for the observed data, which factorizes into several key steps. For each subject $i=1,\dots ,N$:

• Sample View Inclusion Probabilities and Indicators: For each view $k=1,\dots ,K$, we first draw a prior inclusion probability from a Beta distribution and then sample the binary indicator variable:

  $$\begin{array}{cc}\hfill {\mathit{\pi }}_k& \sim \mathrm{Beta}({\alpha }_0,{\beta }_0)\hfill \end{array}$$

  (1)

  $$\begin{array}{cc}\hfill {\gamma }_k\mid {\mathit{\pi }}_k& \sim \mathrm{Bernoulli}\left({\mathit{\pi }}_k\right)\hfill \end{array}$$

  (2)
• Sample Latent Variables: A shared latent vector is drawn from a standard normal prior. The view-specific latent vectors are drawn from a spike-and-slab prior, conditioned on the inclusion variable ${\gamma }_k$. If ${\gamma }_k=0$, the latent vector is fixed at zero (the “spike”); otherwise, it is drawn from a standard normal (the “slab”).

  $$\begin{array}{cc}\hfill {\mathbf{Z}}_i& \sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_L)\hfill \end{array}$$

  (3)

  $$\begin{array}{cc}\hfill p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)& =(1-{\gamma }_k){\delta }_{\mathbf{0}}\left({\mathbf{Z}}_i^{\left(k\right)}\right)+{\gamma }_k\mathcal{N}({\mathbf{Z}}_i^{\left(k\right)}\mid \mathbf{0},{\mathbf{I}}_{L_k})\hfill \end{array}$$

  (4)

  where ${\delta }_{\mathbf{0}}$ is the Dirac delta function at zero.
• Generate Omics Features via Decoders: The observed features for each view are generated from the combined shared and view-specific latent variables. A decoder network $f_{\psi }^{\left(k\right)}$ maps the concatenated latent vector $[{\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)}]$ to the parameters of a zero-inflated Gaussian distribution.

  $$\begin{array}{cc}\hfill ({\mathit{\mu }}_i^{\left(k\right)},{\mathit{\rho }}_i^{\left(k\right)})& =f_{\psi }^{\left(k\right)}\left([{\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)}]\right)\hfill \end{array}$$

  (5)

  $$\begin{array}{cc}\hfill p(X_{ij}^{\left(k\right)}\mid \dots )& ={\rho }_{ij}^{\left(k\right)}{\delta }_0\left(X_{ij}^{\left(k\right)}\right)+(1-{\rho }_{ij}^{\left(k\right)})\mathcal{N}(X_{ij}^{\left(k\right)}\mid {\mathit{\mu }}_{ij}^{\left(k\right)},{\mathit{\sigma }}_k^2)\hfill \end{array}$$

  (6)
  Here, ${\mathit{\mu }}_i^{\left(k\right)}$ is the mean vector, ${\rho }_i^{\left(k\right)}$ is the vector of zero-inflation probabilities (passed through a sigmoid function to ensure they are in $(0,1)$), and ${\mathit{\sigma }}_k^2$ is a view-specific variance parameter.
• Generate Survival Outcome: A latent risk score ${\eta }_i$ for each subject is computed by a survival head network $g_{\theta }$, which takes the shared latent and the gated view-specific latents as input. The observed survival data $(T_i,{\Delta }_i)$ are assumed to follow a Cox proportional hazards model conditioned on this risk score.

  $$\begin{array}{cc}\hfill {\eta }_i& =g_{\theta }\left([{\mathbf{Z}}_i,{\gamma }_1{\mathbf{Z}}_i^{\left(1\right)},\dots ,{\gamma }_K{\mathbf{Z}}_i^{\left(K\right)}]\right)\hfill \end{array}$$

  (7)

  $$\begin{array}{cc}\hfill (T_i,{\Delta }_i)& \sim \mathrm{Cox}\left({\eta }_i\right)\hfill \end{array}$$

  (8)
  The likelihood for the Cox model is given by the partial log-likelihood:

  $$\begin{array}{c}\hfill l\left(\eta \right)=\sum _{i=1}^N{\Delta }_i\left({\eta }_i-log\sum _{j\in R_i}{e}^{{\eta }_j}\right)\end{array}$$

  (9)

  where $R_i=\{j:T_j\ge T_i\}$ is the set of subjects at risk at time $T_i$. This semi-parametric form is powerful because it does not require specifying a baseline hazard function $h_0\left(t\right)$. When multiple events occur at the same time (ties), approximations such as Breslow’s or Efron’s can be used.
  The full joint distribution over all observed and latent variables is given by:

  $$\begin{array}{c}p(\mathbf{X},\mathbf{T},\Delta ,\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })=\hfill \\ \hfill \left[\prod _{i=1}^{N}p(T_i,{\Delta }_i\mid {\eta }_i)\prod _{k=1}^{K}p({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)p\left({\mathbf{Z}}_i\right)\right]\times \hfill \\ \hfill \left[\prod _{k=1}^{K}p({\gamma }_k\mid {\mathit{\pi }}_k)p\left({\mathit{\pi }}_k\right)\right]\end{array}$$

  (10)

2.3. Variational Inference and the ELBO

The posterior distribution over the latent variables, $p(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma \mid \mathbf{X},\mathbf{T},\mathbf{\Delta })$, is intractable due to the non-linearities in the model and the Cox likelihood. We therefore resort to amortized variational inference (VI) [6,7]. We introduce a tractable variational distribution $q(·)$ to approximate the true posterior and maximize the Evidence Lower Bound (ELBO), $\mathcal{L}$, with respect to the variational parameters.

We use a mean-field variational family that factorizes across subjects and latent variables:

$$\begin{array}{c}\hfill q(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )=\prod _{i=1}^{N}q\left({\mathbf{Z}}_i\right)\prod _{k=1}^{K}q\left({\mathbf{Z}}_i^{\left(k\right)}\right)q\left({\gamma }_k\right)\end{array}$$

(11)

where each factor is parameterized by an encoder network:

$$\begin{array}{cc}\hfill q({\mathbf{Z}}_i\mid {\mathbf{X}}_i)& =\mathcal{N}({\mathit{\mu }}_i\left({\mathbf{X}}_i\right),\mathrm{diag}\left({\mathit{\sigma }}_i^2\left({\mathbf{X}}_i\right)\right))\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill q({\mathbf{Z}}_i^{\left(k\right)}\mid {\mathbf{X}}_i^{\left(k\right)},{\mathbf{A}}^{\left(k\right)})& =\mathcal{N}({\mathit{\mu }}_i^{\left(k\right)}({\mathbf{X}}_i^{\left(k\right)},{\mathbf{A}}^{\left(k\right)}),\mathrm{diag}\left({\left({\mathit{\sigma }}_i^{\left(k\right)}\right)}^2({\mathbf{X}}_i^{\left(k\right)},{\mathbf{A}}^{\left(k\right)})\right))\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill q\left({\gamma }_k\right)& =\mathrm{Bernoulli}\left({\tau }_k\right)\hfill \end{array}$$

(14)

where ${\tau }_k$ is the posterior inclusion probability for view k.

The ELBO is defined as:

$$\begin{array}{c}\hfill \mathcal{L}={\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )]-{\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )]\end{array}$$

(15)

This objective can be decomposed into several intuitive terms:

$$\begin{array}{c}\hfill \mathcal{L}={\mathcal{L}}_{\mathrm{recon}}+{\mathcal{L}}_{\mathrm{surv}}-{\mathcal{L}}_{\mathrm{KL}}\end{array}$$

(16)

where:

• Reconstruction Term (${\mathcal{L}}_{\mathbf{recon}}$): The expected log-likelihood of the omics data.

  $$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{recon}}=\sum _{i=1}^N\sum _{k=1}^K{\mathbb{E}}_q[logp({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})]\end{array}$$

  (17)
  The piecewise nature of the zero-inflated log-likelihood makes this term directly computable without explicit sampling of discrete zero-inflation masks.
• Survival Term (${\mathcal{L}}_{\mathbf{surv}}$): The expected Cox partial log-likelihood.

  $$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{surv}}={\mathbb{E}}_q\left[\sum _{i=1}^N{\mathbf{\Delta }}_i\left({\eta }_i-log\sum _{j\in R_i}{e}^{{\eta }_j}\right)\right]\end{array}$$

  (18)
• KL Divergence Regularizer (${\mathcal{L}}_{\mathbf{KL}}$): A sum of Kullback-Leibler (KL) divergences that penalize deviations of the variational posterior from the prior.

  $$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{KL}}& =\sum _{i=1}^N\mathrm{KL}(q\left({\mathbf{Z}}_i\right)\parallel p\left({\mathbf{Z}}_i\right))\hfill \end{array}$$

  (19)

  $$\begin{array}{cc}\hfill & +\sum _{i=1}^N\sum _{k=1}^K{\mathbb{E}}_{q\left({\gamma }_k\right)}\left[\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))\right]\hfill \end{array}$$

  (20)

  $$\begin{array}{cc}\hfill & +\sum _{k=1}^K\mathrm{KL}(q\left({\gamma }_k\right)\parallel p\left({\gamma }_k\right))\hfill \end{array}$$

  (21)
  The KL terms for the Gaussian latents have a closed-form solution. The term involving the discrete ${\gamma }_k$ is approximated as ${\tau }_k·\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel \mathcal{N}(\mathbf{0},\mathbf{I}))$, effectively weighting the KL penalty by the posterior probability that the view is included.

2.4. Graph-Convolutional Encoders

To incorporate prior biological network information (${\mathbf{A}}^{\left(k\right)}$), the encoder networks that parameterize $q\left({\mathbf{Z}}_i^{\left(k\right)}\right)$ are implemented as Graph Convolutional Networks (GCNs) [5]. This allows the model to learn feature representations that are regularized by known biological interactions, improving robustness and interpretability.

For each view k, we first normalize the adjacency matrix to ensure stable learning:

$$\begin{array}{c}\hfill {\widehat{\mathbf{A}}}^{\left(k\right)}={\left({\tilde{\mathbf{D}}}^{\left(k\right)}\right)}^{-1/2}{\tilde{\mathbf{A}}}^{\left(k\right)}{\left({\tilde{\mathbf{D}}}^{\left(k\right)}\right)}^{-1/2}\end{array}$$

(22)

where ${\tilde{\mathbf{A}}}^{\left(k\right)}={\mathbf{A}}^{\left(k\right)}+\mathbf{I}$ (adding self-loops) and ${\tilde{\mathbf{D}}}^{\left(k\right)}$ is the diagonal degree matrix of ${\tilde{\mathbf{A}}}^{\left(k\right)}$.

The GCN takes the full data matrix for view k, ${\mathbf{X}}^{\left(k\right)}\in {\mathbb{R}}^{N\times p_k}$, as input features for the nodes (e.g., genes) of the graph. The GCN then learns feature embeddings by propagating and transforming information between nodes in the graph. The propagation rule for layer l is:

$$\begin{array}{c}\hfill {\mathbf{H}}^{(l+1)}=\mathrm{ReLU}\left({\widehat{\mathbf{A}}}^{\left(k\right)}{\mathbf{H}}^{\left(l\right)}{\mathbf{W}}^{\left(l\right)}\right)\end{array}$$

(23)

where ${\mathbf{H}}^{\left(0\right)}={\left({\mathbf{X}}^{\left(k\right)}\right)}^{\top }\in {\mathbb{R}}^{p_k\times N}$, ${\mathbf{W}}^{\left(l\right)}$ is a learnable weight matrix, and ${\mathbf{H}}^{\left(l\right)}\in {\mathbb{R}}^{p_k\times d_l}$ is the feature embedding matrix at layer l. After L layers, we obtain a final feature embedding matrix ${\mathbf{H}}^{\left(L\right)}\in {\mathbb{R}}^{p_k\times d_L}$.

To obtain sample-specific parameters for $q\left({\mathbf{Z}}_i^{\left(k\right)}\right)$, we apply a readout function. Specifically, the final feature embeddings are used to parameterize a sample-specific network that maps the input data to the latent space parameters:

$$\begin{array}{c}\hfill {\mathbf{h}}_i^{\left(k\right)}={\mathrm{FFN}}_{\mathrm{enc}}({\mathbf{X}}_i^{\left(k\right)},{\mathbf{H}}^{\left(L\right)})\end{array}$$

(24)

This sample-specific embedding ${\mathbf{h}}_i^{\left(k\right)}$ is then mapped to the variational parameters:

$$\begin{array}{cc}\hfill {\mathit{\mu }}_i^{\left(k\right)}& ={\mathbf{W}}_{\mathit{\mu }}{\mathbf{h}}_i^{\left(k\right)}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill log{\mathit{\sigma }}_i^{\left(k\right)}& ={\mathbf{W}}_{\mathit{\sigma }}{\mathbf{h}}_i^{\left(k\right)}\hfill \end{array}$$

(26)

2.5. Optimization and Implementation

The ELBO is optimized with respect to the parameters of the encoders, decoders, and survival head using stochastic gradient ascent. We employ the reparameterization trick for Gaussian latents to obtain low-variance gradient estimates [6]. For the discrete view-selection variables ${\gamma }_k$, we use a hard-concrete relaxation with a straight-through gradient estimator during training [8,9].

Training is performed in mini-batches. For the Cox survival term, which is defined over the full cohort, we use mini-batch approximations, such as constructing risk sets only from subjects within the current batch or using a memory bank. While this introduces a small bias, it is essential for scalability.

Key implementation details for ensuring stable training include:

• Using the log-sum-exp trick for stable computation of the Cox likelihood denominator.
• Applying KL annealing to gradually introduce the KL regularization term, preventing posterior collapse [10].
• Using modern optimizers like AdamW and employing learning rate schedulers [11,12].
• Initializing the view-selection logits (${\tau }_k$) to favor inclusion initially to promote stable learning.

3. Validation and Diagnostics Strategy

To ensure BMGCN is robust, generalizable, and interpretable, we propose a multi-faceted validation strategy.

• Simulation Studies: We will generate synthetic multi-view data with known ground-truth latent structures, pre-defined informative vs. uninformative views, and known survival dependencies. This will allow us to quantitatively assess the model’s ability to recover the true latent variables (cosine similarity), correctly perform view selection (accuracy of ${\tau }_k$ vs. ground truth), and accurately predict survival (Concordance Index).
• Cross-Validation: On real-world datasets, we will use a k-fold cross-validation scheme, stratified by event status, to evaluate predictive performance on held-out data. The primary metric will be the Concordance Index (C-index), which measures the fraction of concordant pairs of subjects. A secondary metric will be the integrated Brier score, which assesses the calibration of survival probability predictions over time.
• Posterior Predictive Checks (PPCs): To assess goodness-of-fit, we will sample from the posterior predictive distribution. We will compare the distributions of simulated data (e.g., zero frequencies, means, variances) against the observed data to detect model misspecification.
• Cox Diagnostics: We will examine Schoenfeld residuals to test the proportional hazards assumption of the Cox model, a critical assumption for the validity of the survival module.
• Ablation Studies: To quantify the contribution of each model component, we will perform ablation studies by systematically removing key features: (1) the graph structure (setting ${\mathbf{A}}^{\left(k\right)}=\mathbf{I}$), (2) the zero-inflation module ($\rho =\mathbf{0}$), and (3) the spike-and-slab selectors (forcing ${\gamma }_k=1$ for all views).
• Latent Space Analysis: We will visualize the learned shared latent space ${\mathbf{Z}}_i$ using techniques like UMAP to check if it clusters subjects by known phenotypes, treatment groups, or survival outcomes, providing a qualitative check on its biological relevance [13].

4. Discussion

BMGCN offers a comprehensive and principled framework for one of the most pressing problems in computational biology: the integrative analysis of multi-omics data for clinical outcome prediction. Its primary strengths lie in its ability to simultaneously address multiple, often co-occurring, challenges.

The Bayesian framework provides not only point estimates but also uncertainty quantification. The probabilistic view selection mechanism offers a path to identifying which data modalities are most informative for a given prediction task, a crucial step for biomarker discovery and for designing future cost-effective studies.

However, the model is not without limitations. The mean-field approximation used for variational inference may not capture complex posterior dependencies between latent variables. The Cox model component relies on the proportional hazards assumption, which may not hold in all biological contexts; diagnostics are essential. Computationally, BMGCN is intensive, particularly with a large number of views or features, though scalability is addressed via mini-batching.

Future work could extend BMGCN in several directions. More flexible variational families, such as normalizing flows, could be explored to improve posterior approximation. The model could be extended to handle other data types (e.g., categorical or count data) and other clinical endpoints. Furthermore, incorporating mechanisms to explicitly handle missing views for certain subjects would greatly enhance its applicability to real-world clinical datasets.

5. Conclusions

We have presented BMGCN, a Bayesian deep generative model for the integrative analysis of multi-omics data with censored survival outcomes. By combining shared and view-specific latent variables, graph-convolutional encoders, a zero-inflated likelihood, and a Cox survival head within a unified probabilistic framework, BMGCN is tailored to the unique challenges of modern biomedical data. Its design prioritizes flexibility, interpretability, and statistical rigor, providing a powerful new tool for uncovering the complex molecular underpinnings of disease and for improving patient outcome prediction.

Appendix A.1. Derivation of the Evidence Lower Bound (ELBO)

The goal of variational inference is to approximate the true, intractable posterior distribution $p(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi }\mid \mathbf{X},\mathbf{T},\mathbf{\Delta })$ with a tractable variational distribution $q(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })$. We assume a mean-field factorization for q as specified in the main text:

$$\begin{array}{c}\hfill q(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })=q\left(\mathbf{Z}\right)q\left(\left\{{\mathbf{Z}}^{\left(k\right)}\right\}\right)q\left(\gamma \right)q\left(\mathit{\pi }\right)=\left[\prod _{i=1}^{N}q\left({\mathbf{Z}}_i\right)\right]\left[\prod _{i=1}^N\prod _{k=1}^{K}q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right]\left[\prod _{k=1}^{K}q\left({\gamma }_k\right)\right]\left[\prod _{k=1}^{K}q\left({\mathit{\pi }}_k\right)\right]\end{array}$$

(A1)

The quality of this approximation is measured by the Kullback-Leibler (KL) divergence between q and p:

$$\begin{array}{c}\hfill \mathrm{KL}(q\parallel p)={\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]-{\mathbb{E}}_q[logp(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi }\mid \mathbf{X},\mathbf{T},\mathbf{\Delta })]\end{array}$$

(A2)

By Bayes’ theorem, the log-posterior is:

$$\begin{array}{c}\hfill logp(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi }\mid \mathbf{X},\mathbf{T},\mathbf{\Delta })=logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })-logp(\mathbf{X},\mathbf{T},\mathbf{\Delta })\end{array}$$

(A3)

Substituting this into the KL divergence:

$$\begin{array}{c}\hfill \mathrm{KL}(q\parallel p)={\mathbb{E}}_q[logq]-{\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]+logp(\mathbf{X},\mathbf{T},\mathbf{\Delta })\end{array}$$

(A4)

Rearranging terms, we get:

$$\begin{array}{c}\hfill logp(\mathbf{X},\mathbf{T},\mathbf{\Delta })={\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]-{\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]+\mathrm{KL}(q\parallel p)\end{array}$$

(A5)

Since $\mathrm{KL}(q\parallel p)\ge 0$, the first two terms on the right-hand side form a lower bound on the log-evidence (marginal likelihood) $logp(\mathbf{X},\mathbf{T},\mathbf{\Delta })$. This is the Evidence Lower Bound (ELBO), denoted by $\mathcal{L}$:

$$\begin{array}{c}\hfill \mathcal{L}={\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]-{\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]\end{array}$$

(A6)

This is Equation (15) from the main text. The ELBO is maximized with respect to the parameters of q to find the best approximation to the true posterior.

Appendix A.2. Decomposition of the ELBO

We now decompose the ELBO $\mathcal{L}$ into its constituent parts: the reconstruction term (${\mathcal{L}}_{\mathrm{recon}}$), the survival term (${\mathcal{L}}_{\mathrm{surv}}$), and the KL divergence regularizer (${\mathcal{L}}_{\mathrm{KL}}$).

Starting from the joint distribution defined in Equation (10):

$$\begin{array}{c}p(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })=\hfill \\ \left[\prod _{i=1}^{N}p(T_i,{\mathbf{\Delta }}_i\mid {\eta }_i)\prod _{k=1}^{K}p({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)p\left({\mathbf{Z}}_i\right)\right]\times \\ \hfill \left[\prod _{k=1}^{K}p({\gamma }_k\mid {\mathit{\pi }}_k)p\left({\mathit{\pi }}_k\right)\right]\end{array}$$

(A7)

Taking the expectation under q:

$$\begin{array}{cc}\hfill {\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]=& {\mathbb{E}}_q\left[\sum _{i=1}^{N}logp(T_i,{\mathbf{\Delta }}_i\mid {\eta }_i)\right]+{\mathbb{E}}_q\left[\sum _{i=1}^N\sum _{k=1}^{K}logp({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})\right]\hfill \\ \hfill & +{\mathbb{E}}_q\left[\sum _{i=1}^N\sum _{k=1}^{K}logp({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)\right]+{\mathbb{E}}_q\left[\sum _{i=1}^{N}logp\left({\mathbf{Z}}_i\right)\right]\hfill \\ \hfill & +{\mathbb{E}}_q\left[\sum _{k=1}^{K}logp({\gamma }_k\mid {\mathit{\pi }}_k)\right]+{\mathbb{E}}_q\left[\sum _{k=1}^{K}logp\left({\mathit{\pi }}_k\right)\right]\hfill \end{array}$$

(A8)

The second term in the ELBO is the entropy of the variational distribution:

$$\begin{array}{cc}\hfill {\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma ,\mathit{\pi })]=& {\mathbb{E}}_q\left[\sum _{i=1}^{N}logq\left({\mathbf{Z}}_i\right)\right]+{\mathbb{E}}_q\left[\sum _{i=1}^N\sum _{k=1}^{K}logq\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right]\hfill \\ \hfill & +{\mathbb{E}}_q\left[\sum _{k=1}^{K}logq\left({\gamma }_k\right)\right]+{\mathbb{E}}_q\left[\sum _{k=1}^{K}logq\left({\mathit{\pi }}_k\right)\right]\hfill \end{array}$$

(A9)

Substituting these into the ELBO $\mathcal{L}$:

$$\begin{array}{cc}\hfill \mathcal{L}=& \underset{{\mathcal{L}}_{\mathrm{recon}}}{\underbrace{{\mathbb{E}}_q\left[{\sum }_{i=1}^N{\sum }_{k=1}^{K}logp({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})\right]}}+\underset{{\mathcal{L}}_{\mathrm{surv}}}{\underbrace{{\mathbb{E}}_q\left[{\sum }_{i=1}^{N}logp(T_i,{\mathbf{\Delta }}_i\mid {\eta }_i)\right]}}\hfill \\ \hfill & -\underset{\mathrm{KL}\mathrm{for}{\mathbf{Z}}_i}{\underbrace{\left({\mathbb{E}}_q\left[{\sum }_{i=1}^{N}logq\left({\mathbf{Z}}_i\right)\right]-{\mathbb{E}}_q\left[{\sum }_{i=1}^{N}logp\left({\mathbf{Z}}_i\right)\right]\right)}}\hfill \\ \hfill & -\underset{\mathrm{KL}\mathrm{for}{\mathbf{Z}}_i^{\left(k\right)}}{\underbrace{\left({\mathbb{E}}_q\left[{\sum }_{i=1}^N{\sum }_{k=1}^{K}logq\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right]-{\mathbb{E}}_q\left[{\sum }_{i=1}^N{\sum }_{k=1}^{K}logp({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)\right]\right)}}\hfill \\ \hfill & -\underset{\mathrm{KL}\mathrm{for}{\gamma }_k}{\underbrace{\left({\mathbb{E}}_q\left[{\sum }_{k=1}^{K}logq\left({\gamma }_k\right)\right]-{\mathbb{E}}_q\left[{\sum }_{k=1}^{K}logp({\gamma }_k\mid {\mathit{\pi }}_k)\right]\right)}}\hfill \\ \hfill & -\underset{\mathrm{KL}\mathrm{for}{\mathit{\pi }}_k}{\underbrace{\left({\mathbb{E}}_q\left[{\sum }_{k=1}^{K}logq\left({\mathit{\pi }}_k\right)\right]-{\mathbb{E}}_q\left[{\sum }_{k=1}^{K}logp\left({\mathit{\pi }}_k\right)\right]\right)}}\hfill \end{array}$$

(A10)

This can be grouped as follows:

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{recon}}& ={\mathbb{E}}_q\left[\sum _{i=1}^N\sum _{k=1}^{K}logp({\mathbf{X}}_i^{\left(k\right)}\mid {\mathbf{Z}}_i,{\mathbf{Z}}_i^{\left(k\right)})\right](\mathrm{Reconstruction}\mathrm{Term})\hfill \end{array}$$

(A11)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{surv}}& ={\mathbb{E}}_q\left[\sum _{i=1}^{N}logp(T_i,{\mathbf{\Delta }}_i\mid {\eta }_i)\right]={\mathbb{E}}_q\left[\sum _{i=1}^N{\mathbf{\Delta }}_i\left({\eta }_i-log\sum _{j\in R_i}{e}^{{\eta }_j}\right)\right](\mathrm{Survival}\mathrm{Term})\hfill \\ \hfill {\mathcal{L}}_{\mathrm{KL}}& =\sum _{i=1}^N\mathrm{KL}(q\left({\mathbf{Z}}_i\right)\parallel p\left({\mathbf{Z}}_i\right))+\sum _{i=1}^N\sum _{k=1}^K{\mathbb{E}}_q\left[\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))\right]\hfill \end{array}$$

(A12)

$$\begin{array}{cc}\hfill & +\sum _{k=1}^K\mathrm{KL}(q\left({\gamma }_k\right)\parallel p({\gamma }_k\mid {\mathit{\pi }}_k))+\sum _{k=1}^K\mathrm{KL}(q\left({\mathit{\pi }}_k\right)\parallel p\left({\mathit{\pi }}_k\right))\hfill \end{array}$$

(A13)

This confirms the decomposition $\mathcal{L}={\mathcal{L}}_{\mathrm{recon}}+{\mathcal{L}}_{\mathrm{surv}}-{\mathcal{L}}_{\mathrm{KL}}$ from Equation (16). For simplicity in the main text, the priors on ${\gamma }_k$ and ${\mathit{\pi }}_k$ are often absorbed or their KL terms are handled implicitly, leading to the simplified ${\mathcal{L}}_{\mathrm{KL}}$ shown in Equations (19-21).

Appendix A.3. Zero-Inflated Gaussian Log-Likelihood and Gradients

The probability density function (PDF) for the zero-inflated Gaussian (ZIG) distribution for a single feature j in view k for subject i is:

$$\begin{array}{c}\hfill p(X_{ij}^{\left(k\right)}\mid {\mathit{\mu }}_{ij}^{\left(k\right)},{\mathit{\sigma }}_k^2,{\rho }_{ij}^{\left(k\right)})={\rho }_{ij}^{\left(k\right)}·{\delta }_0\left(X_{ij}^{\left(k\right)}\right)+(1-{\rho }_{ij}^{\left(k\right)})·\mathcal{N}(X_{ij}^{\left(k\right)}\mid {\mathit{\mu }}_{ij}^{\left(k\right)},{\mathit{\sigma }}_k^2)\end{array}$$

(A14)

where ${\delta }_0\left(x\right)$ is the Dirac delta function, which is infinite at $x=0$ and zero elsewhere, with the property that $\int {\delta }_0\left(x\right)dx=1$.

The log-likelihood for an observation x is piecewise, as it must handle the discrete mass at zero and the continuous density elsewhere:

$$\begin{array}{c}\hfill logp(x\mid \mathit{\mu },{\mathit{\sigma }}^2,\rho )=\left\{\begin{array}{cc}log\left(\rho +(1-\rho )·\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)\right)\hfill & \mathrm{if}x=0\hfill \\ log(1-\rho )+log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)\hfill & \mathrm{if}x\ne 0\hfill \end{array}\right.\end{array}$$

(A15)

To optimize the model, we need the gradients of this log-likelihood with respect to its parameters $\mathit{\mu }$, ${\mathit{\sigma }}^2$, and $\rho$.

Case 1: $x=0$

Let $C=\rho +(1-\rho )·\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)$.

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial \mathit{\mu }}}[logC]& ={\displaystyle \frac{1}{C}}·(1-\rho )·{\displaystyle \frac{\partial }{\partial \mathit{\mu }}}\left[\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)\right]\hfill \\ & ={\displaystyle \frac{1}{C}}·(1-\rho )·\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)·\left({\displaystyle \frac{(0-\mathit{\mu })}{{\mathit{\sigma }}^2}}\right)\hfill \end{array}$$

(A16)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}^2}}[logC]& ={\displaystyle \frac{1}{C}}·(1-\rho )·{\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}^2}}\left[\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{C}}·(1-\rho )·\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)·\left[{\displaystyle \frac{{(0-\mathit{\mu })}^2}{2{\mathit{\sigma }}^4}}-{\displaystyle \frac{1}{2{\mathit{\sigma }}^2}}\right]\hfill \end{array}$$

(A17)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial \rho }}[logC]& ={\displaystyle \frac{1}{C}}·\left[1-\mathcal{N}(0\mid \mathit{\mu },{\mathit{\sigma }}^2)\right]\hfill \end{array}$$

(A18)

Case 2: $x\ne 0$

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial \mathit{\mu }}}[log(1-\rho )+log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)]& ={\displaystyle \frac{\partial }{\partial \mathit{\mu }}}[log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)]={\displaystyle \frac{(x-\mathit{\mu })}{{\mathit{\sigma }}^2}}\hfill \end{array}$$

(A19)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}^2}}[log(1-\rho )+log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)]& ={\displaystyle \frac{\partial }{\partial {\mathit{\sigma }}^2}}[log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)]=\left[{\displaystyle \frac{{(x-\mathit{\mu })}^2}{2{\mathit{\sigma }}^4}}-{\displaystyle \frac{1}{2{\mathit{\sigma }}^2}}\right]\hfill \end{array}$$

(A20)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial \rho }}[log(1-\rho )+log\mathcal{N}(x\mid \mathit{\mu },{\mathit{\sigma }}^2)]& =-{\displaystyle \frac{1}{(1-\rho )}}\hfill \end{array}$$

(A21)

These gradients are used in the backpropagation algorithm to update the parameters of the decoder network $f_{\psi }^{\left(k\right)}$.

Appendix A.4. Gradient of the Cox Partial Log-Likelihood

The Cox partial log-likelihood for the entire dataset is:

$$\begin{array}{c}\hfill l\left(\eta \right)=\sum _{i=1}^N{\mathbf{\Delta }}_i\left({\eta }_i-log\sum _{j\in R_i}{e}^{{\eta }_j}\right)\end{array}$$

(A22)

where $R_i=\{j:T_j\ge T_i\}$ is the risk set for subject i.

We derive the gradient of $l\left(\eta \right)$ with respect to the risk score ${\eta }_m$ for a specific subject m.

The gradient $\partial l/\partial {\eta }_m$ receives contributions from two places in the summation:

• When $i=m$, if subject m experienced the event (${\mathbf{\Delta }}_m=1$), then ${\eta }_m$ appears directly in the first term.
• For any subject i for whom subject m is in the risk set (i.e., $m\in R_i$), ${\eta }_m$ appears in the log-sum-exp term $log\left({\sum }_{j\in R_i}{e}^{{\eta }_j}\right)$.

Therefore:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\eta }_m}}& ={\mathbf{\Delta }}_m·1-\sum _{i:m\in R_i}{\mathbf{\Delta }}_i·\left[{\displaystyle \frac{1}{{\sum }_{j\in R_i}{e}^{{\eta }_j}}}\right]·e^{{\eta }_m}\hfill \\ \hfill & ={\mathbf{\Delta }}_m-\sum _{i:m\in R_i}{\mathbf{\Delta }}_i·\left[{\displaystyle \frac{e^{{\eta }_m}}{{\sum }_{j\in R_i}{e}^{{\eta }_j}}}\right]\hfill \end{array}$$

(A23)

This gradient is used to update the parameters of the survival head network $g_{\theta }$.

Appendix A.5. KL Divergence for Gaussian Distributions

We derive the closed-form solution for the KL divergence between two multivariate Gaussian distributions, which is used for the latent variables ${\mathbf{Z}}_i$ and ${\mathbf{Z}}_i^{\left(k\right)}$.

Let $q\left(\mathbf{z}\right)=\mathcal{N}({\mathit{\mu }}_q,{\Sigma }_q)$ and $p\left(\mathbf{z}\right)=\mathcal{N}({\mathit{\mu }}_p,{\Sigma }_p)$. The KL divergence is:

$$\begin{array}{c}\hfill \mathrm{KL}(q\parallel p)={\mathbb{E}}_q[logq\left(\mathbf{z}\right)]-{\mathbb{E}}_q[logp\left(\mathbf{z}\right)]\end{array}$$

(A24)

The expectation of the log of a multivariate Gaussian $\mathcal{N}(\mathit{\mu },\Sigma )$ is:

$$\begin{array}{c}\hfill {\mathbb{E}}_q[log\mathcal{N}(\mathbf{z}\mid \mathit{\mu },\Sigma )]=-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}log\left|\Sigma \right|-{\displaystyle \frac{1}{2}}{\mathbb{E}}_q\left[{(\mathbf{z}-\mathit{\mu })}^{\top }{\Sigma }^{-1}(\mathbf{z}-\mathit{\mu })\right]\end{array}$$

(A25)

For $q\left(\mathbf{z}\right)$, this is:

$$\begin{array}{cc}\hfill {\mathbb{E}}_q[logq\left(\mathbf{z}\right)]& =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}log\left|{\Sigma }_q\right|-{\displaystyle \frac{1}{2}}{\mathbb{E}}_q\left[{(\mathbf{z}-{\mathit{\mu }}_q)}^{\top }{\Sigma }_q^{-1}(\mathbf{z}-{\mathit{\mu }}_q)\right]\hfill \\ \hfill & =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}log\left|{\Sigma }_q\right|-{\displaystyle \frac{1}{2}}\mathrm{trace}\left({\Sigma }_q^{-1}{\Sigma }_q\right)\hfill \\ \hfill & =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}log\left|{\Sigma }_q\right|-{\displaystyle \frac{d}{2}}\hfill \end{array}$$

(A26)

because ${\mathbb{E}}_q\left[{(\mathbf{z}-{\mathit{\mu }}_q)}^{\top }{\Sigma }_q^{-1}(\mathbf{z}-{\mathit{\mu }}_q)\right]=\mathrm{trace}\left({\Sigma }_q^{-1}{\mathbb{E}}_q\left[(\mathbf{z}-{\mathit{\mu }}_q){(\mathbf{z}-{\mathit{\mu }}_q)}^{\top }\right]\right)=\mathrm{trace}\left({\Sigma }_q^{-1}{\Sigma }_q\right)=\mathrm{trace}\left(\mathbf{I}\right)=d$.

For $p\left(\mathbf{z}\right)$, assuming a standard normal prior $p\left(\mathbf{z}\right)=\mathcal{N}(\mathbf{0},\mathbf{I})$ (so ${\mathit{\mu }}_p=\mathbf{0}$, ${\Sigma }_p=\mathbf{I}$), we have:

$$\begin{array}{cc}\hfill {\mathbb{E}}_q[logp\left(\mathbf{z}\right)]& =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}log\left|\mathbf{I}\right|-{\displaystyle \frac{1}{2}}{\mathbb{E}}_q\left[{\mathbf{z}}^{\top }\mathbf{z}\right]\hfill \\ \hfill & =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}{\mathbb{E}}_q\left[\sum _{l=1}^d{z}_l^2\right]\hfill \\ \hfill & =-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}\sum _{l=1}^d{\mathbb{E}}_q\left[z_l^2\right]\hfill \end{array}$$

(A27)

Since ${\mathbb{E}}_q\left[z_l^2\right]={\mathrm{Var}}_q\left(z_l\right)+{\left({\mathbb{E}}_q\left[z_l\right]\right)}^2={\mathit{\sigma }}_{q,l}^2+{\mathit{\mu }}_{q,l}^2$, we get:

$$\begin{array}{c}\hfill {\mathbb{E}}_q[logp\left(\mathbf{z}\right)]=-{\displaystyle \frac{d}{2}}log\left(2\mathit{\pi }\right)-{\displaystyle \frac{1}{2}}\sum _{l=1}^d({\mathit{\mu }}_{q,l}^2+{\mathit{\sigma }}_{q,l}^2)\end{array}$$

(A28)

Subtracting Equation (A28) from Equation (A26):

$$\begin{array}{cc}\hfill \mathrm{KL}(q\parallel p)& =\left[-{\displaystyle \frac{1}{2}}log\left|{\Sigma }_q\right|-{\displaystyle \frac{d}{2}}\right]-\left[-{\displaystyle \frac{1}{2}}\sum _{l=1}^d({\mathit{\mu }}_{q,l}^2+{\mathit{\sigma }}_{q,l}^2)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}\sum _{l=1}^d({\mathit{\mu }}_{q,l}^2+{\mathit{\sigma }}_{q,l}^2-log\left({\mathit{\sigma }}_{q,l}^2\right)-1)\hfill \end{array}$$

(A29)

where we have used the fact that for a diagonal covariance matrix ${\Sigma }_q=\mathrm{diag}({\mathit{\sigma }}_{q,1}^2,\dots ,{\mathit{\sigma }}_{q,d}^2)$, the determinant $|{\Sigma }_q|={\prod }_{l=1}^d{\mathit{\sigma }}_{q,l}^2$, and thus $log|{\Sigma }_q|={\sum }_{l=1}^{d}log\left({\mathit{\sigma }}_{q,l}^2\right)$.

Appendix A.6. KL Divergence for the Spike-and-Slab Prior

For the view-specific latent ${\mathbf{Z}}_i^{\left(k\right)}$, the prior is a spike-and-slab: $p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)=(1-{\gamma }_k){\delta }_{\mathbf{0}}\left({\mathbf{Z}}_i^{\left(k\right)}\right)+{\gamma }_k\mathcal{N}({\mathbf{Z}}_i^{\left(k\right)}\mid \mathbf{0},\mathbf{I})$

The variational posterior is a Gaussian: $q\left({\mathbf{Z}}_i^{\left(k\right)}\right)=\mathcal{N}({\mathit{\mu }}_q^{\left(k\right)},\mathrm{diag}\left({\left({\mathit{\sigma }}_q^{\left(k\right)}\right)}^2\right))$

The KL divergence $\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))$ is not standard because the prior is a mixture. We handle this by taking the expectation over the discrete variable ${\gamma }_k$.

$$\begin{array}{cc}\hfill {\mathbb{E}}_{q\left({\gamma }_k\right)}\left[\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))\right]& ={\mathbb{E}}_{q\left({\gamma }_k\right)}\left[\int q\left({\mathbf{Z}}_i^{\left(k\right)}\right)log\left({\displaystyle \frac{q\left({\mathbf{Z}}_i^{\left(k\right)}\right)}{p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)}}\right)d{\mathbf{Z}}_i^{\left(k\right)}\right]\hfill \\ \hfill & =q({\gamma }_k=0)·\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel {\delta }_{\mathbf{0}})+q({\gamma }_k=1)·\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel \mathcal{N}(\mathbf{0},\mathbf{I}))\hfill \end{array}$$

(A30)

The term $\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel {\delta }_{\mathbf{0}})$ is problematic because the delta function is not a density with respect to the Lebesgue measure. In practice, as noted in the main text, when ${\gamma }_k=0$, the prior forces ${\mathbf{Z}}_i^{\left(k\right)}=\mathbf{0}$. Therefore, the model approximates this by only applying the KL penalty against the Gaussian slab when the view is included. This leads to the practical approximation:

$$\begin{array}{c}\hfill {\mathbb{E}}_{q\left({\gamma }_k\right)}\left[\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))\right]\approx {\tau }_k·\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel \mathcal{N}(\mathbf{0},\mathbf{I}))\end{array}$$

(A31)

where ${\tau }_k=q({\gamma }_k=1)$. This approximation is used in Equation (20).

Appendix A.7. Derivation of the Variational Objective for the Spike-and-Slab Prior

The exact KL divergence for the view-specific latent, marginalized over ${\gamma }_k$, is:

$$\begin{array}{cc}\hfill \mathrm{KL}\left(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right)& =\int q\left({\mathbf{Z}}_i^{\left(k\right)}\right)log{\displaystyle \frac{q\left({\mathbf{Z}}_i^{\left(k\right)}\right)}{p\left({\mathbf{Z}}_i^{\left(k\right)}\right)}}d{\mathbf{Z}}_i^{\left(k\right)}\hfill \\ \hfill & =\int q\left({\mathbf{Z}}_i^{\left(k\right)}\right)logq\left({\mathbf{Z}}_i^{\left(k\right)}\right)d{\mathbf{Z}}_i^{\left(k\right)}-\int q\left({\mathbf{Z}}_i^{\left(k\right)}\right)log\left[\int p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k)p\left({\gamma }_k\right)d{\gamma }_k\right]d{\mathbf{Z}}_i^{\left(k\right)}\hfill \\ \hfill & =-\mathcal{H}\left[q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right]-{\mathbb{E}}_{q\left({\mathbf{Z}}_i^{\left(k\right)}\right)}\left[log\left((1-{\mathit{\pi }}_k){\delta }_{\mathbf{0}}\left({\mathbf{Z}}_i^{\left(k\right)}\right)+{\mathit{\pi }}_k\mathcal{N}({\mathbf{Z}}_i^{\left(k\right)}\mid \mathbf{0},\mathbf{I})\right)\right]\hfill \end{array}$$

(A32)

where $\mathcal{H}\left[q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\right]$ is the entropy of the variational distribution.

The intractability arises from the log of a sum inside the expectation. The common variational approach is to introduce an auxiliary variational distribution $r\left({\gamma }_k\right)$ for the discrete variable and derive a lower bound. However, in the context of the overall ELBO, this can become very complex.

The approximation used in the paper, ${\mathbb{E}}_{q\left({\gamma }_k\right)}\left[\mathrm{KL}(q\left({\mathbf{Z}}_i^{\left(k\right)}\right)\parallel p({\mathbf{Z}}_i^{\left(k\right)}\mid {\gamma }_k))\right]$, is itself a variational approximation that treats ${\gamma }_k$ as independent. This is known as the “fully-factorized” or “mean-field” assumption across the discrete and continuous latents. The approximation in Equation (A31) (using ${\tau }_k·\mathrm{KL}(q\parallel \mathcal{N})$) is a further simplification that ignores the KL divergence to the spike (delta function) because it is infinite for any non-degenerate q. This is justified by the model’s structure: when ${\gamma }_k=0$, the survival head $g_{\theta }$ receives $\mathbf{0}$, so the model is incentivized to set $q\left({\mathbf{Z}}_i^{\left(k\right)}\right)$ to ${\delta }_{\mathbf{0}}$ to minimize the reconstruction loss, making the infinite KL penalty moot. The ${\tau }_k$ weighting then acts as a regularizer only when the view is likely to be included.

Appendix A.8. Derivation of the Optimal Variational Distribution for π_k

The appendix mentions a variational distribution $q\left({\mathit{\pi }}_k\right)$ but does not derive its optimal form. For a conjugate model, we can find a closed-form update.

From the joint distribution in Equation (A7), the conditional posterior for ${\mathit{\pi }}_k$ (before variational approximation) is proportional to:

$$\begin{array}{cc}\hfill p({\mathit{\pi }}_k\mid {\gamma }_k)& \propto p({\gamma }_k\mid {\mathit{\pi }}_k)p\left({\mathit{\pi }}_k\right)\hfill \\ \hfill & \propto {\mathit{\pi }}_k^{{\gamma }_k}{(1-{\mathit{\pi }}_k)}^{1-{\gamma }_k}·{\mathit{\pi }}_k^{{\alpha }_0-1}{(1-{\mathit{\pi }}_k)}^{{\beta }_0-1}\hfill \\ \hfill & \propto {\mathit{\pi }}_k^{{\alpha }_0+{\gamma }_k-1}{(1-{\mathit{\pi }}_k)}^{{\beta }_0+1-{\gamma }_k-1}\hfill \end{array}$$

(A33)

This shows that the conditional posterior for ${\mathit{\pi }}_k$ is $\mathrm{Beta}({\alpha }_0+{\gamma }_k,{\beta }_0+1-{\gamma }_k)$.

In variational inference, the optimal $q^{\star }\left({\mathit{\pi }}_k\right)$ that maximizes the ELBO is the one that minimizes the KL divergence to the true conditional posterior. Under the mean-field assumption, the optimal $q^{\star }\left({\mathit{\pi }}_k\right)$ is:

$$\begin{array}{cc}\hfill q^{\star }\left({\mathit{\pi }}_k\right)& \propto exp\left({\mathbb{E}}_{q\left({\gamma }_k\right)}\left[logp({\mathit{\pi }}_k\mid {\gamma }_k)\right]\right)\hfill \\ \hfill & \propto exp\left({\mathbb{E}}_{q\left({\gamma }_k\right)}\left[log\mathrm{Beta}({\mathit{\pi }}_k\mid {\alpha }_0+{\gamma }_k,{\beta }_0+1-{\gamma }_k)\right]\right)\hfill \end{array}$$

(A34)

Since the expectation of ${\gamma }_k$ under q is ${\tau }_k$, we can substitute to find:

$$\begin{array}{cc}\hfill q^{\star }\left({\mathit{\pi }}_k\right)& \propto {\mathit{\pi }}_k^{{\alpha }_0+{\tau }_k-1}{(1-{\mathit{\pi }}_k)}^{{\beta }_0+1-{\tau }_k-1}\hfill \\ \hfill & =\mathrm{Beta}({\mathit{\pi }}_k\mid {\alpha }_0+{\tau }_k,{\beta }_0+1-{\tau }_k)\hfill \end{array}$$

(A35)

Therefore, the optimal variational distribution for ${\mathit{\pi }}_k$ is $q\left({\mathit{\pi }}_k\right)=\mathrm{Beta}({\alpha }_0+{\tau }_k,{\beta }_0+1-{\tau }_k)$. The corresponding KL divergence term $\mathrm{KL}(q\left({\mathit{\pi }}_k\right)\parallel p\left({\mathit{\pi }}_k\right))$ in Equation (A13) can then be computed in closed form using the KL divergence between two Beta distributions.

Appendix A.9. Graph Convolutional Layer: Signal Processing Interpretation

The propagation rule ${\mathbf{H}}^{(l+1)}=\mathrm{ReLU}\left({\widehat{\mathbf{A}}}^{\left(k\right)}{\mathbf{H}}^{\left(l\right)}{\mathbf{W}}^{\left(l\right)}\right)$ can be analyzed from a signal processing perspective on graphs.

The normalized adjacency matrix ${\widehat{\mathbf{A}}}^{\left(k\right)}$ can be decomposed via its eigenvectors: ${\widehat{\mathbf{A}}}^{\left(k\right)}=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^{\top }$, where $\mathbf{\Lambda }$ is a diagonal matrix of eigenvalues and $\mathbf{U}$ is the matrix of eigenvectors. This is the graph Fourier transform.

A graph convolution in the spatial domain (like the GCN layer) is equivalent to a multiplication in the spectral domain. The operation ${\widehat{\mathbf{A}}}^{\left(k\right)}{\mathbf{H}}^{\left(l\right)}$ is a form of low-pass filter. It smooths the feature representation of each node by averaging it with the features of its neighbors. The eigenvalues in $\mathbf{\Lambda }$ are in the range $[-1,1]$, and the largest eigenvalues correspond to the lowest frequencies (most global, smoothest signals on the graph).

The weight matrix ${\mathbf{W}}^{\left(l\right)}$ then acts as a learnable filter that can amplify or attenuate different frequency components. The ReLU non-linearity introduces the capacity to learn more complex, non-linear feature interactions.

This mathematical interpretation justifies why GCNs are effective for biological data: they promote smoothness of the latent representations with respect to the prior biological network, meaning that genes/proteins that are known to interact will have similar representations in the latent space, leading to more biologically plausible and robust models.

Appendix A.10. Hard-Concrete Relaxation for γ_k

The main text mentions using a “hard-concrete relaxation with a straight-through gradient estimator” for ${\gamma }_k$. Here, we provide the mathematical formulation of this trick.

The binary variable ${\gamma }_k\in \{0,1\}$ is replaced during training with a continuous relaxation ${\tilde{\gamma }}_k\in (0,1)$ sampled from a Hard-Concrete distribution. The sampling process is:

$$\begin{array}{cc}\hfill u& \sim \mathrm{Uniform}(0,1)\hfill \\ \hfill s& =\mathit{\sigma }\left({\displaystyle \frac{1}{\beta }}\left(logu-log(1-u)+log\alpha \right)\right)\hfill \\ \hfill \tilde{s}& =s·({\gamma }_{\max }-{\gamma }_{\min })+{\gamma }_{\min }\hfill \\ \hfill {\tilde{\gamma }}_k& =min(1,max(0,\tilde{s}))\hfill \end{array}$$

(A36)

where $\alpha$ is the location parameter (often derived from the logit of ${\tau }_k$), $\beta$ is a temperature parameter, and $[{\gamma }_{\min },{\gamma }_{\max }]$ is a small interval (e.g., $[-0.1,1.1]$) to stretch the output.

During the forward pass, ${\tilde{\gamma }}_k$ is used. For the backward pass (gradient computation), the “straight-through” estimator is used: the gradient is computed as if the operation was an identity function. That is, if $y=\mathrm{hard\_concrete}\left(x\right)$, then $\partial y/\partial x=1$.

This allows gradients to flow through the discrete sampling process, enabling end-to-end training. As training progresses, the temperature $\beta$ is annealed to zero, causing the Hard-Concrete distribution to collapse to a Bernoulli distribution, effectively discretizing ${\tilde{\gamma }}_k$ back to ${\gamma }_k$.

Appendix A.11. Derivation of Gradients for Model Parameters

The appendix derives the ELBO with respect to the variational parameters. It is also crucial to derive the gradients with respect to the model parameters $\psi$ (decoder) and $\theta$ (survival head).

The ELBO is:

$$\begin{array}{c}\hfill \mathcal{L}={\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )]-{\mathbb{E}}_q[logq(\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )]\end{array}$$

(A37)

The second term (the entropy of q) does not depend on the model parameters $\psi ,\theta$. Therefore, the gradient with respect to $\psi$ and $\theta$ comes only from the first term, the expected joint log-likelihood.

$$\begin{array}{c}\hfill {\nabla }_{\psi ,\theta }\mathcal{L}={\nabla }_{\psi ,\theta }{\mathbb{E}}_q[logp(\mathbf{X},\mathbf{T},\mathbf{\Delta },\mathbf{Z},\left\{{\mathbf{Z}}^{\left(k\right)}\right\},\gamma )]\end{array}$$

(A38)

This expectation can be broken down as shown in Equation (A8). The gradients for the reconstruction term ${\mathcal{L}}_{\mathrm{recon}}$ with respect to $\psi$ are derived in subapp:ziggradients. The gradients for the survival term ${\mathcal{L}}_{\mathrm{surv}}$ with respect to $\theta$ are derived in Appendix A.4.

The key point is that these gradients can be estimated using Monte Carlo sampling from the variational distribution q, combined with the reparameterization trick for the Gaussian latents and the straight-through estimator for the discrete latents, making the entire model trainable via standard stochastic gradient descent.

Appendix A.12. Reparameterization Trick for Gaussian Latents

To obtain low-variance, unbiased gradient estimates for the Gaussian latent variables, we use the reparameterization trick.

For a latent variable $\mathbf{z}\sim \mathcal{N}(\mathit{\mu },\Sigma )$, instead of sampling $\mathbf{z}$ directly, we sample an auxiliary noise variable $ϵ\sim \mathcal{N}(\mathbf{0},\mathbf{I})$ and compute $\mathbf{z}=\mathit{\mu }+{\Sigma }^{1/2}ϵ$, where ${\Sigma }^{1/2}$ is the matrix square root. For a diagonal covariance $\Sigma =\mathrm{diag}({\mathit{\sigma }}_1^2,\dots ,{\mathit{\sigma }}_d^2)$, this simplifies to element-wise multiplication: $\mathbf{z}=\mathit{\mu }+\mathit{\sigma }\odot ϵ$, where ⊙ denotes the Hadamard (element-wise) product.

The gradient of an expectation ${\mathbb{E}}_{q\left(\mathbf{z}\right)}\left[f\left(\mathbf{z}\right)\right]$ with respect to $\mathit{\mu }$ and $\mathit{\sigma }$ can then be computed as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial \mathit{\mu }}}{\mathbb{E}}_{q\left(\mathbf{z}\right)}\left[f\left(\mathbf{z}\right)\right]& ={\mathbb{E}}_{ϵ}\left[{\displaystyle \frac{\partial f\left(\mathbf{z}\right)}{\partial \mathbf{z}}}·{\displaystyle \frac{\partial \mathbf{z}}{\partial \mathit{\mu }}}\right]={\mathbb{E}}_{ϵ}\left[{\displaystyle \frac{\partial f\left(\mathbf{z}\right)}{\partial \mathbf{z}}}\right]\hfill \end{array}$$

(A39)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}{\mathbb{E}}_{q\left(\mathbf{z}\right)}\left[f\left(\mathbf{z}\right)\right]& ={\mathbb{E}}_{ϵ}\left[{\displaystyle \frac{\partial f\left(\mathbf{z}\right)}{\partial \mathbf{z}}}·{\displaystyle \frac{\partial \mathbf{z}}{\partial \mathit{\sigma }}}\right]={\mathbb{E}}_{ϵ}\left[{\displaystyle \frac{\partial f\left(\mathbf{z}\right)}{\partial \mathbf{z}}}·ϵ\right]\hfill \end{array}$$

(A40)

This allows gradients to flow through the sampling process, enabling efficient optimization via stochastic gradient descent.