Year: 2021

Abstract:
The study of complex systems requires the integration of multiple heterogeneous and high-dimensional data types (e.g. multi-omics). However, previous generative approaches for multi-modal inputs suffer from two shortcomings. First, they are not stochastic processes, leading to poor uncertainty estimations over their predictions. This is mostly due to the computationally intensive nature of traditional stochastic processes, such as Gaussian Processes (GPs), that makes their applicability limited in multi-modal learning frameworks. Second, they are not able to effectively approximate the joint posterior distribution of multi-modal data types with various missing patterns. More precisely, their model assumptions result in miscalibrated precisions and/or computational cost of sub-sampling procedure. In this paper, we propose a class of stochastic processes that learns a graph of dependencies between samples across multi-modal data types through adopting priors over the relational structure of the given data modalities. The dependency graph in our method, multi-modal Relational Neural Process (mRNP), not only posits distributions over the functions and naturally enables rapid adaptation to new observations by its predictive distribution, but also makes mRNP scalable to large datasets through mini-batch optimization. We also introduce mixture-of-graphs (MoG) in our model construction and show that it can address the aforementioned limitations in joint posterior approximation. Experiments on both toy regression and classification tasks using real-world datasets demonstrate the potential of mRNP for offering higher prediction accuracies as well as more robust uncertainty estimates compared to existing baselines and state-of-the-art methods.

heading: 1 INTRODUCTION
text: Many prominent methods for multi-modal learning (Argelaguet et al., 2018; Andrew et al., 2013; Klami et al., 2013; Zhao et al., 2016) are based on the canonical correlation analysis (CCA), which extracts shared components across multiple views. The main idea is that, given two vectors of random variables, the method finds the linear projections in a shared latent space, in which the projected vectors are maximally correlated (Thompson, 1984). This can help to understand the overall dependency structure between these two vectors. However, the classical CCA cannot handle non-linearity (Andrew et al., 2013) and suffers from a lack of probabilistic interpretation when applied to high dimensional data (Klami et al., 2013). To address these issues, probabilistic CCA (PCCA) has been proposed and extended to non-linear settings using kernel methods and neural networks (Bach & Jordan, 2005).
More recently, Wu & Goodman (2018) introduced the multimodal variational autoencoder (mVAE) that models the joint posterior as a product-of-experts (PoE) over the marginal posteriors, enabling cross-modal generation at test time without requiring additional inference networks and multi-stage training regimes. However, training a PoE is difficult and the following techniques are needed in order to ensure that the individual views are learnt faithfully; i) using artificial sub-sampling of the observed views (Wu & Goodman, 2018), ii): applying variants of contrastive divergence (Hinton, 2002), and iii): utilizing information bottleneck on the marginal representations of each view (Lee & Schaar, 2021). Furthermore, each expert holds the power of veto and low density of only one marginal posteriors among a given set of observations leads to the low density of joint distribution (Shi et al., 2019). In the case of Gaussian experts, different levels of complexity of modalities or
sub-optimal initialisation might result in the miscalibrated precisions that potentially lead to a biased overall mean prediction (Shi et al., 2019).
Despite the success of the existing multi-modal generative methods, they suffer from inducing stochasticity in their inferred functions. This leads to; i) poor adaptation to new observations1, ii) inefficient generalizability in few-shot settings and, iii) inaccurate uncertainty estimation, all of which are important in many safety critical applications such as healthcare. On the other hand, stochastic processes, such as Gaussian processes (GPs), are exchangeable models that posit distributions over possible functions, and are updated in light of data through the probabilistic inference. Despite these advantages, GPs are computationally expensive and their underlying model is not flexible for high-dimensional inputs, making them infeasible for multi-modal learning settings with multiple heterogeneous and high-dimensional data types.
Apart from the aforementioned issues, existing generative models for multi-modal learning focus on latent representation, but do not fully incorporate the label information. In order to facilitate predictions on the target labels, they need to apply two-step procedures, in which the shared latent spaces are used for downstream classification. While information relevant to the reconstruction of views will be well-captured in the shared representations, discriminative information relevant to the target task may be discarded, leading to poor prediction performance of such models (Lee & Schaar, 2021). Furthermore, the shared representations of these methods in high-dimensional space are often difficult to interpret and need an extra downstream analyses such as PCA or t-SNE (Van Der Maaten, 2014) that may cause biased results (Lötsch & Ultsch, 2020).
In this paper, instead of following the existing multi-modal generative methods that rely on deterministic predictive functions, we propose the first multi-modal stochastic processes family that learns distributions over functions for inputs with any missing patterns. Like GPs, it naturally provides robust uncertainty estimates and can encode inductive biases. Unlike GPs, though, mRNP is computationally efficient during training and evaluation since it can take advantage of mini-batch optimization, and learns to adapt their priors to data. The inferred dependency graph in our model construction can be used to visualise the high dimensional latent representation, obviating the need for ad-hoc post-processing steps as required in most of the multi-modal learning methods.
The presented work makes four major contributions: 1) We develop a novel multi-modal Relational Neural Process, mRNP, that defines a distribution over functions of multiple data types by employing local latent variables, and learns a dependency structure among the samples of the given modalities. 2) We theoretically prove exchangeability and consistency of mRNP, two necessary conditions that have to be satisfied during the construction of such a model, showing that mRNP is a valid stochastic process. 3) We further show that the local latent variable structure in mRNP is able to encode inductive biases and demonstrate this by designing an mRNP model that behaves similarly to a GP with an RBF kernel (an ablation study). 4) We introduce mixture-of-graphs (MoG) in our model construction that can address the issues such as computational complexity and miscalibrated precisions observed in the previous approaches in multi-modal learning.

heading: 2 METHOD
text: We propose a new graph-structured multi-modal learning method, referred as multi-modal Relational Neural Process (mRNP), that combines the benefits of neural networks with that of stochastic processes. In the supervised multi-modal setting, our dataset contains sets of feature-label pairs, {x̄, y} from |V| different modalities with various missing patterns, where |V| is number of views, x̄ = {xv|vth modality present}v∈V ∈ X̄ is the input modalities and y ∈ Y is the given label. One of the key motivation of mRNP is that it gives us the ability to define a distribution over functions rather than learning a single static function when fitting the data. In this framework, inspired by the idea of inducing points in sparse Gaussian processes (Titsias, 2009; Damianou & Lawrence, 2013), we first select a reference set of samples that consists of functions f : x̄→ y that are sampled from some underlying distributionH. We then establish a probability distribution over fh around those samples. We define the reference set, XR = {XR,1, . . . ,XR,V } ⊂ X̄ , where XR,v ∈ RNv×Dv , and
1Rapid adaption to new observations is what distinguishes a model that can learn a whole distribution over functions by seeing different functions in each batch from the one that learns a single underlying function by observing several batches during the training. At the test time, the former will narrow down what the current function is by seeing a few context observations.
Nv and Dv are the number of samples and dimension of features in the domain v, respectively. We also define O = X̄ \XR, which is the set of all possible samples that are not in XR. Please note that two sets XR and O includes both labeled and unlabeled samples. We also denote our labeled samples as Dx = {x̄1, . . . , x̄N}, a set of any finite random set from X̄ , as well as defining sets XM = Dx\XR and XB = XR ∪ XM . The Venn diagram of the sets are provided in Figure 1. We infer the parameters of our probabilistic model using variational inference. In the following, we first introduce each of the latent variables in our model as well as their corresponding prior and posterior distributions, and then prove that it corresponds to an infinitely exchangeable stochastic process. The graphical model of mRNP is illustrated in Figure 1.
Modality-specific latent spaces. The first step is to construct modality-specific latent representations, denoted by U = {Uv}v∈V , a set of Nv ×Du matrices, which independently embed samples of each modality (i.e. Nv) into a Du dimensional latent space as follows:∫ pθ(UB , XB) dX =
∏ v∈V ∫ pθ(UB,v , XB,v) dXB,v = ∏ v∈V ∫ pθ(UB,v |XB,v) p(XB,v) dXB,v.
We further parameterize the distribution over the samples ui,v independently: pθ(UB,v |XB,v) = ∏ i∈Bv pθ(ui,v |xi,v), (1) where pθ(ui,v |xi,v) can be any distribution that is derived from the observed input data. For simplicity, we use a diagonal Gaussian, the parameters of which are a function of the observed input. More specifically, we use two functions denoted by ϕemb,µv (XB,v) and ϕ emb,σ v (XB,v) to infer the mean and variance of the distribution for each modality. Depending on the nature of the modality, these functions can be implemented using any highly expressive functions such as many variants of neural networks (NNs) and graph NNs.
A directed graph among the reference samples. The next step is to construct a directed graph across the samples inXR using different modalities. Given the latent embedding UR, we first construct a set of directed acyclic graphs (DAGs), G = {Gv}v∈V , where Gv is a random binary adjacency matrix between reference samples present in the vth modality. Then we combine the graphs Gv from multiple views to construct a common graph G. Inspired by the concept of stochastic orderings (Shaked & Shanthikumar, 2007), we impose a topological ordering over the vectors in UR to avoid cycles in each individual Gv . The distribution of the adjacency matrices are defined as follows:
p(Gv |UR,v) = ∏ i∈Rv ∏ j∈Rv,j 6=i Bernoulli ( Gij,v | I[t(ui,v) > t(uj,v)]ϕsim(ui,v , uj,v) ) . (2)
Similar to Louizos et al. (2019), we employ a parameter free scalar projection t(ui,v) = ∑ k tk(uik,v), where tk(·) is a monotonic function. In practice, we use the log cumulative distribution function of a standard normal distribution. Given the graphs Gv, we construct a directed graph G, in which the weight of each edge is average of its weights in different modalities, i.e. G = ∑ v∈V βvGv/ ∑ v∈Vmvm T v , where mv is a binary vector and mi,v = 1 when sample i is observed in the modality v, and βv = 1/|V |, assuming that the all modalities are equally important. The proposed formulation, referred as mixture-of-graphs (MoG), is used to learn the shared structured representation of multi-modal inputs by using the graphs of individual modalities.
A relational graph of dependencies from reference samples to XM . In addition to the directed graph G used for the reference samples, we apply MoG to construct a bipartite graph
A = ∑ v∈V αvAv that is a bi-adjacency matrix from XR to XM . Given the latent embedding UB as described previously, we model the elements of each individual bi-adjacency matrix by Bernoulli distribution as follows:
p(Av |UM,v , UR,v) = ∏ i∈Mv ∏ j∈Rv Bernoulli ( Aij |ϕsim(ui,v , uj,v) ) , (3)
where ϕsim(·, ·) is a score function measuring the similarity between the latent representations of the input samples in each modality v, and αv = 1/|V |, assuming the all modalities equally important. Depending on the desired relational inductive biases, we can appropriately define the function ϕsim(·, ·). The Bernoulli-Poisson link ϕsim(ui,v,uj,v) = 1− exp(− ∑Du k=1 τk uik,v ujk,v) and ϕsim(ui,v,uj,v) = exp(− τ2 ||ui,v − uj,v|| 2) are two examples of such potential score functions. Please note that we use the latter one throughout this paper.
For both A = {Av}v∈V and G, we use the concrete relaxation (Maddison et al., 2016; Gal et al., 2017) during training while we sample from Bernoulli distributions in the testing phase.
A shared latent variable. Having obtained the sample representations U and the dependency graphs A and G, we construct a shared latent variable, denoted by N ×Dz matrix Z, which can be used to predict target variable distributions yi. We parameterize the distributions of the predictive targets as follows:∫ pθ(yB ,ZB |A,G,XR)dZB = ∫ ∫ pθ(yR,ZR |G,XR) pθ(yM ,ZM |A,XR,yR) dZR dZM
= ∏ i∈R ∫ pθ(zi |parGi(XR,yR)) pθ(yi | zi) dzi ∏ j∈M ∫ pθ(zj |parAj (XR,yR)) pθ(yj | zj) dzj ,
where parGi(·), parAj (·) are functions that return the parents of the points i and j according to G and A, respectively. We summarize the information from the parent samples and their targets in XR through the local latent variable zi, allowing the distribution yi to be explicitly dependent on the available graph dependency structures A and G. We set the distribution over zi as an independent Gaussian distribution whose parameters are a function of either A or G although any distribution with a permutation invariant probability density with respect to the parents can be used. More formally, the following distribution is adopted:
pθ(zi |parAi(XR,yR)) = Dz∏ k=1 pθ(zik |parAi(XR,yR)) = Dz∏ k=1 N (µpriorik , σ prior ik ), (4)
µpriorik = ci ∑ j∈XR Aijϕ prior,µ k (xj,R , yj,R), σ prior ik = exp(ci ∑ j∈XR Aijϕ prior,σ k (xj,R , yj,R)),
where xj,R = {xj,R,v}v∈V , ϕprior,µ and ϕprior,σ are transform functions with a co-domain R|z|, and ci = ( ∑ j Aij + )
−1 is a normalization constant. Given input {XB,v}v∈V , we factorize the variational posterior distribution as
qφ(Z | XB) = ∏ i∈B qφ(zi |xi,B), (5)
In practice, we also parameterize the priors over the latent zi in terms of the posterior distribution of the reference samples, i.e. qφ(zi |xi,R). More precisely, we define ϕprior,µ and ϕprior,σ in equation (4) as ϕpost,µ(xi,R) +ϕ emb,µ label (yi,R) and ϕ post,σ(xi,R) +ϕ emb,σ label (yi,R), respectively; where ϕpost and ϕemblabel provide the means and variances of qφ(zi |xi,R) and the linear embedding of the labels of reference samples, respectively.
We should point out that in this formulation, the prediction of yi indirectly depends on the input samples xi through the graphs G and A, which are functions of ui. This leads to an uninformative standard normal prior over zi for the samples with very small probability of being connected to the reference set through A, and hence the prediction of yi will be constant. This can be seen as encoding inductive bias similar to a GP with an RBF kernel.
MoG vs PoE. PoE is a prime choice to learn the variational joint posteriors; it can learn under any combination of missing modalities (Shi et al., 2019). Despite this advantage, it also has two main limitations: 1) the underlying model suffers from the overconfident experts since experts with
greater precision will have more influence over the combined prediction than experts with lower precision; 2) training and inference can be costly due to its artificial sub-sampling procedure. Given the aforementioned limitations of PoE, we propose to build a graph of dependencies among local latent variables as a mixture of modality-specific graphs, i.e. mixture-of-graph (MoG), and then parameterize the joint latent space Z conditioned on those graph. Unlike PoE, MoG effectively spreads its density over all the individual modalities by imposing a weight to each individual graphs. When prior knowledge is available, terms αv and βv in the MoG construction can be used to encourage certain modalities, which is the case in multi-omics data integration where individual modalities often contain complementary information of target task. A more detailed discussion can be found in Section D of the supplement.
We also should point out that mRNP is not a graph learning model. Rather, we propose a novel multi-view NPs that is able to solve a pitfall in utilizing stochastic processes for multi-view setting by learning a graph of dependencies. We can also consider graph learning in mRNP as a kind of cross-attention in the form of a dependency graph among local latent variables. From this viewpoint, mRNP is learning attention mechanism to combine modalities. The current framework in which mRNP learns the individual graphs might not be computationally efficient, especially if training set is large. We may improve the scalability issue of graph construction part of our mRNP by e.g. considering hierarchical structure similar to hFGW (Xu et al., 2020). We leave this for a future study.
Overall likelihood and learning. Putting everything together, the marginal likelihood is pθ(yB | XB) = ∑ G,A ∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR) dUB dZB . (6) We should point out that, for the sake of simplicity, we choose the reference set to be part of the training set Dx throughout this paper. More specifically, we assume XB = Dx. In case of existing reference samples that are not part of the training set, we need to marginalize out over them in equation (6) in addition to the marginalizations over the latent variables and structured graphs. We can also use the marginalization technique to incorporate unlabelled data in order to learn a better representation and/or to impute the missing labels. We leave this for a future study.
We deploy variational inference to optimize the model parameters θ and variational parameters φ by minimizing the following derived Evidence Lower Bound (ELBO): L = Eqφ(ZR | XR) p(G | UR) ∏v∈V pθ(UR,v |XR,v) [log pθ(yR,ZR |G,XR)− log qφ(ZR | XR)] +
Eqφ(ZM | XM ) p(A | UB) ∏ v∈V pθ(UB,v |XB,v)[log pθ(yM |ZM )+
log pθ(ZM |parA(XR,yR))− log qφ(ZM | XM )].
Please note that we can use mini-batches for the second expectation term where the size of the batches scale according to the size of the reference set XR. Due to the DAG structure in the first expectation term, we cannot decompose it to independent sums.
Predictive distribution. After optimizing the model parameters, we derive the predictive distribution for unseen samples {x∗v}v∈V with missing modalities based on Bayes theorem as follows:
pθ(y ∗ | {x∗v}v∈V ,XB ,yB) = ∑ a∗ ∫ ∏ v∈V pθ(u ∗ v,UR,v |x∗v,XR,v) p(a∗ | {u∗v}v∈V ,UR) pθ(z∗ |para∗(XR,yR)) pθ(y∗ | z∗) dUR du∗v∈V dz∗, where u∗v is the embedding representation of the observed modality x ∗ v through the neural network ϕembv and a ∗ is a binary vector similar to a row of A that denotes reference parents of the new samples. This is similar to the predictive distribution in few-shot learning (Sung et al., 2018).
Proposition I. mRNP corresponds to Bayesian models as the distributions defined in Equation (6) are valid, permutation invariant stochastic processes.
proof. The proof is deferred to the supplement.

heading: 3 RELATED WORKS
text: Multi-modal VAE. There has been a long line of research in using variants of VAEs to explicitly model the joint distribution over latents and data (Suzuki et al., 2016; Vedantam et al., 2017; Wu
& Goodman, 2018; Shi et al., 2019). Most of these works are not able to efficiently explore more than two modalities and/or suffer from multi-stage training. Using PoE, mVAE (Wu & Goodman, 2018) addresses the aforementioned problems and learns a joint embedding under any combination of missing modalities. However, 1) it suffers from the miscalibration of the precision of the experts due to differences in complexity of input modalities, and 2) given a complete dataset with no missing modalities, it requires artificial sub-sampling in order to faithfully learn the individual modalities.
Closely related to the mVAE, DeepIMV (Lee & Schaar, 2021) applies PoE in a supervised framework for multi-omics data integration with missing views. Our mRNP is similar to this method in terms of using supervised framework. But, unlike mRNP, DeepIMV predictive function is deterministic and therefore is not suitable for out-of-distribution prediction and/or few-shot settings. It also includes a set of |V| view-specific predictors in its construction and apply information bottleneck (IB) principle on the marginal representations of individual views in order to allow each encoder to build view-specific expertise for predicting the target. Depending on the number of input modalities, computational complexity might be an issue.
Neural process (NP). NPs (Garnelo et al., 2018b;a) have been recently proposed to learn an approximation of a stochastic process by using a neural network-based formulation. Traditional stochastic processes such as GPs are usually computationally expensive and the available kernels are usually restricted in their functional form and are not very flexible for high dimensional problems (Garnelo et al., 2018b). On the other hand, NPs are data-driven models that impose stochasticity in the function realizations and are able to define a flexible class of stochastic processes well suited for highly non-trivial functions(Lee et al., 2020). The main idea behind NPs is to define explicit global latent variables in their construction which can capture functional uncertainty(Garnelo et al., 2018b). Closely related to our mRNP, Louizos et al. (2019) discards the idea of the global latent variables and instead builds a graph of dependencies among local latent variables, making it more suitable for modeling high-dimensional data. There are other NP variants (Kim et al., 2019; Lee et al., 2020) to improve the performance of NPs. However, none of these works addresses multi-modal learning problem. To the best of our knowledge, mRNP is the first NP-based model to explore multi-modal in a computationally efficient way.

heading: 4 EXPERIMENTS
text: We evaluate the effectiveness of mRNP based on two different multi-modal experiments: 1) On two real-world datasets, we measure the performance of mRNP by using its prediction scores and the quality of uncertainty quantification in multi-modal classification tasks under any combination of missing modalities; 2) We investigate the inductive biases that we can encode in mRNP in an ablation study by visualizing the predictive distributions in two toy examples of 1-dimensional multi-modal regression tasks. We compare its performance with those of three state-of-the-art methods; MOFA (Argelaguet et al., 2018), mVAE (Wu & Goodman, 2018), and DeepIMV (Lee & Schaar, 2021). We should emphasize that MOFA and mVAE models are unsupervised embedding models and suffer from task specific training. Specifically, these models derive low-dimensional (non-)linear embedding of the input samples, and hence a classification model has to be trained for a downstream analysis task. We should also point out that only autoencoder-based methods are able to be trained using the observed samples, for which the optimized parameters can be used for unseen samples. So, for MOFA, we need to combine training and test data in order to learn the low-dimensional embedding.
We implemented our model in PyTorch (Paszke et al., 2019). Throughout the classification experiments, we use two fully connected layers with ReLu as the activation function at each hidden layer to obtain an intermediate hidden representation hv of the inputs xv, and then parameterize a linear output layer to learn the parameters of pθ(uv |xv), which we consider as Gaussian. We also factorize the posterior qφ(Z | XB) = ∏ i∈B ( p(zi) ∏ v∈V q (v) φ (zi |xi,B,v) ) , where p(zi) = N (µ0,Σ0), and
q (v) φ (zi |xi,B,v) = N (µv,Σv) is the underlying inference network for the modality v. We use a linear layer as an encoder for each modality. For fair comparison, the results of the competing methods are obtained based on their original implementations with the same architectures as mRNP. For the regression tasks, we used one fully connected layer with ReLu as the activation function to transfer xv to hv, followed by two separate linear layers to learn the parameters of pθ(uv |xv), and qφ(zv | X ), respectively. We train mRNP for a maximum of 50 (5000) epochs using Adam (Kingma
& Ba, 2014) with a learning rate of 0.001 for classification (regression). More implementation details are included in the supplement.

heading: 4.1 PREDICTION PERFORMANCE AND UNCERTAINTY QUALITY
text: For the classification task, we consider two publicly available real-world, multi-modal datasets; the HandWriting (HW) (Dua & Graff, 2017) and the Cancer Cell Line Encyclopedia (CCLE) (Barretina et al., 2012). Throughout the experiments, the dimensions {Du, Dz, Dh} are {32, 64, 200} and {32, 32, 32} for HW and CCLE, respectively. We use 300 (150) random samples from the training set as XR for HW (CCLE). We train mRNP for 50 (10) epochs for HW (CCLE) and use the validation set for the early stopping. All of our results are averaged over 100 runs, ten randomly drawn dataset splits and ten runs with different random seeds on each split.
In order to capture uncertainty, we choose the predictive entropy that combines both epistemic and aleatoric uncertainties (Mukhoti & Gal, 2018). In particular, we assume that epistemic uncertainty of our Bayesian model will increase in areas where we have no data, i.e. out of distribution (o.o.d.) datasets, while the predictive entropy of in-distribution can be used to compare different models in terms of combination of both uncertainties. In order to determine whether a sample is in or out of distribution based on its predictive entropy, we also report the area under both the receiver operating characteristic (ROC) and precision-recall (PR) curves. We should point out that, by using ROC and PR, we can make sure that the improvement in predictive entropy is not due to a trivial model/learning since the model must have low entropy on the in-distribution test set but high entropy on the o.o.d. datasets in order to perform well in both ROC and PR. For the HW dataset, we consider Gaussian N (0, 1), uniform U [0, 1], and Poisson(λ = max(x)) noises as o.o.d. datasets. For CCLE, we use three real-word cancer datasets COAD, KIRC, and SKCM from The Cancer Genome Atlas (TCGA) (Tomczak et al., 2015) as well as Poisson noise as o.o.d datasets.
Multi-view handwriting data with missing views. First, we consider the classification of multiview handwriting data. This dataset consists of features of handwritten numerals, i.e. {0, . . . , 9}, extracted from a collection of Dutch utility maps (van Breukelen et al., 1998). The dataset includes 76 Fourier coefficients of the character shapes, 47 Zernike moments and 6 morphological features for a total of 2,000 samples (200 samples per class). We follow the same data pre-processing as in Lee & Schaar (2021), and choose 20% and 16% of samples as test and validation sets.
The summary of the results are reported in Table 1. The proposed method outperforms the state-of-the-art (SOTA) methods by a significant margin in terms of both classification and o.o.d prediction accuracies. Our Bayesian mRNP model performs significantly better than DeepIMV, which is the only supervised baseline, in terms of uncertainty quantification. This might be
due to the fact that the softmax entropy in DeepIMV is only able to capture aleatoric uncertainty while features with high epistemic uncertainty are ignored, resulting in a relatively high number of inaccurate but certain samples. mRNP, on the other hand, captures both aleatoric and epistemic uncertainties.
We also compare the shared (structured) representations of different modalities from mRNP with the joint latent representations learned by DeepIMV and mVAE. Figure 2 shows an example of the graph G that mRNP learns, as well as PCA projections of latent representations from DeepIMV and mVAE. First, it shows that mRNP learns a meaningful G by connecting samples with same class to each other. Second, a few samples that are not in the same community as their labeled class, have
lower degrees compared to the other nodes in the same community. This indicates that the model is less confident about the class of such a sample. The proposed mRNP provides more interpretable shared representation than traditional non-structured representations while obviating the need for post-processing step to visualize a high-dimensional embedding space in a 2- or 3-dimensional one. Furthermore, we evaluate the joint representations of these models on the HW dataset through two downstream tasks: clustering and classification. For the first task, we apply K-means clustering and for the second one we train a linear SVM. The normalized mutual information (NMI) and classification accuracies are reported in Table 2.
We further perform an ablation study on HW dataset and provide some additional insights regarding the sensitivity of mRNP to the number of reference samples. Table 3 shows the test accuracy, entropy, o.o.d entropy and AUC for o.o.d prediction. Note that the mean and standard error across all of the o.o.d. datasets are reported for both o.o.d. entropy and AUC. We observe that the choice of samples as a reference set is important for mRNP
to fit the data well and that the performance does not always improve with more samples. While the accuracy of mRNP reduced less than 1% in very high sample sizes, its performance degraded more than 30% with small sample size of reference set, e.g. 50 samples. Moreover, the model is less robust and shows higher in-distribution entropy in this case. We should also point out that we only explore the effect of reference size, not the quality of the samples. Investigating the automatic selection of reference samples may improve scalability and alleviate the dependence of mRNP on acquiring a reasonable XR. We leave this for a future study.
Multi-view CCLE dataset. We further apply mRNP to explore sensitivity of 472 heterogeneous cell lines to the drug Panobinostat. We study the integration of three omics, microRNA sequence (mRNA-seq), reverse phase protein array (RPPA), and metabolites. Following Lee & Schaar (2021), we use ActArea as an indicator of drug sensitivity and then label a response as sensitive when this value is higher than median of all responses, and otherwise assign it to the nonsensitive group. Note that for TCGA datasets, which we used as o.o.d samples, only two of these modalities are available, i.e. mRNA-seq and RPPA, so we consider metabolite as a missing modality. The summary of the results are reported in Table 4. mRNP outperforms SOTAs in terms of F1 score while having higher average AUC and AP on the o.o.d. datasets. Lower in-distribution predictive entropy of mRNP compared to DeepIMV and mVAE indicates that mRNP is a better
choice to deal with high-dimensional data of small sample size, which is often the case in multi-omics data integration when studying complex disease. We further compare the models, based on their robustness to the missing views in the test set. In this ablation study, we train mRNP and two other baselines with training set of CCLE, in which only 3% of the views was missing in the main dataset. At the test time, we artificially introduce missing views by defining a probability of not observing an individual view for each sample. Please note that each view of the samples are treated independently. Then, we increase the probability of views missing from 10% to 50% and measure how different models perform in terms of prediction accuracy. The classification errors of the baselines increase monotonically as the percentage of missing views is increased. mVAE has the worst overall per-
formance, followed by DeepIMV. Using 70% of all the views, the average F1 score of mRNP reduces less than 0.5%, while DeepIMV and mVAE performances degraded around 3%. In the worst case scenario (missing 50% of the views), the mRNP outperforms the DeepIMV and mVAE by 6% and 7%, respectively. This clearly shows the robustness of mRNP to the missing views at test time, and proves that it can faithfully learn individual modalities without any sub-sampling procedure.

heading: 4.2 ABLATION STUDY TO EXPLORE THE INDUCTIVE BIASES IN MULTI-MODAL REGRESSION
text: To validate whether mRNP can encode the inductive biases, we consider two toy multi-modal 1-d regression tasks. Following the single-view regression tasks in Louizos et al. (2019) and Osband et al. (2016), we first draw 12 samples from U [0, 0.6] and 8 samples from U [0.8, 1] as the common space between two modalities, and then parameterize the target by yi = si+ +sin(4(si+ ))+sin(13(si+ )), where ∼ N (0, 0.032). Next, we generate observations for the first and second modality as xv1 = s and xv2 = 6s + , respectively. Please note that we impose 10% chance that the views are missing. Figure 4 shows the results of both modalities, in which mRNP has a very similar behaviour to the GP. Similar to GP, it reports high uncertainty in the areas where there are no observed samples. For the second regression study, we draw 20 samples fromU [−4, 4] for the first modality and simulate the target labels as yi = x3i,v1 + , where ∼ N (0, 92). Then, we generate samples for
the second modality as a nonlinear function xi,v2 = 2xi,v1 + η if xi,v1 ≥ 0 else 0.7xi,v1 , where η ∼ N (0, 0.12). Similar to GP, mRNP shows the tendency to quickly move towards a flat prediction outside the areas with observed samples (section D in the supplement).

heading: 5 CONCLUSION
text: We have proposed a novel Bayesian relational generative method, tailored for multi-modal neural processes with any missing patterns in their input modalities — mRNP. The mRNP learns the dependency graph among the samples of the given modalities and takes advantage of it to define distributions over the functions of multi-modal data types. Unlike traditional stochastic processes, mRNP can be optimized through mini-batch. By introducing the mixture-of-graphs (MoG), we also could address the issues such as e.g. precision miscalibration of experts and computational complexity observed in the available multi-modal generative models. We have evaluated mRNP on two different real-world multi-modal classification tasks, which demonstrate that not only mRNP substantially outperforms SOTA methods in terms of both prediction accuracy and uncertainty quantification, but also it captures a meaningful joint (structured) representation across views. Furthermore, our visualization results on two toy multi-modal regression tasks show that the predictive distributions of mRNP is similar to a GP with an RBF kernel.

heading: A DETAILS ON THE EXPERIMENTAL SETUPS, RESOURCES, AND RUNTIME
text: In both classification tasks, we used early stopping, based on the accuracy on the validation set and no other regularization was employed. In order to avoid sub-optimal initialization in CCLE dataset, we also used the validation set for initialization tuning of the parameters.We also employed a soft-free bits (Chen et al., 2016) modification of the bound to help with the optimization. More precisely, we allowed one free bit on average across all dimensions and batch elements throughout training of mRNP.
During the regression experiments, we randomly select 10 samples from training set as XR and consider the dimensions {Du, Dz} as {3, 80} and {3, 10} for the first and second tasks, respectively. In both cases we use one fully connected layer with ReLu as the activation function to transfer xv to hv , following with a linear layer to learn the parameters of pθ(uv |xv), and a linear layer to learn the parameters of qφ(zv | X ). The dimension of the hidden layer is set to be 100. All of the experiments are run on a single GPU node Tesla K80. Each training epoch of mRNP for HW, CCLE, and toy regression took 1.08, 0.53, and 0.02 seconds on average, respectively.

heading: B OTHER RELATED WORKS
text: Multi-modal VAE (continued). Shi et al. (2019) introduces mixture-of-experts multimodal VAE (mmVAE) that does not suffer from potentially overconfident expert. However, it is not able to handle missing modalities during the training, and thus limited in its applicability, especially for biomedical data of small sample sizes with missing views. Our mRNP combines benefits of both mVAE and mmVAE, which is especially critical to deal with high-dimensional data integration. Compared to these approaches, mRNP can in general be more scalable since it does not require sub-sampling as mVAE does. Instead, it effectively takes a vote amongst the graphs and spreads its density over all the individual graph through αv and βv. While the aforementioned models attempt to learn low-dimensional latent variables for multiple views, the focus of mRNP is to introduce distribution over the functions in a supervised setting.
Conditional VAE (CVAE). CVAE (Sohn et al., 2015) are conceptually similar to NPs and has been used to train multi-modal generative models. Unlike the global latent variables that are present in the NPs, the latent variables of CVAEs are local and their decoder functions are deterministic. Hence, it needs to be separately sampled for each of the target predictions y∗i . This is in contrast to the latent variable of an NP that is only sampled once and used to predict multiple values of y∗. More precisely, conditioning on the labels of training samples is done by adding the dependency both in the prior p(z∗ | yc) and decoder p(y∗ | z∗, yc) leading to a deterministic function of the training samples.
CCA based models. There is a rich literature on the extension of CCA to generative modelling (Bach & Jordan, 2005; Virtanen et al., 2011; Klami et al., 2013; Wang et al., 2015; Argelaguet et al., 2018). These extensions can deal with high dimensional data of small sample size (Argelaguet et al., 2018) and outperform classical CCA in terms of interpretability and expressive power. Some of them, such as (Bach & Jordan, 2005; Klami et al., 2013), are generic factor analysis models that decompose the data into shared and view-specific components, and impose an additional constraint to extract the statistical dependencies between the views. Most of the generative methods such as MOFA (Argelaguet et al., 2018) retain the linear nature of CCA while providing more robust
inference compared to the classical solutions. There are also a number of recent models based on VAE that incorporate non-linearity and keep the probabilistic interpretability of CCA (Virtanen et al., 2011; Gundersen et al., 2019). However, they require additional computation to handle missing modalities(Wu & Goodman, 2018).

heading: C MRNP AS STOCHASTIC PROCESS
text: Proposition I. mRNP corresponds to Bayesian models as the distributions defined in Equation (7) are valid, permutation invariant stochastic processes.
proof. Similar to Garnelo et al. (2018b) and Louizos et al. (2019), we rely on de Finetti’s and Kolmogorov Extension Theorems Klenke (2013) that states exchangeability and consistency as two sufficient conditions to define a stochastic process.
Exchangeability. This condition requires p(yB | XB) to be permutation invariant. More precisely, if each probability density of the model be permutation invariant, then the overall probability will be permutation equivariant. Based on equation (1) in the main paper, we can consider ūi = fu(x̄i), where ūi = {ui,v|vth modality present}v∈V and x̄i = {xi,v|vth modality present}v∈V . Therefore, we can show that the latent variables are permutation equivariant with respect to XB as
fu(σ(XB)) = σ(fu(XB)),
where σ(·) is a permutation function. Then, mRNP defines each elements of the shared graphs A and G, denoted as ai,j and gi,j , as a function of (ūi, ūj). Therefore, we also have permutation equivariance for the rows and columns of these two graphs. Based on equation (5) in the main paper, mRNP defines the distribution over zi as an independent distribution whose parameters are a function of either A or G, therefore, the distribution of zi is invariant to its parents’ permutations. Therefore, the permutation of XB will result in the same re-ordering of the ZB as follows
σ(ZB) = fz(σ(XB)), (7)
where fz : XB −→ ZB . Let’s also define fy : zi −→ yi. As a result, we can also show σ(yB) = fy(σ(XB)). The product, integral, and summation operators are all permutation invariant. And since we show all random variables of our mRNP model are permutation equivariant to XB , mRNP model is permutation invariant.
Consistency. This condition requires to show that if we marginalise out a part of p(yB | XB), the resulting marginal distribution is the same as that defined on the original space. To this end, we follow Louizos et al. (2019) and define X̃B = XB ∪ {x̄0, y0}. Then, we need to show p(yB | XB) =∫ p(yB̃ | X̃B) dy0 in two cases; where {x̄0, y0} belongs to either the training set Dx or the reference set that are not part of the training set, i.e. XR\Dx. In the first case, where {x̄0, y0} ∈ Dx, the new sample will be added to a leaf in the dependency graph. Therefore, it will not affect any of the samples in XR, hence we can marginalize it out as follow:∫
p(yB̃ | X̃B) dy0 = ∑
G,A,a0
∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR)
pθ(u0,v |x0,v)p(a0 | ū0,UR) pθ(z0 |para0(XR,yR)) (∫ pθ(y0 | z0)dy0 ) dUB dZB dū0 dz0. As ∫ pθ(y0 | z0)dy0 = 1, we can re-write above equation as:∫ p(yB̃ | X̃B) dy0 =
∑ G,A,a0 ∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR)
pθ(u0,v |x0,v)p(a0 | ū0,UR) (∫ pθ(z0 |para0(XR,yR)) dz0 ) dUB dZB dū0, (8)
where ∫ pθ(z0 |para0(XR,yR)) dz0 = 1. We further can summarize the equation (8):∫
p(yB̃ | X̃B) dy0 = ∑ G,A ∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR)
pθ(u0,v |x0,v) (∑ a0 p(a0 | ū0,UR) ) dUB dZB dū0. (9)
Hence, we can further marginalize over dū0 as:∫ p(yB̃ | X̃B) dy0 = ∑ G,A ∫ ∏ v∈V
pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR)(∫ pθ(u0,v |x0,v)dū0 ) dUB dZB . (10)
Finally, as ∫ pθ(u0,v |x0,v)dū0 = 1, we can drive equation (7) in the main paper. More precisely, we will have∫ p(yB̃ | X̃B) dy0 =
∑ G,A ∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR) dUB dZB
= p(yB | XB).
In the second care, where {x̄0, y0} ∈ XR\Dx, {x̄0, y0} belongs to the reference set, but not included in the training set. As we discussed in the main paper, the sets XM and Dx will be same across XB and X̃B . Hence we need to marginalize out in terms of all samples in XR\Dx as follow∫
p(yB̃ | X̃B) dy0 =∑ G,A ∫ ∏ v∈V pθ(UB,v |XB,v) p(G,A | UB) pθ(yB ,ZB |A,G,XR) dUB dZBdyi∈XR\Dx
= p(yB | XB). (11)
Therefore, mRNP is consistent under marginalization in both scenarios.

heading: D ADDITIONAL EXPERIMENTS
text: In addition to Figure 4 in the main paper, we further demonstrate that mRNP shows similar behaviour as a GP to encode the inductive biases. In this regression study, we draw 20 samples from U [−4, 4] for the first modality and simulate the target labels as yi = x3i,v1 + , where ∼ N (0, 9
2). Then, we generate samples for the second modality as a nonlinear function xi,v2 = 2xi,v1 + η if xi,v1 ≥ 0 else 0.7xi,v1 , where η ∼ N (0, 0.12). Figure S1 shows the results of both modalities. Similar to GP, mRNP shows the tendency to quickly move towards a flat prediction outside the areas with observed samples.
For this regression task, we train mRNP for a maximum of 1500 epochs using Adam (Kingma & Ba, 2014) with a learning rate of 0.001. We would also like to point out that we separately train an RBF kernel GP for each modality.

heading: E MOG VS POE
text: As we stated overconfident predictions by one expert in PoE can be detrimental to the whole model. By contrast, MoG does not suffer from potentially overconfident graph, since it effectively takes a vote amongst the graphs, and spreads its density over all the individual graphs. Let us assume that we have two input modalities, and the prediction of the first expert (graph) is noisy and overconfident. We would like to compare the effect of this on two modeling perspectives.
Modality I Modality II
Figure S1: Predictive distributions for the second toy multi-modal regression task based on mRNP (top) and GP (bottom).
Effect of overconfident predictions of expert in PoE model. PoE assumes p(z |x1,x2) = p(z |x1)× p(z |x2). Considering we have normal distribution in both experts, i.e. p(z |x1) = N (µ0, σ0) and p(z |x2) = N (µ,σ1). Therefore, p(z |x1,x2) = N (µz, σz), where µz = µ0σ 2 1+µ1σ 2 0
σ20+σ 2 1
and σ2z = σ20σ 2 1
σ20+σ 2 1 . Let us consider that the first expert is overconfident, i.e. σ0 → 0. As a result, p(z |x1,x2) = N (µ0, σ20). This means that the second expert does not matter in this situation and if the first expert be potentially overconfident due to difference in complexity of input modalities or initialisation conditions, second one does not have a chance to correct it.
Effect of overconfident learning of graph in MoG model. MoG assumes p(z |A1,A2) = p(z |αA1 + (1− α)A2), where A1 and A2 are the graphs in two view.
According to equation (4), p(zi |αA1 + (1 − α)A2) = ∏Dz k=1 p(zik |αA1 + (1 −
α)A2) = ∏Dz k=1N (µzik, σzik) , where µzik = ∑ j (αA1,ij + (1− α)A2,ij)ϕ prior,µ k and σ z ik =
exp (∑
j (αA1,ij + (1− α)A2,ij)ϕ prior,σ k ) . If we assume A1 is completely noisy, still A2 is
able to correct it. In addition, we can incorporate prior knowledge through α. If we have a priori information that the first view is noisy, we can consider α (1−α) that leads to reduce the negative effect of the first view.

Review summaries:
This paper proposed to use neural processes for supervised multi-modal learning, which have the abilities to estimate the uncertainty of prediction and to handle missing modalities. Specifically, a directed acyclic graph is learned for each modality in the neural process, which are then used to construct a mixture-of-graphs (MoG) to sidestep the modalities missing problem. The introduction of inducing points makes the predictive distribution tractable. Experimentally, on the label prediction and uncertainty estimation tasks, the proposed model performs well compared to the recent multimodal learning methods.

This paper tries to use Bayesian relational generative model for scalable multi-modal learning. They propose a class of stochastic processes that learns a graph of dependencies between samples across multi-modal data types through adopting priors over the relational structure of the given data modalities. The so-called mRNP method can address the limitations in joint posterior approximation. 

The paper proposes a Bayesian model that learns a predictive distribution given multimodal data with labels. The novelty of the proposed method is that it learns a dependency graph between samples from different modalities and introduces a new mixture-of-graphs (MoG) method to aggregate different dependency graphs.  Experimental results on several datasets demonstrate that the proposed method, which is fully supervised, improves the predictive performance compared to two unsupervised approaches and one supervised baseline. Qualitative results corroborate the quality of the learned representations, which can be visualized as relational graphs, whereas the baselines require additional dimensionality reduction techniques for visualization.