# A Bayesian Gaussian Process-Based Latent Discriminative Generative Decoder
(LDGD) Model for High-Dimensional Data

Navid Ziaei<EMAIL_ADDRESS>
Behzad Nazari<EMAIL_ADDRESS>
Department of Electrical & Computer Engineering
Isfahan University of Technology
Uri T. Eden<EMAIL_ADDRESS>
Department of Mathematics & Statistics
Boston University
Alik<EMAIL_ADDRESS>
Department of Psychiatry & Behavioral Sciences
University of Minnesota
Ali Yousefi<EMAIL_ADDRESS>
Department of Computer Science
Worcester Polytechnic University

###### Abstract

Extracting meaningful information from high-dimensional data poses a
formidable modeling challenge, particularly when the data is obscured by noise
or represented through different modalities. This research proposes a novel
non-parametric modeling approach, leveraging the Gaussian process (GP), to
characterize high-dimensional data by mapping it to a latent low-dimensional
manifold. This model, named the latent discriminative generative decoder
(LDGD), employs both the data and associated labels in the manifold discovery
process. We derive a Bayesian solution to infer the latent variables, allowing
LDGD to effectively capture inherent stochasticity in the data. We demonstrate
applications of LDGD on both synthetic and benchmark datasets. Not only does
LDGD infer the manifold accurately, but its accuracy in predicting data
points’ labels surpasses state-of-the-art approaches. In the development of
LDGD, we have incorporated inducing points to reduce the computational
complexity of Gaussian processes for large datasets, enabling batch training
for enhanced efficient processing and scalability. Additionally, we show that
LDGD can robustly infer manifold and precisely predict labels for scenarios in
that data size is limited, demonstrating its capability to efficiently
characterize high-dimensional data with limited samples. These collective
attributes highlight the importance of developing non-parametric modeling
approaches to analyze high-dimensional data.

Keywords: High-dimensional data analysis, Gaussian process models, Latent
variable, Variational inference, Decoding

##  1 Introduction

Dealing with high-dimensional data presents significant challenges, often
resulting in analytical complications and increased computational demands.
Dimensionality reduction techniques aim to transform complex, high-dimensional
datasets into more manageable and insightful lower-dimensional representations
(Zebari et al., 2020). These techniques enhance data understanding and
interpretation across various scientific domains, including but not limited to
neuroscience, finance, and biology (Koh et al., 2023; Ma and Dai, 2011; Ali et
al., 2017; Zhong and Enke, 2017). For instance, in neural data analysis,
dimensionality reduction simplifies the vast amount of data from brain imaging
studies or neural recordings, aiding in the identification of patterns of
activity underlying cognitive processes or disorders (Cunningham and Yu, 2014;
Wang et al., 2014; Heller and David, 2022).

Dimensionality reduction techniques fall into two main categories: parametric
and non-parametric. Non-probabilistic methods include principal component
analysis (PCA) (Bro and Smilde, 2014), linear discriminant analysis (LDA)
(Balakrishnama and Ganapathiraju, 1998), and more recent approaches including
t-distributed stochastic neighbor embedding (t-SNE) (Cieslak et al., 2020) and
uniform manifold approximation and projection (UMAP) (McInnes et al., 2018).
These methods are widely used for data simplification and visualization.
Despite their significance in advancing science and their popularity, they
often lack the ability to provide a sophisticated understanding, such as
stochasticity in the inferred projections, and are highly sensitive to noise
and outlier data points. On the other hand, probabilistic models, including
probabilistic PCA (PPCA) (Tipping and Bishop, 1999), variational auto-encoders
(VAE) (Kingma et al., 2019), and Gaussian process latent variable model
(GPLVM) (Lawrence and Hyvärinen, 2005) are gaining prominence as dominant
approaches in data analysis. The GPLVM and PPCA address major shortcomings of
non-probabilistic methods, as their inference of low-dimension representation
is probabilistic rather than deterministic. GPLVM is built upon PPCA, with the
added advantage of utilizing non-linear mappings to find lower-dimensional
representations. Despite significant progress in these models, challenges
persist, and there are opportunities to enhance their ability to characterize
complex and high-dimensional data. A key area for improvement is the
integration of label and category information, which these models often
overlook due to their predominantly unsupervised nature. An important
challenge facing Gaussian process-based models like GPLVM is their excessive
computational cost, especially as the data size increases. Our research seeks
to tackle these computational challenges and employ label and category
information to enhance the scalability and utility of models like GPLVM.

To address the high computational cost associated with non-parametric models,
several approaches have been developed. The Nyström approximation (Williams
and Seeger, 2000) offered an initial method to reduce computational
requirements by sampling a subset of the data to approximate the Gaussian
process covariance matrix. This was followed by the introduction of the
pseudo-inputs concept by Snelson and Ghahramani (2005), which parameterizes
its covariance using the locations of $M$ pseudo-input points, determined
optimally through gradient-based optimization. Building upon these
foundations, Titsias (2009) introduced variational inducing points, utilizing
variational inference for the optimal selection of inducing points.
Furthermore, stochastic gradient descent-based variational inducing points,
developed by Hensman et al. (2013), employed stochastic optimization to
further enhance scalability for large datasets. These advancements lay the
foundation for further enhancements in this category of models.

Alongside advancements in enhancing the scalability of GPLVMs, there is
growing research aimed at integrating labels or categories into the
dimensionality reduction process. This integration allows for the inference of
both the underlying structure and its relevance within the mapping.
Incorporating the labels of data points into the dimensionality reduction
process has been explored in approaches such as supervised PPCA (SPPCA) (Yu et
al., 2006). This approach employs two linear mappings: one mapping from the
latent space to the input space, representing the original high-dimensional
data, and another mapping from the latent space to the output space,
corresponding to the data points’ label or category or other supervised
information. This concept was extended to nonlinear dimensionality reduction
techniques, especially with the GPLVM approach. The discriminative GPLVM
(D-GPLVM) applies principles of generalized discriminant analysis (GDA)
(Baudat and Anouar, 2000) to modify latent positions within the data, aiming
to bring data points of the same class closer and distance those of differing
classes. Here, latent positions denote the model’s estimated, unobserved
variables that capture the underlying structure of the data. Other approaches,
such as the supervised GPLVM and shared GPLVM (SGPLVM) (Gao et al., 2010; Ek
and Lawrence, 2009; Jiang et al., 2012), employ two distinct GPLVMs to uncover
nonlinear relationships among latent variables, data, and labels. The
integration of SGPLVM with SPPCA, also known as supervised latent linear GPLVM
(SLLGPLVM), was presented by Jiang et al. (2012). This approach employs two
mapping functions to link the input and output spaces via a latent variable
space. To minimize computational costs, SLLGPLVM employs a linear mapping from
latent variables to high-dimensional data, similar to the method used by Zhang
et al. (2017). While this modification reduces the computational cost, it
sacrifices the flexibility to capture complex nonlinear relationships inherent
in the data, potentially limiting the model’s ability to properly model highly
intricate patterns found in more complex datasets. A further extension to
SGPLVM is the supervised GPLVM based on a Gaussian mixture model (SGPLVM-GMM)
(Zhang et al., 2017). SGPLVM-GMM assumes that latent variables follow a
Gaussian mixture distribution, where mixture components are conditioned on
class labels. Thus, the model simultaneously learns to reduce data dimension
and classify data samples by predicting class probabilities. These models,
including SGPLVM and its extensions such as SLLGPLVM and SGPLVM-GMM, have
successfully utilized latent variables and class labels to enhance feature
extraction and data classification. However, they primarily rely on point
estimates of the latent process, thereby increasing the risk of overfitting,
especially in scenarios where linear relationships do not fully capture data
structure.

Although these models have made progress in incorporating training labels into
latent process inferencing, a gap remains in their ability to manage the
complexities introduced by noise and limited sample sizes. To address these
challenges, our research aims to develop a fully Bayesian dimensionality
reduction solution that integrates the samples’ labels. Another challenge with
previous GP-based models is scalability, which we address in this work through
variational inducing points. The distinctiveness of our model lies in
employing a Bayesian approach for both training and inference while assuming a
prior over latent space. The Bayesian approach helps us better manage data
complexity introduced by noise and missing points; it also lets the framework
be applied to datasets that have a limited sample size. Our approach, which
can be considered an extension of SGPLVM, presents a robust and scalable
method for dimensionality reduction as we incorporate the inducing points.
Through multiple examples, we demonstrate that our proposed method enhances
the feature extraction process by mapping data points to a lower dimension and
ensures the retention of essential data patterns. We illustrate that our
proposed approach can be used as a robust decoder model to draw low-
dimensional representations and labels for new data points. In deriving the
model, we use Gaussian processes (GPs) and Bayesian inference, turning the
model into a non-parametric approach. As such, the solution prevents
overfitting and is best suited for datasets with a limited number of samples.
Our proposed approach can be readily applied to datasets from various fields.
As such, we focus on highlighting the framework application and evaluating its
performance through synthetic and benchmark datasets.

In this research, we introduce our method which is called Latent
Discriminative Generative Decoder (or, briefly LDGD) model. In LDGD, a low-
dimensional random variable defines the underlying generative process
responsible for generating both the continuous values and label data. This
random variable defines the covariance structures of two Gaussian processes
(GPs)—one for the continuous measurement and another for the label. We then
employ a Bayesian inference solution integrated with the variational inducing
points approach to train the model. The usage of inducing points enables
efficient management of computational costs during training and inference.
Specifically, we draw two sets of inducing points, one set for each GP. With
two sets of inducing points, we show that LDGC achieves a balanced prediction
accuracy and data generation capability. Additionally, we will introduce a
batch training pipeline for LDGD, which significantly aids in faster training
and enhances the scalability of the model for larger datasets. LDGD also
serves as an adaptive feature extraction and classification solution, allowing
for the prediction of data point labels from continuous values. Moreover, it
functions as a generative model capable of producing data in high-dimensional
space. To facilitate further research and development, the code for LDGD
implementation is openly available at our GitHub repository:
https://github.com/Navid-Ziaei/LDGD. While some of these attributes are shared
with variational auto-encoders, we argue that this framework represents a
substantial advancement over variational auto-encoders (VAE).

A key enhancement of LDGD over VAEs is its fully Bayesian inference process,
which quantifies prediction uncertainty as the posterior’s variance. Another
notable advancement is LDGD’s ability to partially optimize the dimensions of
latent variables during the training phase, enhancing both feature extraction
and model interpretability. Additionally, LDGD’s Bayesian framework inherently
protects against overfitting, making it well-suited for complex datasets with
limited sample sizes. These advantages collectively position LDGD not just as
a mere alternative but as a significant advancement over traditional encoder-
decoder models, especially in high-dimensional data analysis.

The subsequent sections will explore the technical and applied aspects of the
LDGD modeling framework. We first delve into the intricacies of the LDGD
framework, covering its formulation, model training, and inference procedures
in detail. Additionally, we introduce a variant of LDGD, termed ”fast LDGD”,
which integrates a neural network into the LDGD framework, enhancing its
computational efficiency and predictive power. To assess the performance and
attributes of our proposed model, we have conducted evaluations of LDGD on a
synthetic dataset. Following this, we demonstrate its applicability across
several benchmark datasets, including Oil Flow, Iris, and MNIST. Subsequently,
we compare the distinctive features of LDGD via other methods, such as SLL-
GPLVM and SGPLVM, to highlight its unique advantages and improvements. The
final part of our analysis is a comparison between LDGD and variational auto-
encoders. This thorough review highlights the significant contributions of the
LDGD model to machine learning, especially in efficiently reducing dimensions
and accurately decoding information.

##  2 Materials and Methods

###  2.1 Gaussian process Regression

Gaussian processes (GP) offer a powerful non-parametric approach for
regression and classification tasks (Seeger, 2004). A Gaussian process is a
collection of random variables, with the property that any finite subset of
these variables has a joint Gaussian distribution. GP defines a distribution
over a function defined by:

$p(\mathbf{f}\mid\mathbf{X})=\mathcal{N}(\bm{\mu},\mathbf{K}_{NN}),$

where $\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{N}]$ are input vectors
with $\mathbf{x}_{i}\in\mathbb{R}^{d}$,
$\mathbf{f}=[f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{N})]$ represents latent
variables, called function values in GP context, at the input points with mean
$\bm{\mu}=[m(\mathbf{x}_{1}),\ldots,m(\mathbf{x}_{N})]$ and a $N\times N$
covariance matrix $\mathbf{K}_{NN}$. Each element in the $i$th row and $j$th
column is given by $K_{ij}=k_{\theta}(\mathbf{x}_{i},\mathbf{x}_{j})$ for
$i,j\in\\{1,\ldots,N\\}$, where $k$ is a positive definite kernel constructing
elements of the covariance matrix $\mathbf{K}_{NN}$ and $\theta$ defines the
kernel’s free parameters. In practice, with no prior observation, the mean of
the function is assumed to be zero, i.e., $m(\mathbf{x})=\mathbf{0}$.

Consider a dataset $d=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{N}$, where
$y_{i}\in\mathbb{R}$ represents the noisy observations at each input location
$\mathbf{x}_{i}$ within the domain of interest, encapsulating the underlying
process we aim to model using a Gaussian process. The generation of $y_{i}$ is
modeled as follows:

$y_{i}=f(\mathbf{x}_{i})+\epsilon,$

where $\epsilon$ is considered as Gaussian noise, defined by
$\epsilon\sim\mathcal{N}(0,\sigma^{2}_{y})$. The likelihood of observing
$y_{i}$ given $\mathbf{x}_{i}$ is defined by

$p(y_{i}\mid\mathbf{x}_{i})=\mathcal{N}(y_{i}\mid
f(\mathbf{x}_{i}),\sigma^{2}_{y}).$ (1)

Furthermore, given the GP model, the joint distribution of
$\mathbf{y}=[y_{1},\dots,y_{N}]$ given $\mathbf{X}$ is defined by

$p(\mathbf{y}\mid\mathbf{X})=\mathcal{N}(\mathbf{y}\mid\mathbf{0},\mathbf{K}_{NN}+\sigma^{2}_{y}\mathbb{I}_{N}).$

To train the model, we need to estimate the kernel parameters ($\theta$),
which can be found by solving

$\arg\max_{\theta}\log p(\mathbf{y}\mid\mathbf{X}).$

For a set of test points
$\mathbf{X}^{*}\equiv[\mathbf{x}^{*}_{1},\ldots,\mathbf{x}^{*}_{N^{*}}]$, we
are interested in finding the predictive distribution
$p(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*})=\mathcal{N}(\bm{\mu}^{*},\bm{\Sigma}^{*})$,
where the mean and covariance are defined by

$\displaystyle\bm{\mu}^{*}$
$\displaystyle=\mathbf{K}_{*N}\mathbf{K}_{NN}^{-1}\mathbf{y},$
$\displaystyle\bm{\Sigma}^{*}$
$\displaystyle=\mathbf{K}_{**}-\mathbf{K}_{*N}(\mathbf{K}_{NN}+\sigma^{2}_{y}\mathbb{I}_{N})^{-1}\mathbf{K}_{N*}.$

Here, $\mathbf{K}_{**}$ is the covariance matrix between the test points, and
$\mathbf{K}_{*N}$ is the cross-covariance matrix between the test points and
the training points.

The formulation requires the inversion of an $N\times N$ matrix to derive the
predictive distribution, resulting in a computational complexity of
$\mathcal{O}(N^{3})$. This computation becomes expensive in datasets even with
a moderate number of data points. A solution that addresses this challenge is
based on variational inducing points proposed by Titsias (2009). In this
approach, $M$ data points $\mathbf{Z}\equiv\\{\mathbf{z}_{i}\\}_{i=1}^{M}$,
where $\mathbf{z}_{i}\in\mathbb{R}^{d}$, is introduced in the input space.
These $M$ data points, where $M\ll N$, replace the $N$ original data points
used in deriving the covariance matrix. Inducing points $\mathbf{Z}$ can be
chosen as fixed points or model-free parameters (Titsias, 2009). This research
considers $\mathbf{Z}$ as free parameters with random initial values. We
further assume for each inducing point, there is a corresponding inducing
variable $u_{i}$. Induced variables
$\mathbf{u}\equiv\left[u_{1},...,u_{M}\right]$ follow the distribution
$p(\mathbf{u}\mid\mathbf{Z})=\mathcal{N}(\mathbf{u}\mid\mathbf{0},\mathbf{K}_{MM})$,
where $\mathbf{K}_{MM}$ is an $M\times M$ matrix. This matrix is generated
using the same kernel function $k_{\theta}(z_{i},z_{j})$ as in the GP model,
applied to the inducing points. The conditional distribution of $\mathbf{f}$
given $\mathbf{u}$ is:

$p(\mathbf{f}\mid\mathbf{u};\mathbf{Z},\mathbf{X})=\mathcal{N}(\mathbf{f}\mid\bm{\mu_{f}},\bm{\Sigma_{f}}),$
(2)

where $\bm{\mu_{f}}=\mathbf{K}_{NM}\mathbf{K}_{MM}^{-1}\mathbf{u}$ and
$\bm{\Sigma_{f}}=\mathbf{K}_{NN}-\mathbf{K}_{NM}\mathbf{K}_{MM}^{-1}\mathbf{K}_{MN}$.
Using the inducing points, we can draw predictions for new points with a much
smaller dataset and significantly reduce the computational cost for the
covariance inverse. However, to use $\mathbf{Z}$ for the prediction step, we
first need to estimate these points. With the inducing points, the free
parameters of the GP model increase to include both $\theta$ and $\mathbf{Z}$.
We find the maximum likelihood estimate of $\theta$ and $\mathbf{Z}$ through
the marginal distribution of $\mathbf{y}$, defined by

$\displaystyle p(\mathbf{y}\mid\mathbf{X};\mathbf{Z},\theta)$
$\displaystyle=\int
p(\mathbf{y},\mathbf{f},\mathbf{u}\mid\mathbf{X};\mathbf{Z},\theta)\,d\mathbf{f}\,d\mathbf{u}$
$\displaystyle=\int
p(\mathbf{y}\mid\mathbf{f};\theta)p(\mathbf{f}\mid\mathbf{u},\mathbf{X};\mathbf{Z},\theta)p(\mathbf{u};\mathbf{Z},\theta)\,d\mathbf{f}\,d\mathbf{u},$

where $\mathbf{f}$ and $\mathbf{u}$ are set of latent variables of length $N$
and $M$, respectively. To determine $\theta$ and $\mathbf{Z}$ without needing
to invert the large matrix $\mathbf{K}_{NN}$, we employ a variational approach
by introducing the variational distribution $q(\mathbf{f},\mathbf{u})$. We
define $q(\mathbf{f},\mathbf{u})$ as
$q(\mathbf{u})q(\mathbf{f}\mid\mathbf{u})$, where
$q(\mathbf{u})=\prod_{j=1}^{M}\mathcal{N}(u_{j};m_{j},s_{j})$ and
$q(\mathbf{f}\mid\mathbf{u})=p(\mathbf{f}\mid\mathbf{u},\mathbf{X};\mathbf{Z},\theta)$,
as already defined in Equation (2). With the $q$ distribution, we can derive a
lower bound on the marginal distribution of $\mathbf{y}$ using Jensen
inequality, which becomes:

$\displaystyle\log p(\mathbf{y}\mid\mathbf{X};\mathbf{Z},\theta)$
$\displaystyle\geq\mathbb{E}_{q(\mathbf{f},\mathbf{u})}[\log
p(\mathbf{y}\mid\mathbf{f};\theta)]-\text{KL}[q(\mathbf{u})\parallel
p(\mathbf{u};\mathbf{Z},\theta)],$ (3)

where we find $\mathbf{Z}$ and $\theta$, along with $m_{j},s_{j}$ for
$j=1,\ldots,M$ to maximize the lower bound. The first term on the right-hand
side of the Equation 3 encourages the model to fit the data well, while the
second term regularizes the model by penalizing deviations of the variational
distribution from the prior. It is worth noting that the inversion of the
$\mathbf{K}_{NN}$ matrix is no longer needed in prediction and parameter
estimation steps.

For $N_{t}$ test points
$\mathbf{X}^{*}\equiv\left[x^{*}_{1},...,x^{*}_{N_{t}}\right]$, the predictive
distribution is defined by
$p(\mathbf{f}^{*}\mid\mathbf{X},\mathbf{y},\mathbf{X}^{*})$. Given the
inducing points, we can write this distribution as:

$p(\mathbf{f}^{*}\mid\mathbf{X},\mathbf{y},\mathbf{X}^{*})=\int
p(\mathbf{f}^{*},\mathbf{u}\mid\mathbf{X},\mathbf{y},\mathbf{X}^{*})\,d\mathbf{u}.$

Note that in the equation, for clarity, we dropped its parameters - i.e.,
$\mathbf{Z}$ and $\theta$. We estimate the exact posterior
$p(f^{*},\mathbf{u}\mid\mathbf{X},\mathbf{y},\mathbf{X}^{*})$ with the updated
variational distribution $q(f^{*},\mathbf{u})$. Thus, the predictive
distribution can be approximated as:

$\displaystyle p(\mathbf{f}^{*}\mid\mathbf{X},\mathbf{y},\mathbf{X}^{*})$
$\displaystyle\approx\int q(\mathbf{f}^{*},\mathbf{u})\,d\mathbf{u}=\int
p(\mathbf{f}^{*}\mid\mathbf{u})q(\mathbf{u})\,d\mathbf{u}.$ (4)

We have already derived these distributions during the model training. They
are defined by:

$\displaystyle q(\mathbf{u})$
$\displaystyle\sim\mathcal{N}(\mathbf{u}\mid\mathbf{m},\mathbf{S}),$
$\displaystyle p(\mathbf{f}^{*}\mid\mathbf{u})$
$\displaystyle\sim\mathcal{N}\left(\mathbf{f}^{*}\mid\mathbf{K}_{*M}\mathbf{K}_{MM}^{-1}\mathbf{u},\mathbf{K}_{**}-\mathbf{K}_{*M}\mathbf{K}_{MM}^{-1}\mathbf{K}_{M*}\right).$

where $\mathbf{m}=[m_{1},\ldots,m_{M}]$ and $\mathbf{S}$ is a $M\times M$
diagonal matrix. The diagonal elements of $\mathbf{S}$ are
$\\{s_{i}\\}_{i=1}^{M}$, which are the variational free parameters that are
optimized during training. The approximated integral in Equation (4) has a
closed form, which is a Gaussian with the following mean and covariance:

$\displaystyle\bm{\mu}^{*}$
$\displaystyle=\mathbf{K}_{*M}\mathbf{K}_{MM}^{-1}\mathbf{m},$
$\displaystyle\bm{\Sigma}^{*}$
$\displaystyle=\mathbf{K}_{**}-\mathbf{K}_{*M}\mathbf{K}_{MM}^{-1}(\mathbf{K}_{MM}-S)\mathbf{K}_{MM}^{-1}\mathbf{K}_{M*}.$

###  2.2 Gaussian Process Latent Variable Model

Gaussian process latent variable models emerge from the necessity of
understanding high-dimensional data by projecting it into a simpler and much
lower-dimensional latent space. At its core, a GPLVM is built upon the
hypothesis that even though we observe data $\mathbf{Y}\in\mathbb{R}^{N\times
D}$, there exists a latent representation $\mathbf{X}\in\mathbb{R}^{N\times
Q}$, with $Q\ll D$, that captures the essential features or structures present
in the high-dimensional data. This idea was inspired by PPCA, evolving into a
nonlinear dimensionality reduction method that extends the capabilities of its
precursor (Tipping and Bishop, 1999). The relationship between GPLVM and PPCA
is explained in Appendix A. In GPLVM, the conditional distribution of
$\mathbf{Y}$ given $\mathbf{X}$ is assumed to be:

$p(\mathbf{Y}\mid\mathbf{X})=\prod_{d=1}^{D}p(\mathbf{y}_{:,d}\mid\mathbf{X})\text{,}$

where $\mathbf{y}_{:,d}$ is $d$th column of $\mathbf{Y}$. Each element of
$\mathbf{y}_{:,d}$ is a noisy realization of the function
$\mathbf{f}_{d}\in\mathbb{R}^{N}$ at location $\mathbf{x}_{i}$ defined by

$y_{i,d}=f_{d}(\mathbf{x}_{i})+\epsilon_{d},$

where $\epsilon_{d}\sim\mathcal{N}(0,\sigma^{2}_{d})$ and
$\mathbf{f}_{d}\equiv\left[f_{d}(\mathbf{x}_{1}),...,f_{d}(\mathbf{x}_{N})\right]$
is function values with GP prior. This leads to the following conditional
distribution:

$p(\mathbf{y}_{:,d}\mid\mathbf{X})=\mathcal{N}(\mathbf{y}_{:,d}\mid\mathbf{0},\mathbf{K}_{NN}+\sigma_{d}^{2}\mathbb{I}_{N}),$

with $\mathbf{K}_{NN}$ being an $N\times N$ covariance matrix defined using a
kernel function as a function of the latent space $\mathbf{X}$.

We use the ARD kernel (Seeger, 2004) for the kernel. This kernel is generally
used for determining the relevance of dimensions in $\mathbf{X}$ by tuning
${\alpha}_{q}$. A typical ARD kernel is defined by:

$k(\mathbf{\mathbf{x}},\mathbf{\mathbf{x}^{\prime}})=\sigma_{f}^{2}\exp\left(-\frac{1}{2}\sum_{q=1}^{Q}\alpha_{q}(x_{q}-x^{\prime}_{q})^{2}\right),$

where $\mathbf{x}$ and $\mathbf{x^{\prime}}$ represent two input vectors in
the latent space. The term $\sigma_{f}^{2}$ denotes the variance parameter of
the kernel, dictating the scale of the output space. Each $\alpha_{q}$ is an
inverse length-scale parameter associated with the $q$th dimension of the
latent space, which influences how variations in that particular dimension
affect the kernel computation. These length-scale parameters allow the kernel
to adjust the relevance of the feature dimensions.

The GPLVM can be interpreted as a variant of a multi-output Gaussian process.
Every dimension of the observed data, $\mathbf{Y}$, is a separate realization
of a Gaussian process, all derived by a common latent structure (Liu et al.,
2018). To fit GPLVM to a dataset, we require finding a suitable representation
of the latent variables that best explain the observed data. This can be done
by estimating ARD kernel-free parameters. Three primary methods used for
finding suitable representation include point estimate, Maximum A Posteriori
(MAP) (Lawrence and Hyvärinen, 2005), and a fully Bayesian approach (Titsias
and Lawrence, 2010). Whereas the point estimate method can be more susceptible
to overfitting, the Bayesian approach typically offers a better generalization
to data due to its integration over uncertainties. The payoff is that the
Bayesian approach is computationally more expensive and often requires
sophisticated inference techniques such as Markov Chain Monte Carlo (MCMC)
(Gilks et al., 1995) or variational inferencing (Blei et al., 2017). This
research focuses on Bayesian methods to derive a more robust inference. In the
Bayesian framework, it is imperative to specify a prior distribution for the
latent variables. We assume the prior on the latent variables are Gaussian
defined by:

$p(\mathbf{X})=\prod_{n=1}^{N}\mathcal{N}(\mathbf{x}_{n}\mid\mathbf{0},\mathbb{I}_{Q}),$

where $\mathbf{x}_{n}$ denotes the $n$th row of the latent variables matrix
$\mathbf{X}$. With this prior, the marginal probability of $\mathbf{Y}$
involves the integral, which is defined by:

$p(\mathbf{Y})=\int{p(\mathbf{Y}\mid\mathbf{X})p(\mathbf{X})d\mathbf{X}}=\int{\prod_{d=1}^{D}p(\mathbf{y}_{:,d}\mid\mathbf{X})p(\mathbf{X})d\mathbf{X}}.$
(5)

This integral corresponds to the expectation of conditional probability of
$\mathbf{Y}$ over $\mathbf{X}$. The integral does not have a closed-form
solution, as $\mathbf{K}_{NN}$ is a function of $\mathbf{X}$. To find the
marginal distribution, we apply Jensen’s inequality, which gives a variational
lower bound on the log of likelihood:

$\log{p(\mathbf{Y})}\geq\sum_{d=1}^{D}{\mathbb{E}_{q(\mathbf{X})}[\log
p(\mathbf{y}_{:,d}\mid\mathbf{X})]}-\text{KL}(q(\mathbf{X})\parallel
p(\mathbf{X}\mid\mathbf{Y})).$

It is necessary to pick a proper $\mathbf{q}$ distribution to achieve a tight
bound here. The number of samples in the training set introduces additional
complexities to the model; thus, we will use the inducing points to derive a
more tractable solution to calculate the bound. Different strategies have been
proposed to calculate the bound, including the solution proposed by Hensman et
al. (2013), which creates a closed-form solution for the bound. Salimbeni and
Deisenroth (2017) employed the variational distribution to simplify the
problem by enabling lower-dimensional sampling conducive to convergence. Here,
we discuss our approach, which integrates the labels or categories attached to
each data point into the formation of our low-dimensional representation,
thereby facilitating the classification of these labels. This has been
accomplished through the doubly variational method, upon which we have
developed the LDGD.

###  2.3 Latent Discriminative Generative Decoder Model

Here, we introduce the LDGD framework. While the idea behind LDGD is the same
as GPVLM, our objective is different. In LDGD, we aim to infer the latent
process that can capture essential features of the observed data and separate
the classes in the data, contrasting with GPLVM, which primarily focuses on
dimensionality reduction without explicitly optimizing for class separability.
Essentially, we are identifying a set of features in the data that can be used
for both data generation and achieving high accuracy in the data class or
category prediction. We note that LDGD inference will differ from those
derived by generative auto-encoders, as the inferred latent process needs to
capture those essential features of the data that are expressive of labels or
categories of data. The following section will delve into the framework
definition and its different modeling steps, including training and
prediction.

####  2.3.1 Observed Data

In LDGD, we assume the observed data consists of continuous measures,
represented as $\mathbf{Y}^{r}\in\mathbb{R}^{N\times D}$, and corresponding
classes of data or categories, represented as
$\mathbf{Y}^{c}\in\mathbb{R}^{N\times K}$. $N$ denotes the number of samples,
$D$ is the dimension of continuous measures per sample, and $K$ is the number
of distinct categories or classes per sample. Each row of $\mathbf{Y}^{r}$,
denoted as $\mathbf{y}^{r}_{i}$, corresponds to the continuous measurement for
the $i$th sample. Each column, denoted as $\mathbf{y}^{r}_{:,d}$, represents
the data for the $d$th dimension across all samples. Similarly, each row of
$\mathbf{Y}^{c}$, denoted as $\mathbf{y}^{c}_{i}$, corresponds to the one-hot
encoded class label for each sample.

####  2.3.2 Latent Variables

LDGD incorporates latent variables $\mathbf{X}\in\mathbb{R}^{N\times Q}$, to
capture the essential features driving both $\mathbf{Y}^{r}$ and
$\mathbf{Y}^{c}$. Here, $Q$ represents the latent space’s dimensionality,
which is generally much smaller than the dimension of the observed data. We
assume a Gaussian prior for these latent variables, expressed as:

$p(\mathbf{X})=\prod_{i=1}^{N}\mathcal{N}(\mathbf{x}_{i}\mid\mathbf{0},\mathbb{I}_{Q}),$
(6)

where $\mathbf{0}$ is a zero vector and $\mathbb{I}_{Q}$ is the identity
matrix of size $Q\times Q$.

####  2.3.3 Gaussian process Priors

To model the relationships between the latent space and the observed data, we
impose two Gaussian process (GP) priors on function values $\mathbf{F}^{r}$
and $\mathbf{F}^{c}$. Specifically, we define GP priors over
$\mathbf{F}^{r}\in\mathbb{R}^{N\times D}$, for the continuous measurements,
and $\mathbf{F}^{c}\in\mathbb{R}^{N\times K}$ for the categorical measurements
or labels. These prior distributions are defined by:

$\displaystyle p(\mathbf{F}^{c}\mid\mathbf{X})$
$\displaystyle=\prod_{k=1}^{K}\mathcal{N}(\mathbf{f}_{:,k}^{c}\mid\mathbf{0},\mathbf{K}_{NN}^{c})$
(7) $\displaystyle p(\mathbf{F}^{r}\mid\mathbf{X})$
$\displaystyle=\prod_{d=1}^{D}\mathcal{N}(\mathbf{f}_{:,d}^{r}\mid\mathbf{0},\mathbf{K}_{NN}^{r}),$
(8)

where $\mathbf{K}_{NN}^{c}$ and $\mathbf{K}_{NN}^{r}$ are the $N\times N$
covariance matrices for the GPs corresponding to class labels and observed
measurements, respectively. The element in the $i$th row and $j$th column of
these matrices are defined by kernel functions
$k^{c}_{i,j}=k_{\psi}(\mathbf{x}_{i},\mathbf{x}_{j})$ and
$k^{r}_{i,j}=k_{\theta}(\mathbf{x}_{i},\mathbf{x}_{j})$, with $\theta$ and
$\psi$ as the kernel parameters. Furthermore, $\mathbf{f}_{:,k}^{c}$
represents the $k$th column of $\mathbf{F}^{c}$, and $\mathbf{f}_{:,d}^{r}$
denotes the $d$th column of $\mathbf{F}^{r}$. The independence assumption
embedded in Equation (7) and Equation (8) is a direct consequence of the dual
formulation in PPCA. A detailed explanation of duality is provided in the
Appendix A.

####  2.3.4 Observed Data Model

LDGD delineates the relationship between the features and the observed data.
This relationship is modeled by:

$\displaystyle p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})$
$\displaystyle=\prod_{i=1}^{N}\prod_{d=1}^{D}\mathcal{N}(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}),\sigma_{d}^{2})$ (9) $\displaystyle
p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})$
$\displaystyle=\prod_{i=1}^{N}\prod_{k=1}^{K}\text{Bernoulli}(y_{i,k}^{c}\mid
g(f_{k}^{c}(\mathbf{x}_{i}))),$ (10)

where $\sigma_{d}^{2}$ is the variance of the Gaussian noise and $g(\cdot)$ is
a squash function, such as the Sigmoid or Probit, ensuring the output
represents a probability. For the $i$th element in the $d$th column of
$\mathbf{F}^{r}$, we define $f^{r}_{i,d}\triangleq f^{r}_{d}(\mathbf{x}_{i})$.
Similarly, for the $i$th element in the $k$th column of $\mathbf{F}^{c}$, we
define $f^{c}_{i,k}\triangleq f^{c}_{k}(\mathbf{x}_{i})$. Figure 1 shows the
graphical model behind the LDGD framework.

$\mathbf{X}$$\mathbf{F}^{c}_{:,k}$$\mathbf{F}^{r}_{:,d}$$\mathbf{y}^{c}_{:,k}$$\mathbf{y}^{r}_{:,d}$$k=1,\ldots,K$$d=1,\ldots,D$
(a)

$\mathbf{x}_{i}$$\mathbf{F}^{c}_{:,k}$$\mathbf{F}^{r}_{:,d}$$y^{c}_{i,k}$$y^{r}_{i,d}$$\mathbf{u}^{c}_{:,k}$$\mathbf{u}^{r}_{:,d}$$\mathbf{Z}^{r}$$\mathbf{Z}^{c}$$\mathbf{\mu_{i}}$$\mathbf{s}_{i}$$\mathbf{m}^{c}_{k}$$\mathbf{S}^{c}_{k}$$\mathbf{m}^{r}_{d}$$\mathbf{S}^{r}_{d}$$k=1,\ldots,K$$d=1,\ldots,D$$i=1,\ldots,N$
(b)

Figure 1: Graphical Models Depicting LDGD. (a) Exact inference, (b)
Variational inference

####  2.3.5 Marginal Likelihood Lower Bound

Equation sets (6) to (10) define the LDGD model. We need to find the marginal
distribution of $Y^{r}$ and $Y^{c}$ to fit the model to the dataset. The
marginal distribution is defined by:

$p(\mathbf{Y}^{r},\mathbf{Y}^{c})=\int
p(\mathbf{Y}^{r},\mathbf{Y}^{c},\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{X})\,d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{X}.$
(11)

Similar to Equation (5), this integral has no closed-form solution.
Additionally, calculating the integral numerically is impractical due to the
high dimensionality of the latent elements over which the integral is taken.
Given there is no solution to find this integral, maximizing the evidence by
adjusting its free parameter becomes impossible. Instead, we derive a lower
bound for the marginal likelihood, similar to what we showed in Equation (3).
To achieve this, we will be using inducing points and variational
distributions. As we show later, our solution will have a modest computational
cost and can be applied to high-dimensional datasets. Besides, the model
training can be done through batches of datasets, making it appropriate for
iterative training.

##### Inducing Points and Variational Distributions

In Section 2.1, we have explained the utilization of inducing points for GP
regression. Similarly, in LDGD, we apply the same concept, with the
distinction that we now need to address both discrete and continuous
observation processes. Attributes of these processes impose different
sensitivities to the choice of inducing points; we pick two sets of inducing
points, one for the discrete and one for the continuous observations. For the
label data in the classification component, we define the inducing points and
variables as

$p(\mathbf{U}^{c}\mid\mathbf{Z}^{c})=\prod_{k=1}^{K}\mathcal{N}(\mathbf{u}_{:,k}^{c}\mid\mathbf{0},k^{c}(\mathbf{Z}^{c},\mathbf{Z}^{c})),$
(12)

where
$\mathbf{U}^{c}\equiv\left[\mathbf{u}_{:,1}^{c},\ldots,\mathbf{u}_{:,K}^{c}\right]$,
$\mathbf{u}_{:,k}^{c}\in\mathbb{R}^{M}$, and
$\mathbf{K}_{M_{c}M_{c}}^{c}\triangleq k^{c}(\mathbf{Z}^{c},\mathbf{Z}^{c})$.
Here, $\mathbf{Z}^{c}$ represents the set with $M_{c}$ inducing points
associated with the $\mathbf{F}^{c}$. Analogously, for the continuous
measurement of the regression component, we have:

$p(\mathbf{U}^{r}\mid\mathbf{Z}^{r})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{u}_{:,d}^{r}\mid\mathbf{0},k^{r}(\mathbf{Z}^{r},\mathbf{Z}^{r})),$
(13)

where
$\mathbf{U}^{r}\equiv\left[\mathbf{u}_{:,1}^{r},\ldots,\mathbf{u}_{:,D}^{r}\right]$,
$\mathbf{u}_{:,d}^{r}\in\mathbb{R}^{M}$, and
$\mathbf{K}_{M_{r}M_{r}}^{r}\triangleq k^{r}(\mathbf{Z}^{r},\mathbf{Z}^{r})$.
Here, we use a new set of $M_{r}$ inducing points, $Z^{r}$, associated with
the $\mathbf{F}^{r}$. The choice of two disjoint sets of inducing points
addresses different sensitivities we have in the regression and
classification.

##### Doubly Stochastic Variational GP and Variational Posterior Distribution.

The true posterior distribution is given by

$p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})=\frac{p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})}{p(\mathbf{Y}^{r},\mathbf{Y}^{c})}.$
(14)

Calculating the denominator involves the computation of Equation (11), which
is intractable. We propose variational distributions as a part of our approach
to approximate the posterior distribution in LDGD. Our methodology
incorporates concepts from the doubly stochastic variational Gaussian process
(DSVGP) outlined by Salimbeni and Deisenroth (2017), and integrates the
scalable variational GP framework for classification proposed by Hensman et
al. (2015). This integration is achieved through the utilization of two
distinct sets of inducing variables, as delineated in Equations (12) and (13),
which are applied respectively for classification and regression tasks.
Training involves maximizing the marginal likelihood (evidence), and our
method provide a means to find a lower bound for the otherwise intractable
marginalization process. To achieve this lower bound, we begin by
marginalizing the joint distribution of LDGD, taking into account the
variational inducing points introduced:

$\displaystyle p(\mathbf{Y}^{r},\mathbf{Y}^{c})=$ $\displaystyle\int
p(\mathbf{Y}^{r},\mathbf{Y}^{c},\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})\,d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}.$
(15)

The joint distribution of LDGD, given the observed data $\mathbf{Y}^{r}$ and
$\mathbf{Y}^{c}$, the function values $\mathbf{F}^{r}$ and $\mathbf{F}^{c}$,
the inducing variables $\mathbf{U}^{r}$ and $\mathbf{U}^{c}$, and the latent
features $\mathbf{X}$, can be factorized as:

$\displaystyle
p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})=$
$\displaystyle
p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{F}^{c},\mathbf{U}^{c}\mid\mathbf{X})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{F}^{r},\mathbf{U}^{r}\mid\mathbf{X})p(\mathbf{X}).$
(16)

Solving this integral is not straightforward. A common technique used in
variational inference is multiplying and dividing the terms inside the
integral by a proposal distribution, which is an approximation of the joint
posterior. The joint posterior
$p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})$
can be factorized into three separate posteriors:

$p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})=p(\mathbf{F}^{c},\mathbf{U}^{c}\mid\mathbf{Y}^{c},\mathbf{X})p(\mathbf{F}^{r},\mathbf{U}^{r}\mid\mathbf{Y}^{r},\mathbf{X})p(\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c}).$
(17)

where $p(\mathbf{F}^{c},\mathbf{U}^{c}\mid\mathbf{Y}^{c},\mathbf{X})$ and
$p(\mathbf{F}^{r},\mathbf{U}^{r}\mid\mathbf{Y}^{r},\mathbf{X})$ are the joint
posterior of the Gaussian process function values and their inducing variables
for regression and classification respectively.
$p(\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})$ shows the posterior over the
latent features. To derive a tractable solution using posterior, we
approximate posterior $p(\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})$ with
variational distribution $q_{\phi}(\mathbf{X})$,
$p(\mathbf{F}^{c},\mathbf{U}^{c}\mid\mathbf{Y}^{c},\mathbf{X})$ with
$q(\mathbf{F}^{c},\mathbf{U}^{c})$ and
$p(\mathbf{F}^{r},\mathbf{U}^{r}\mid\mathbf{Y}^{r},\mathbf{X})$ with
$q(\mathbf{F}^{r},\mathbf{U}^{r})$:

$\displaystyle q_{\phi}(\mathbf{X})$
$\displaystyle=\prod_{i=1}^{N}\mathcal{N}(\mathbf{x}_{i}\mid\bm{\mu}_{i},\mathbf{s}_{i}\mathbb{I}_{Q}),$
(18) $\displaystyle q(\mathbf{F}^{c},\mathbf{U}^{c})$
$\displaystyle=p(\mathbf{F}^{c}\mid\mathbf{U}^{c},\mathbf{X})q_{\lambda}(\mathbf{U}^{c}),\quad
q_{\lambda}(\mathbf{U}^{c})=\prod_{k=1}^{K}\mathcal{N}(\mathbf{u}_{k}^{c}\mid\mathbf{m}_{k}^{c},\mathbf{S}_{k}^{c}),$
(19) $\displaystyle q(\mathbf{F}^{r},\mathbf{U}^{r})$
$\displaystyle=p(\mathbf{F}^{r}\mid\mathbf{U}^{r},\mathbf{X})q_{\gamma}(\mathbf{U}^{r}),\quad
q_{\gamma}(\mathbf{U}^{r})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{u}_{d}^{r}\mid\mathbf{m}_{d}^{r},\mathbf{S}_{d}^{r}),$
(20)

where $\phi=\\{\bm{\mu}_{i},\mathbf{s}_{i}\\}_{i=1}^{N}$,
$\lambda=\\{\mathbf{m}_{k}^{c},\mathbf{S}_{k}^{c}\\}_{k=1}^{K}$,and
$\gamma=\\{\mathbf{m}_{d}^{r},\mathbf{S}_{d}^{r}\\}_{d=1}^{D}$ are the
variational parameters. Using the variational approximations, the joint
posteriors can be expressed as:

$p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})\approx
q(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})\equiv
q(\mathbf{F}^{c},\mathbf{U}^{c})q(\mathbf{F}^{r},\mathbf{U}^{r})q_{\phi}(\mathbf{X}).$
(21)

Employing the approximated posterior, we reformulate the evidence introduced
in Equation (15) and apply Jensen’s inequality on the logarithm of the
evidence to establish the evidence lower-bound (ELBO):

$\displaystyle\log p(\mathbf{Y}^{r},\mathbf{Y}^{c})\geq$ $\displaystyle
E_{q(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})}\left[\log\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}\right]\triangleq\text{ELBO}.$

We can further breakdown this ELBO in 5 terms

ELBO
$\displaystyle\equiv\text{ELL}^{\text{reg}}+\text{ELL}^{\text{cls}}-\text{KL}_{u}^{c}-\text{KL}_{u}^{r}-\text{KL}_{X},$
(22)

where the $\text{ELL}^{reg}$ and $\text{ELL}^{cls}$ are expected log-
likelihood (ELL) for regression and classification paths, respectively, and
are defined by:

$\displaystyle\text{ELL}^{\text{reg}}$ $\displaystyle\equiv
E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{r}\mid\mathbf{U}^{r},\mathbf{X})q_{\gamma}(\mathbf{U}^{r})}\left[\log
p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})\right]\right],$
$\displaystyle\text{ELL}^{\text{cls}}$ $\displaystyle\equiv
E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{c}\mid\mathbf{U}^{c},\mathbf{X})q_{\lambda}(\mathbf{U}^{c})}\left[\log
p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})\right]\right].$

The KL terms are defined as:

$\displaystyle\text{KL}_{u}^{c}\equiv\textit{KL}(q_{\lambda}(\mathbf{U}^{c})\parallel
p(\mathbf{U}^{c})),$
$\displaystyle\text{KL}_{u}^{r}\equiv\textit{KL}(q_{\gamma}(\mathbf{U}^{r})\parallel
p(\mathbf{U}^{r})),$
$\displaystyle\text{KL}_{X}\equiv\textit{KL}(q_{\phi}(\mathbf{X})\parallel
p(\mathbf{X})).$

$\text{ELL}^{reg}$ and $\text{ELL}^{cls}$ contribute to data fitting, whereas
$\text{KL}_{u}^{c}$, $\text{KL}_{u}^{r}$, and $\text{KL}_{X}$ terms serve as
regularizers and help prevent overfitting. Minimizing the KL divergence
between the posterior defined in Equation (17) and its variational estimation
defined in Equation (21) is equivalent to maximizing the ELBO. A detailed
explanation is provided in Appendix B. A closed form for this lower bound is
desirable for optimizing the variational parameters, inducing points, and
kernel parameters. While KL terms have a closed form due to the Gaussian
distribution characteristics, $\text{ELL}^{cls}$ and $\text{ELL}^{reg}$ do not
have a closed form in general. The details of calculating these two terms are
explained in the subsequent sections.

####  2.3.6 Classification Expected Log-Likelihood ($\text{ELL}^{\text{cls}}$)

In this section, we delve into the detailed analysis of
$\text{ELL}^{\text{cls}}$, defined as

$\displaystyle\text{ELL}^{\text{cls}}$
$\displaystyle=E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{c}\mid\mathbf{U}^{c},\mathbf{X})q_{\lambda}(\mathbf{U}^{c})}\left[\log
p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})\right]\right]$ $\displaystyle=\int
q_{\phi}(\mathbf{X})p(\mathbf{F}^{c}\mid\mathbf{U}^{c},\mathbf{X})q_{\lambda}(\mathbf{U}^{c})\sum_{i,k}\log
p(y_{i,k}^{c}\mid
f_{k}^{c}(\mathbf{x}_{i}))\,d\mathbf{F}^{c}\,d\mathbf{U}^{c}\,d\mathbf{X}$
$\displaystyle=\sum_{i,k}\int
q_{\phi}(\mathbf{x}_{i})p(f_{k}^{c}\mid\mathbf{u}_{k}^{c},\mathbf{x}_{i})q_{\lambda}(\mathbf{u}_{k}^{c})\log
p(y_{i,k}^{c}\mid
f_{k}^{c}(\mathbf{x}_{i}))\,df_{k}^{c}\,d\mathbf{u}_{k}^{c}\,d\mathbf{x}_{i}$
$\displaystyle=\sum_{i,k}ELL_{i,k}^{\text{cls}}.$ (23)

We rewrite the $ELL_{i,k}^{\text{cls}}$:

$\displaystyle ELL_{i,k}^{\text{cls}}$ $\displaystyle=\int
q_{\phi}(\mathbf{x}_{i})p(f_{k}^{c}\mid\mathbf{u}_{k}^{c},\mathbf{x}_{i})q_{\lambda}(\mathbf{u}_{k}^{c})\log
p(y_{i,k}^{c}\mid
f_{k}^{c}(\mathbf{x}_{i}))\,df_{k}^{c}\,d\mathbf{u}_{k}^{c}\,d\mathbf{x}_{i}$
$\displaystyle=\int
q_{\phi}(\mathbf{x}_{i})q_{\lambda}(f_{k}^{c}\mid\mathbf{x}_{i})\log
p(y_{i,k}^{c}\mid f_{k}^{c}(\mathbf{x}_{i}))\,df_{k}^{c}\,d\mathbf{x}_{i}$
$\displaystyle=E_{q_{\phi}(\mathbf{x}_{i})q_{\lambda}(f_{k}^{c}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,k}^{c}\mid f_{k}^{c}(\mathbf{x}_{i}))\right],$

where we define the predictive distribution as:

$\displaystyle q_{\lambda}(f_{k}^{c}\mid\mathbf{x}_{i})$
$\displaystyle\triangleq\int p(f_{k}^{c}\mid
u_{k}^{c},\mathbf{x}_{i})q_{\lambda}(\mathbf{u}_{k}^{c})\,d\mathbf{u}_{k}^{c}=\mathcal{N}\left(\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i}),\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})\right)$
$\displaystyle\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i})$
$\displaystyle=\mathbf{k}_{iM_{c}}^{c}(\mathbf{K}_{M_{c}M_{c}}^{c})^{-1}\mathbf{m}_{k}^{c}$
$\displaystyle\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})$
$\displaystyle=k_{ii}^{c}+\mathbf{k}_{iM_{c}}^{c}(\mathbf{K}_{M_{c}M_{c}}^{c})^{-1}(\mathbf{S}_{d}^{c}-\mathbf{K}_{M_{c}M_{c}}^{c})(\mathbf{K}_{M_{c}M_{c}}^{c})^{-1}\mathbf{K}_{M_{c}i}^{c},$

with $\mathbf{k}^{c}_{iM_{c}}=k^{c}(\mathbf{x}_{i},\mathbf{Z}^{c})$ being the
$i$th row of $K^{c}_{NM_{c}}$ and
$k^{c}_{ii}=k^{c}(\mathbf{x}_{i},\mathbf{x}_{i})$. The details are provided in
Appendix C. For ease of use, we define $B(\mathbf{x}_{i})\triangleq
E_{q_{\gamma}(\mathbf{f}_{k}^{c}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,k}^{c}\mid f_{k}^{c}(\mathbf{x}_{i}))\right]$. For the Bernoulli
likelihood with Probit likelihood function, the inner expectation
$B(\mathbf{x}_{i})$ is:

$\displaystyle B(\mathbf{x}_{i})$
$\displaystyle=E_{q_{\gamma}(\mathbf{f}_{k}^{c}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,k}^{c}\mid f_{k}^{c}(\mathbf{x}_{i}))\right]$ $\displaystyle=\int
q_{\gamma}(\mathbf{f}_{k}^{c}\mid\mathbf{x}_{i})\log\Phi((2y_{i,k}^{c}-1)f_{k}^{c}(\mathbf{x}_{i}))\,df_{k}^{c}$
$\displaystyle=\int\mathcal{N}\left(\mathbf{f}_{k}^{c}\mid\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i}),\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})\right)\log\Phi((2y_{i,k}^{c}-1)f_{k}^{c}(\mathbf{x}_{i}))\,d\mathbf{f}_{k}^{c}$
$\displaystyle=\int\frac{1}{\sqrt{2\pi\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}}\exp\left(-\frac{(\mathbf{f}_{k}^{c}-\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i}))^{2}}{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}\right)\log\Phi((2y_{i,k}^{c}-1)f_{k}^{c}(\mathbf{x}_{i}))\,d\mathbf{f}_{k}^{c}.$

To use Gauss-Hermite quadrature, we need to transform this integral into the
standard form of $\int_{-\infty}^{\infty}e^{-x^{2}}g(x)\,dx$. We can perform a
change of variable
$\mathbf{z}=\frac{\mathbf{f}_{k}^{c}-\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i})}{\sqrt{2\bm{\Sigma}_{\mathbf{f}_{k}^{c}}(\mathbf{x}_{i})}}$
so that
$\mathbf{f}_{k}^{c}=\mathbf{z}\sqrt{2\bm{\Sigma}_{\mathbf{f}_{k}^{c}}(\mathbf{x}_{i})}+\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i})$
and $d\mathbf{f}_{k}^{c}=\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}\,dz$.
The integral becomes:

$B(\mathbf{x}_{i})=\int_{-\infty}^{\infty}\frac{e^{-z^{2}}}{\sqrt{\pi}}\log\Phi\left((2y_{i,k}^{c}-1)\left(z\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}+\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i})\right)\right)\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}\,dz.$

Now, the integral is in a form suitable for Gauss-Hermite quadrature, where

$\sqrt{\bm{\Sigma}_{\mathbf{f}_{k}^{c}}(\mathbf{x}_{i})}\log\Phi((2y_{i,k}^{c}-1)(\mathbf{z}\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}+\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i}))).$

Let $z^{(j)}_{i}$ and $w^{(j)}$ denote the nodes and weights of the Gauss-
Hermite quadrature. The integral $B(\mathbf{x}_{i})$ can be approximated as:

$B(\mathbf{x}_{i})\approx\frac{1}{\sqrt{\pi}}\sum_{l=1}^{L}w^{(l)}\log\Phi((2y_{i,k}^{c}-1)(z^{(l)}_{i}\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})}+\bm{\mu}_{f_{k}^{c}}(\mathbf{x}_{i})))\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}_{i})},$

where $L$ is the number of quadrature points. The $z^{(l)}_{i}$ are the roots
of the physicists’ version of the Hermite polynomial $H_{n}(z_{i})$, and the
associated weights $w^{(l)}$ are given by

$w^{(l)}=\frac{2^{L-1}L!\sqrt{\pi}}{L^{2}[H_{L-1}(z_{i})]^{2}}.$

This approximation allows for the numerical evaluation of the integral using
the Gauss-Hermite quadrature method (Abramowitz et al., 1988). to calculate
the $E_{q_{\phi}(\mathbf{x}_{i})}\left[B(\mathbf{x}_{i})\right]$ we can
draw$J$ sample from $\mathbf{x}_{i}^{(j)}\sim q_{\phi}(\mathbf{x}_{i})$ and
estimate the expectation.

$\displaystyle
E_{q_{\phi}(\mathbf{x}_{i})}\left[B(\mathbf{x}_{i})\right]\approx$
$\displaystyle\frac{1}{\sqrt{\pi}}\sum_{l=1}^{L}\sum_{j=1}^{J}w^{(l)}\log\Phi((2y_{i,k}^{c}-1)(z^{(l)}_{i}\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}^{(j)}_{i})}+$
(24)
$\displaystyle\bm{\mu}_{f_{k}^{c}}(\mathbf{x}^{(j)}_{i})))\sqrt{2\bm{\Sigma}_{f_{k}^{c}}(\mathbf{x}^{(j)}_{i})}.$
(25)

####  2.3.7 Regression Expected Log Likelihood ($ELL^{\text{reg}}$)

In this section, we address the regression component of our model. With a
similar calculation we have the:

$\displaystyle\text{ELL}^{\text{reg}}$
$\displaystyle=E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{r}\mid\mathbf{U}^{r},\mathbf{X})q_{\lambda}(\mathbf{U}^{r})}\left[\log
p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})\right]\right]$ $\displaystyle=\int
q_{\phi}(\mathbf{X})p(\mathbf{F}^{r}\mid\mathbf{U}^{r},\mathbf{X})q_{\gamma}(\mathbf{U}^{r})\sum_{i,d}\log
p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))\,d\mathbf{F}^{r}\,d\mathbf{U}^{r}\,d\mathbf{X}$
$\displaystyle=\sum_{i,d}\int
q_{\phi}(\mathbf{x}_{r})p(\mathbf{f}_{d}^{r}\mid\mathbf{u}_{d}^{r},\mathbf{x}_{i})q_{\gamma}(\mathbf{u}_{d}^{r})\log
p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))\,d\mathbf{f}_{d}^{r}\,d\mathbf{u}_{d}^{r}\,d\mathbf{x}_{i}$
$\displaystyle=\sum_{i,d}ELL_{i,d}^{\text{reg}}.$

Here are the computation details for calculating the ELL:

$\displaystyle ELL_{i,d}^{\text{reg}}$ $\displaystyle=\int
q_{\phi}(\mathbf{x}_{i})p(\mathbf{f}_{d}^{r}\mid\mathbf{u}_{d}^{r},\mathbf{x}_{i})q_{\gamma}(\mathbf{u}_{d}^{r})\log
p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))\,d\mathbf{f}_{d}^{r}\,d\mathbf{u}_{d}^{r}\,d\mathbf{x}_{i}$
$\displaystyle=\int q_{\phi}(\mathbf{x}_{i})\int
q_{\gamma}(\mathbf{u}_{d}^{r})p(\mathbf{f}_{d}^{r}\mid\mathbf{u}_{d}^{r},\mathbf{x}_{i})\,d\mathbf{u}_{d}^{r}\,\log
p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))d\mathbf{f}_{d}^{r}\,d\mathbf{x}_{i}$
$\displaystyle=\int
q_{\phi}(\mathbf{x}_{i})q_{\gamma}(\mathbf{f}_{d}^{r}\mid\mathbf{x}_{i})\log
p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))\,d\mathbf{f}_{d}^{r}\,d\mathbf{x}_{i}$
$\displaystyle=E_{q_{\phi}(\mathbf{x}_{i})}\left[E_{q_{\gamma}(\mathbf{f}_{d}^{r}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,d}^{r}\mid f_{d}^{r}(\mathbf{x}_{i}))\right]\right]$
$\displaystyle=E_{q_{\phi}(\mathbf{x}_{i})}\left[A(\mathbf{x}_{i})\right].$
(26)

We already know that $y_{i,d}^{r}=f_{d}^{r}(\mathbf{x}_{i})+\epsilon_{d}$,
therefore the internal term has a Gaussian distribution

$p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))=\mathcal{N}(f_{d}^{r}(\mathbf{x}_{i}),\bm{\Sigma}_{d}^{2}).$

from Appendix C we know the predictive distribution.

$\displaystyle q_{\gamma}(f_{d}^{r}\mid\mathbf{x}_{i})$
$\displaystyle=\mathcal{N}\left(f_{d}^{r}\mid\bm{\mu}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i}),\bm{\Sigma}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i})\right)$
$\displaystyle\bm{\mu}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i})$
$\displaystyle=\mathbf{k}_{iM_{r}}^{r}(\mathbf{K}_{M_{r}M_{r}}^{r})^{-1}\mathbf{m}^{r}_{d}$
$\displaystyle\bm{\Sigma}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i})$
$\displaystyle=k_{ii}^{r}+\mathbf{k}_{iM_{r}}^{r}(\mathbf{K}_{M_{r}M_{r}}^{r})^{-1}(\mathbf{S}_{d}^{r}-\mathbf{K}_{M_{r}M_{r}}^{r})(\mathbf{K}_{M_{r}M_{r}}^{r})^{-1}\mathbf{k}_{M_{r}i}^{r},$

where $\mathbf{k}^{r}_{iM_{r}}=k^{r}(\mathbf{x}_{i},\mathbf{Z}^{r})$,
$k_{ii}^{r}=k^{r}(\mathbf{x}_{i},\mathbf{x}_{i})$, and $K_{M_{r}M_{r}}^{r}$ is
the covariance matrix evaluated at inducing points.

$\displaystyle A(\mathbf{x}_{i})$
$\displaystyle=E_{q_{\gamma}(f_{d}^{r}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,d}^{r}\mid f_{d}^{r}(\mathbf{x}_{i}))\right]$ $\displaystyle=\int
q_{\gamma}(f_{d}^{r}\mid\mathbf{x}_{i})\log p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{x}_{i}))\,df_{d}^{r}$ $\displaystyle=-0.5\left[\log
2\pi+\log\bm{\Sigma}_{d}^{2}+\frac{(y_{i,d}^{r}-\bm{\mu}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i}))^{2}+\bm{\Sigma}_{\mathbf{f}_{d}^{r}}(\mathbf{x}_{i})}{\bm{\Sigma}_{d}^{2}}.\right].$

We can easily calculate this by sampling.

$\displaystyle ELL_{i,d}^{\text{reg}}$
$\displaystyle=E_{q_{\phi}(\mathbf{x}_{i})}\left[E_{q_{\gamma}(f_{d}^{r}\mid\mathbf{x}_{i})}\left[\log
p(y_{i,d}^{r}\mid f_{d}^{r}(\mathbf{x}_{i}))\right]\right]$
$\displaystyle=-\frac{1}{2J}\sum_{j=1}^{J}\left[\log
2\pi+\log\bm{\Sigma}_{d}^{2}+\frac{(y_{i,d}^{r}-\bm{\mu}_{\mathbf{f}_{d}^{r}}(x^{(j)}_{i}))^{2}+\bm{\Sigma}_{\mathbf{f}_{d}^{r}}(x^{(j)}_{i})}{\bm{\Sigma}_{d}^{2}}\right].$
(27)

It is worth mentioning that for some certain kernel functions, finding the
analytic form for the expectation introduced in Equation (26) is possible for
certain kernel functions. This concludes our analysis and approximation of the
Regression expected log-likelihood for the model.

####  2.3.8 Training procedure

Learning occurs through optimizing the kernel hyperparameters, inducing
points, and variational parameters by maximizing the ELBO as the objective
function.

$\operatorname*{arg\,max}_{\mathbf{\theta},\mathbf{\phi},\mathbf{\lambda},\mathbf{\gamma},\mathbf{Z}^{r},\mathbf{Z}^{c}}\text{ELBO}$

Where $\mathbf{\theta}$ represents the kernel hyperparameters,
$\mathbf{\phi}=\\{(\bm{\mu}_{i},\mathbf{s}_{i}),i=1,\ldots,N\\}$ are the
variational parameters for latent features,
$\mathbf{\lambda}=\\{(\mathbf{m}_{k}^{c},\mathbf{S}_{k}^{c}),k=1,\ldots,K\\}$
and
$\mathbf{\gamma}=\\{(\mathbf{m}_{d}^{r},\mathbf{S}_{d}^{r}),d=1,\ldots,D\\}$
are the variational parameters for inducing variables. In the previous
section, we discussed how the expected log-likelihood for regression
($\text{ELL}^{reg}$) and classification ($\text{ELL}^{cls}$) can be broken
down and calculated for each sample in the dataset. This formulation
simplifies the calculations and enables batch training and sampling methods.
To optimize the variational parameters for $q_{\phi}(\mathbf{X})$ during
sampling, We generate samples $\mathbf{x}_{i}^{(j)}\sim
q_{\phi}(\mathbf{x}_{i})$ using the reparameterization technique , where
$\epsilon^{(j)}\sim\mathcal{N}(0,\mathbb{I}_{Q})$ and
$\mathbf{x}_{i}^{(j)}=\bm{\mu}_{i}+\mathbf{s}_{i}\odot\epsilon^{(j)}$. Based
on Equation (22), (26) , and (23), the training procedure is simplified to
this maximization problem:

$\displaystyle\operatorname*{arg\,max}_{\mathbf{\theta},\mathbf{\phi},\mathbf{\lambda},\mathbf{\gamma},\mathbf{Z}^{r},\mathbf{Z}^{c}}$
$\displaystyle\sum_{i,d}\text{ELL}_{i,d}^{\text{reg}}+\sum_{i,k}\text{ELL}_{i,k}^{\text{cls}}-\sum_{i=1}^{N}\text{KL}(q_{\phi}(\mathbf{x}_{i})\parallel
p(\mathbf{x}_{i}))$
$\displaystyle-\sum_{k=1}^{K}\text{KL}(q_{\lambda}(\mathbf{u}_{k}^{c})\parallel
p(\mathbf{u}_{k}^{c}\mid\mathbf{Z}^{c}))-\sum_{d=1}^{D}\text{KL}(q_{\gamma}(\mathbf{u}_{d}^{r})\parallel
p(\mathbf{u}_{d}^{r}\mid\mathbf{Z}^{r}))$ (28)

where the KL divergence terms have closed forms due to the conjugacy of the
Gaussian distributions.$\text{ELL}_{i,d}^{\text{reg}}$ and
$\text{ELL}_{i,k}^{\text{cls}}$ was introduced in Equation (2.3.7) and (24).
The Algorithm 1 details the training procedure. Detailed implementations of
the predictive distribution, along with the calculations for the expected log-
likelihood for regression and classification, are provided in Appendices C and
D, respectively.

Algorithm 1 Training the Gaussian process Regression Model

1:

2:Batch size ($B$), Learning rate ($lr$), Number of iterations ($N_{iter}$)

3:Training continuous measures $\mathbf{Y}^{r}\in\mathbb{R}^{N\times D}$

4:Training categorical observation or class labels
$\mathbf{Y}^{c}\in\mathbb{R}^{N\times K}$

5:Optimized model parameters
$\\{\theta^{*},\psi^{*},\phi^{*},\lambda^{*},\gamma^{*},\mathbf{Z}^{*}\\}$

6:procedure Train
Model($\mathbf{Y}^{r},\mathbf{Y}^{c},\mathbf{Z}^{c},\mathbf{Z}^{r},lr,N_{iter}$)

7: Initialize of inducing points $\mathbf{Z}^{c}\in\mathbb{R}^{M_{c}\times Q}$
and $\mathbf{Z}^{r}\in\mathbb{R}^{M_{r}\times Q}$

8: Initialize kernel parameters $\mathbf{\psi}$ and $\mathbf{\theta}$

9: Initialize likelihood noise variance $\bm{\Sigma}_{y}$

10: Initialize
$\phi=\\{(\bm{\mu}_{i},\mathbf{s}_{i})\mid\bm{\mu}_{i}\in\mathbb{R}^{Q},\mathbf{s}_{i}\in\mathbb{R}^{Q}\\}_{i=1}^{N}$
for latent variables $\mathbf{X}$

11: Initialize
$\lambda=\\{(\mathbf{m}_{d}^{r},\mathbf{S}_{d}^{r})\mid\mathbf{m}_{d}^{r}\in\mathbb{R}^{M_{r}},\mathbf{S}_{d}^{r}\in\mathbb{R}^{M_{r}\times
M_{r}}\\}_{d=1}^{D}$ for inducing variables $\mathbf{U}^{r}$

12: Initialize
$\gamma=\\{(\mathbf{m}_{k}^{c},\mathbf{S}_{k}^{c})\mid\mathbf{m}_{k}^{c}\in\mathbb{R}^{M_{c}},\mathbf{S}_{k}^{c}\in\mathbb{R}^{M_{c}\times
M_{c}}\\}_{k=1}^{K}$ for inducing variables $\mathbf{U}^{c}$

13: Initialize $ELBO$ to empty list

14: for $i=1$to$N_{iter}$ do

15: Sample $\bm{\epsilon}_{i}\sim\mathcal{N}(\mathbf{0},\mathbb{I}_{Q})$ for
$i=1,\ldots,N$

16:
$\mathbf{x}_{i}\leftarrow\bm{\mu}_{i}+\sqrt{\mathbf{s}_{i}}\bm{\epsilon}_{i}$
for $i=1,\ldots,n$$\triangleright$ Re-parameterization trick

17: $\mathbf{X}\leftarrow\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\\}$

18: Get a random batch $\mathbf{X}_{b}$

19: Calculate $K^{r}_{NN}$, $\mathbf{K}_{M_{r}M_{r}}^{r}$,
$\mathbf{K}_{M_{r}N}^{r}$, $\mathbf{K}^{c}_{NN}$,
$\mathbf{K}^{c}_{M_{c}M_{c}}$, $\mathbf{K}^{c}_{M_{c}N}$

20:
$\bm{\mu}_{\mathbf{f}^{r}},\bm{\Sigma}_{\mathbf{f}^{r}}\leftarrow$PredictiveDistribution($\mathbf{X}_{b}$,
$\mathbf{K}^{r}_{NN}$, $\mathbf{K}^{r}_{M_{r}M_{r}}$,
$\mathbf{K}^{r}_{M_{r}N}$, $\mathbf{m}^{r}$, $\mathbf{S}^{r}$)

21:
$\bm{\mu}_{\mathbf{f}^{c}},\bm{\Sigma}_{\mathbf{f}^{c}}\leftarrow$PredictiveDistribution($\mathbf{X}_{b}$,
$\mathbf{K}^{c}_{NN}$, $\mathbf{K}^{c}_{M_{c}M_{c}}$,
$\mathbf{K}^{c}_{M_{c}N}$, $\mathbf{m}^{c}$, $\mathbf{S}^{c}$)

22:
$\text{ELL}^{\text{reg}}\leftarrow$ExpectedLogLikelihoodRegression($\mathbf{Y}^{r}$,
$\bm{\mu}_{\mathbf{f}^{r}},\bm{\Sigma}_{\mathbf{f}^{r}}$,
$\bm{\Sigma}_{d}^{2}$, $\mathbf{X}$)

23:
$\text{ELL}^{\text{cls}}\leftarrow$ExpectedLogLikelihoodClassification($\mathbf{Y}^{c}$,
$\bm{\mu}_{\mathbf{f}^{c}},\bm{\Sigma}_{\mathbf{f}^{c}},\mathbf{X}$)

24: $\text{KL}_{U}\leftarrow\text{KL}(q_{\lambda}(\mathbf{U}^{c})\parallel
p(\mathbf{U}^{c}\mid\mathbf{Z}^{c}))+\text{KL}(q_{\gamma}(\mathbf{U}^{r})\parallel
p(\mathbf{U}^{r}\mid\mathbf{Z}^{r}))$

25:
$\text{ELBO}\leftarrow\text{ELL}^{\text{reg}}+\text{ELL}^{\text{cls}}-\text{KL}_{U}-\text{KL}(q_{\phi}(\mathbf{X})\parallel
p(\mathbf{X}))$

26: Compute gradient of $-ELBO$ w.r.t
$\\{\theta,\psi,\phi,\lambda,\gamma,\mathbf{Z}\\}$

27: Update $\\{\theta,\psi,\phi,\lambda,\gamma,\mathbf{Z}\\}$ using gradient
descent with learning rate $lr$

28: Append $-ELBO$to$ELBO$list

29: end for

30:
return$\\{\theta^{*},\psi^{*},\phi^{*},\lambda^{*},\gamma^{*},\mathbf{Z}^{*}\\}$

31:end procedure

####  2.3.9 Prediction

After optimizing the hyperparameters using a training dataset, we aim to
employ LDGD to predict the labels for a new set of measurements, denoted as
${\mathbf{Y}^{r}}^{*}$. This set of measurements constitutes the test set. We
assume the absence of corresponding labels ${\mathbf{Y}^{c}}^{*}$ for these
new data points. In essence, we employ the test data points to infer their
latent features and subsequently attempt to decode these latent features to
determine the distribution over the corresponding labels.

The variational posterior for the latent variables $\mathbf{X}$ was introduced
in Equation (18). Given a new test point ${\mathbf{Y}^{r}}^{*}$, we seek to
compute the posterior distribution
$p(\mathbf{X}^{*}\mid{\mathbf{Y}^{r}}^{*},\mathbf{Y}^{r},\mathbf{Y}^{c})$,
which represents our belief about the latent variables after observing the new
test data. This can be approximated as
$q(\mathbf{X}^{*})=\prod_{i=1}^{N_{test}}\mathcal{N}(\mathbf{x}_{i}^{*}\mid\bm{\mu}_{*},\mathbf{s}_{*}\mathbb{I}_{Q})$.

To find the optimum $\bm{\mu}_{*}$ and $\mathbf{s}_{*}$, we aim to minimize
the negative of the lower bound concerning these new variational parameters.
During this process, the inducing points and inducing variables are held
constant. Therefore, the objective function for optimizing the variational
parameters associated with a test point is formulated as follows:

$\arg\max_{\bm{\mu}_{*},\mathbf{s}_{*}}\left\\{\sum_{d}\text{ELL}_{*,d}^{\text{reg}}-\text{KL}\left(q_{\phi}(\mathbf{X}^{*})\parallel
p(\mathbf{X}^{*})\right)\right\\}.$

Using this result, we can effectively decode and determine
$p({\mathbf{Y}^{c}}^{*}\mid\mathbf{X}^{*})$.

####  2.3.10 Fast LDGD

In LDGD, we need to run an iterative optimization in the prediction or test
phase, making it less suitable in a real-time scenario. In this section, we
introduce Fast LDGD, in which we employ a neural network to estimate the
variational parameters for latent space directly. To this end, we modify the
variational distributions introduced in Equation (20) and estimate the
posterior $p(\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})$ with

$q_{\phi}(\mathbf{X})=\prod_{i=1}^{N}\mathcal{N}(\mathbf{x}_{i}\mid\mu(\mathbf{y}_{i}^{r}),s(\mathbf{y}_{i}^{r})\mathbb{I}_{Q}),$
(29)

Where $\mu(\mathbf{y}_{i}^{r})$ and $s(\mathbf{y}_{i}^{r})$ are neural network
outputs with trainable parameters $\mathbf{\phi}$. The architecture of the
neural network can be chosen based on the data. In this approach, the number
of the model parameters depends on the continuous data dimensionality and does
not depend on the number of data points. In contrast, in LDGD, the number of
variational parameters in Equation (20) grows linearly with the number of
training samples. Thus, this approach is beneficial when dealing with a
decoding problem (classification) with a relatively larger number of data
points. Our goal in the training phase is to minimize the negative of the ELBO
by optimizing neural network parameters, inducing points, inducing variables’
variational parameters, kernel parameters, and model parameters. By achieving
this, we ensure the neural network can estimate a suitable approximation to
the true posterior given the continuous data. Once trained, the model can
instantly provide estimates of the posterior’s mean and covariance for a test
point without additional optimization. This capability facilitates real-time
inference.

This approach has several advantages. First, by directly mapping from
observations to latent variables’ variational parameters, we bypass the need
for iterative optimization to determine $\mathbf{\mu}_{i}$ and
$\mathbf{s}_{i}$ (introduced in Equation (20)) for each test data point.
Consequently, this significantly reduces computational time, enabling real-
time predictions, which is crucial in time-sensitive applications.

##  3 Experiments and Results

In this section, we study different attributes of LDGD by applying it to
various datasets. In particular, we examine its capability to draw low-
dimensional data representations, decode accurately (classification), and
generate data points. We investigate these attributes using both synthetic and
real-world datasets. By utilizing synthetic data, we can conduct a
comprehensive investigation of the model attributes under controlled
conditions. The real-world datasets include the Oil Flow (Bishop and James,
1993), Iris (Fisher, 1988), and MNIST datasets (Deng, 2012). Through these
datasets, we further explore how the model’s attributes can extend to handling
complex datasets, thereby underscoring the unique benefits of LDGD. These
benefits include its ability to infer meaningful features from high-
dimensional data, address variability present in the data, and achieve optimal
low-dimensional data representation.

We also compare LDGD modeling results with different variants of GPLVM.
Through these comparisons, we aim to underscore LDGD’s distinctive advantages.
In applying this model to classification tasks, we evaluate it across three
datasets using classification metrics to demonstrate the model’s robustness in
class label prediction. Additionally, we compared LDGD classification accuracy
with two other supervised GPLVM variants, SLLGPLVM (Jiang et al., 2012) and
SGPLVM (Gao et al., 2010), to illustrate its superiority in classification. We
explore the model’s data generation capability using the MNIST dataset,
showing that LDGD can effectively replicate the complex patterns of
handwritten digits. We also compare LDGD and fast LDGD with variational auto-
encoders, where we highlight LDGD and fast LDGD’s key features over
variational auto-encoders, including optimal discovery of low dimension
representation, robustness with limited dataset and inference process. These
analyses will illustrate LDGD’s advantages in analyzing complex and high-
dimensional datasets.

###  3.1 Datasets

####  3.1.1 Synthetic data

For the synthetic data, we utilize a ”moon-like” dataset, which is a two-
dimensional dataset resembling a pair of crescent moons depicted in Figure 2.
We transfer this two-dimensional data into much higher dimensional space. To
do this, we use two different transformations. The first transformation
involves a linear transformation, where the original data is multiplied by a
randomly generated matrix of size $D\times 2$, where $D$ is set to 5, 15, or
20. In the second transformation, we double the dimension of the data by
adding white noise data dimensions. For instance, if the data dimension is 5,
we generate five channels of independent white noise with unit variance and
append it to the data. The added dimensions do not carry class-specific
information, making inference and classification tasks hard. For each
specified dimensionality, we generated 500 synthetic data samples, divided
them into two classes with 250 samples per class, and labeled them
synthetic-10, synthetic-30, and synthetic-40.

####  3.1.2 Oil Flow Dataset

The Oil Flow dataset is a high-dimensional dataset consisting of
12-dimensional feature vectors. These data points are recorded to capture
various factors representing oil flow dynamics within a pipeline (Bishop and
James, 1993). The dataset embodies three distinct phases of flow: horizontally
stratified, nested annular, and homogeneous mixture flow. These three phases
represent the conditions under which oil, gas, and water can coexist and move
through a pipeline. Despite the high dimensionality of the dataset, the
crucial information about the flow conditions can be effectively described by
two main variables: the fraction of water and the fraction of oil in the flow.
This property of the dataset makes it an extensively used benchmark for
evaluating the effectiveness of dimensionality reduction techniques,
showcasing their ability to simplify and interpret complex, high-dimensional
data in fluid dynamics and pipeline management. A total of 1,000 samples are
included in the dataset.

####  3.1.3 Iris Dataset

The Iris dataset is widely used in machine learning, covering morphological
variations across three categories of iris flowers, namely Setosa, Versicolor,
and Virginica (Fisher, 1988). It comprises 150 data points distributed equally
among the species, each detailing four critical features reflecting the
flowers’ physical characteristics, specifically petal and sepal sizes. The
dataset’s clear class distinctions and simple structure make it an ideal
example of how our dimensionality reduction technique finds meaningful low-
dimensional representation.

####  3.1.4 MNIST Dataset

The MNIST (Modified National Institute of Standards and Technology) dataset is
popularly used in machine learning and image processing (Deng, 2012). The
dataset comprises 70,000 images of handwritten digits (0 through 9), and it
serves as a fundamental benchmark for evaluating the performance of various
machine-learning algorithms. Each image is represented in a gray-scale format
with a resolution of 28x28 pixels, equating to 784 dimensions when flattened
into a vector form.

###  3.2 Performance Assessment

####  3.2.1 Synthetic Data

##### Lower Dimension Representation.

We argue that two key advantages of our model are its ability to infer
meaningful information in high-dimensional data and, more importantly,
quantify the validity of its probabilistic inference. To demonstrate this,
LDGD is applied on the synthetic-10, synthetic-30, and synthetic-40 datasets,
described in Section 3.1.2. We run LDGD with two different settings. In the
first setting, we assumed $Q$, the latent dimension, is 2. This choice of
dimension is helpful to properly visualize the inferred $\mathbf{X}$. Figure
2-B, 2-C, and 2-D respectively show inferred $\mathbf{X}$ for the
synthetic-10, synthetic-30, and synthetic-40 dataset. The inferred points
(posterior mean of $\mathbf{X}$) represent a rotated version of the initial
data shown in Figure 2-A. Here, our pipeline converges to a similar data
representation independent of the observed data dimension. The variance of the
predictive distribution $p(\mathbf{Y}^{c}\mid X)$ serves as a measure of
uncertainty in label prediction. The heatmap, employing a binary colormap,
visualizes the uncertainty levels across different regions. The inferred
posterior variance reflects higher confidence in areas with denser data
concentrations (segments) and lower confidence in regions where these
concentrations approach one another.

We use the second setting to highlight LDGD’s capability to infer the latent
dimension representing observed data and corresponding labels optimally. This
can be done by adjusting ARD coefficients in the training step. For this
setting, we use the synthetic-10 dataset and set LDGD’s latent dimensions to
10, equal to the data dimension, with 25 inducing points for each
classification and regression path. The model was trained using 80% of the
dataset, selected randomly. The training process aimed to minimize Equation
(22), as described in Section 2.3.5. Figure 2-E shows values of ELBO loss over
training iterations. For the classification, the ARD coefficients indicate
LDGD converges to 2-D representation aligned with what is embedded in the data
(Figure 2-F). For the regression, the ARD coefficients suggest that all
dimensions are employed for data generation (Figure 2-G). The ARD coefficients
in each of the two paths, along with the inferred state, demonstrate how the
framework tries to reconstruct both paths through the same latent structure.
The inferred $\mathbf{X}$ shows what the model requires for classification and
what is needed for reconstructing data. The regression ARD coefficients
display larger values for the dimensions identified by the classification ARD
coefficients, suggesting that these dimensions play a significant role in data
reconstruction. However, additional dimensions are also necessary. This is
anticipated, as the expanded data in each dimension is a weighted sum of the
original data, necessitating the involvement of other latent dimensions as
well.

Subsequently, with the optimized hyperparameters, we optimized the partial
ELBO introduced in Equation (2.3.9) to find the estimated posterior over
latent variables for the remaining 20% of the data without any knowledge of
the labels. Figure 2-I demonstrates that, even for test data with unknown
labels, the model could discern the intrinsic pattern hidden within the high
dimensions even for test data with unknown labels.

Figure 2: Visualization of Dimensionality Reduction on Synthetic Data. (A)
Displays the initial two dimensions of the moon-like dataset. (B-D)
Illustrates the latent space heatmap representation across different synthetic
data dimensions: 10 dimensions (B), 20 dimensions (C), and 40 dimensions (D).
Red and blue points represent class 1 samples ($y^{c}=0$) and class 2 samples
($y^{c}=1$), respectively. Green crosses indicate classification-inducing
points, while yellow crosses denote regression-inducing points, with five
inducing points used for each. The green data points are more uniformly
distributed over the space as it is constructing the $\mathbf{Y}^{r}$. The
yellow ones are aligned on two sides of the decision boundary as they are more
needed for a correct classification. The heatmap visualizes the model’s
uncertainty level (posterior variance). (E) Shows the training curve (ELBO
loss) for synthetic data with ten dimensions, where the latent space is also
set to 10 dimensions. (F-G) Depicts the ARD coefficients for the
classification kernel (F) and the regression kernel (G). The trained
coefficients highlight that the model selects two dimensions to represent a
10-dimensional space in a lower-dimensional setting for decoding labels and
employs almost all dimensions to reconstruct data in the original space. (H)
displays a scatter plot of the training points in the lower-dimensional space
for the two most dominant dimensions, while (I) shows the scatter plot for
test points where the labels are unknown.

##### Classification Results.

To assess the model’s performance in classification tasks, LDGD was evaluated
using 5-fold cross-validation. In each fold, we optimized the parameters using
80% of the synthetic dataset’s data points, with the remaining 20% serving as
the test set, where labels were assumed to be unknown. Twenty-five inducing
points within a 10-dimensional latent manifold were utilized for LDGD.
Subsequently, the model’s precision, recall, and F-measure were evaluated on
the test set. These results are detailed in Table 1. The model demonstrated
accurate label prediction on the test set, showcasing another critical feature
of LDGD: its capability to decode labels.

Dataset | Precision | Recall | F1 Score
---|---|---|---
Synthetic-10 | 1.0 | 1.0 | 1.0
Synthetic-20 | 1.0 | 1.0 | 1.0
Synthetic-40 | 1.0 | 1.0 | 1.0
Table 1: Classification performance of LDGD for Synthetic data

####  3.2.2 Oil Flow, Iris, and MNIST Dataset

##### Lower Dimension Representation.

To investigate the LDGD application further in inferring low-dimensional
representation of data, we applied it to the Oil Flow, Iris, and MNIST
datasets. We used an ARD kernel with seven dimensions for the Oil Flow and
Iris datasets. The model exhibits low sensitivity to the choice of latent
dimension size; therefore, we have selected 7, which exceeds the
dimensionality of the Iris dataset but remains below that of the Oil Flow
dataset. We also assume we have ten inducing points for the classification
path and the same number for the regression path. We optimized the parameters
using 80% of the data along with their labels and inferred the latent space
for the remaining 20% without label information. The results indicate that the
trained ARD coefficients suggest an optimal use of about one dimension for the
Iris dataset and two dimensions for the Oil Flow dataset (see Figure 3). We
can see there are different levels of certainty in each axis; for instance, in
Figure 3-B, we can see much higher confidence in x-axis, corresponding to
$5^{th}$ element in ARD space, indicating its significant role in the model’s
predictions. On the other hand, the second largest ARD coefficient is less
significant, leading to lower confidence in y-axis. For Iris, the dimension of
$\mathbf{X}$ is larger than the observed data. Though this is not the case the
framework is designed for, we observed that ARD only keeps two coefficients
and pushes other coefficients toward zero. Though this does not reflect any
attributes of the LDGD, it shows the robustness of the pipeline and its
principled representation of the latent data.

Given the MNIST dataset’s larger feature space, we used 20 dimensions for the
latent space and 150 inducing points for each path. The model utilized almost
all the available dimensions for this dataset, highlighting the necessity of
high dimensions to classify and reconstruct such a large dataset. LDGD’s
representation in 2-D may not exhibit as much distinct separation as T-SNE due
to its generative nature and probabilistic embeddings that prioritize
capturing the global data structure over local neighborhood relationships.

In parallel, we applied PCA, t-SNE, GPLVM, Bayesian GPLVM, FGPLVM, SGPLVM, and
SLLGPLVM to the Oil Flow and Iris datasets. The latent variable’s dimensions
were fixed to two across all models for fair comparisons. SGPLVM, SLLGPLVM,
and LDGD utilize labels during the training process, while others are
unsupervised and agnostic to the labels. Figure 4 illustrates the latent
representations generated by LDGD compared to those obtained using other
models. The Iris dataset, being less complex, allowed almost all models to
find a lower-dimensional representation that demonstrates a clear distinction
between species. So, LDGD’s performance is similar to other models. For the
Oil Flow dataset, LDGD and SLLGPLVM effectively separate the three flow
phases, while some other methods exhibit partial overlap between different
categories. Unlike SLLGPLVM, LDGD incorporates inducing points within its
pipeline, which allows LDGD to scale better than SLLGPLVM. Consequently, a
direct comparison for the MNIST dataset is not feasible since SLLPLVM is ill-
equipped to handle larger datasets.

Figure 3: Comparative analysis of latent space representation in Iris and Oil
Flow dataset. (A) reveals the most dominant latent dimension for the Iris
dataset, identified through ARD coefficients. (B) illustrates the scatter plot
of the two dominant dimensions in the latent space for the Iris dataset,
highlighting the data’s intrinsic clustering. The vertical and horizontal
error bars show the variance at each data point in latent space. (C) displays
the most dominant latent dimension for the Oil dataset, as determined by ARD
coefficients. (D) presents the scatter plot of the two dominant dimensions in
the latent space for the Oil Flow dataset, showcasing its unique distribution.
Figure 4: Two-dimensional visualization of dataset projections using various
dimensionality reduction techniques. The top row displays projections of the
Oil Flow dataset, while the bottom row shows the Iris dataset. Techniques used
include (A) PCA, (B) t-SNE, (C) GPLVM, (D) Bayesian GPLVM, (E) FGPLVM, (F)
SGPLVM, (G) SLLGPLVM, and (H) LDGD. We observed SLLGPLVM and LDGD inference in
low dimension reflecting the class label, ending to separate regions for
different data classes. With LDGD, we can even have a one-dimensional
representation of data. Thus, we might even get a better reduction rate using
LDGD.

##### Classification Results.

Here, we aim to demonstrate the model’s ability to decode or classify real-
world datasets. We have applied it to the Iris and Oil Flow with seven latent
dimensions and ten inducing points. Then, we applied LDGD with 20 latent
variables and 150 inducing points for classification and regression for the
MNIST dataset. The evaluation criteria included classification accuracy,
precision, recall, and F1 score. The results are provided in Table 2. As we
see, the model performance is as good as the state of the arts in the iris and
oil dataset, and the performance in the MNIST dataset is not far from the
other image classifier methods like convolutional neural network (CNN) methods
but comparable to them Alvear-Sandoval et al. (2019). The classification
results in datasets with different input dimensions and different data points
further demonstrate the model potential in the feature extraction and
classification step.

Dataset | Accuracy | Precision | Recall | F1 Score
---|---|---|---|---
Iris | 1.00 | 1.00 | 1.00 | 1.00
Oil Flow | 0.99 | 1.00 | 0.99 | 0.99
MNIST | 0.95 | 0.94 | 0.93 | 0.94
Table 2: Classification performance of LDGD for Iris, Oil Flow, and MNIST
dataset

To quantitatively compare the classification capabilities of our model, we
benchmarked its performance against several state-of-the-art models. While our
model and others perform well on the simpler Iris dataset, LDGD stands out for
its unique ability to handle the complexity of the MNIST dataset.
Consequently, our comparison focuses on the Oil Flow dataset to provide a fair
comparison of LDGD’s performance with its competitors.

Table 3 summarizes the classification results, where LDGD consistently
outperforms other models across all measures. This confirms LDGD’s
effectiveness in extracting latent features and accurately predicting labels.
The model’s superior classification performance results from its adaptive
feature extraction capabilities. These empirical findings highlight LDGD’s
dual role as both an efficient tool for dimensionality reduction and a
powerful discriminative model. The model’s achievement of high classification
metrics across various datasets underscores its adaptability and potential for
use in numerous practical machine-learning scenarios, even with larger
datasets. Our results for MNIST, a comparatively large dataset for Gaussian
process-based models, demonstrate LDGD’s scalability. The LDGD model surpasses
both SLLGPLVM and SGPLVM in terms of scalability. The scalability is
attributed mainly to integrating batch training, which enables efficient
handling of large datasets without compromising the model’s performance or
accuracy.

Model | Accuracy | Precision | Recall | F1 Score
---|---|---|---|---
LDGD | 0.99 | 1.00 | 0.99 | 0.99
SGPLVM | 0.93 | 0.92 | 0.91 | 0.92
SLLGPLVM | 0.95 | 0.94 | 0.93 | 0.94
Table 3: Classification performance of LDGD for the Oil Flow dataset compared
to other models

##### LDGD Data Generation.

So far, our focus has been on inferring the latent space and then utilizing
LDGD for label prediction. It is noteworthy that LDGD is capable of generating
continuous samples given a label. Figure 5 displays three sample data points,
each selected randomly near a class in the latent space within the Iris and
Oil Flow datasets. As observed, the reconstructed samples in high-dimensional
space lie close to the actual data points. This demonstrates how LDGD training
effectively balances the label and continuous paths, ensuring that the
generated data accurately represents its respective classes. Additionally,
Figure 5-C illustrates nine sample digits generated randomly from the MNIST
dataset, highlighting the model’s ability to produce samples from more complex
data.

Figure 5: LDGD Data Generation Analysis. (A-D) These figures present scatter
plots of data points in the 2D latent space for the Iris (A) and the Oil Flow
dataset (C). For each dataset, a sample point near each cluster is randomly
selected (marked with a cross) to illustrate the model’s generative
capabilities. The corresponding high-dimensional reconstructions of these
starred points are shown for the Iris dataset (B) and the Oil Flow dataset
(D), alongside the real values of other points. (E) displays nine randomly
chosen test points in the latent space for the MNIST dataset and their
reconstructed images through the model’s generative path. The generated images
properly reconstruct corresponding digits.

###  3.3 LDGD, Fast LDGD, and Variational Auto-Encoders

While LDGD might resemble a variational auto-encoder (VAE), fundamental
differences between our proposed framework and VAE models must be highlighted
further. In VAEs, the latent manifold has Gaussian priors, similar to LDGD and
fast LDGD, where the posterior over latent variables are approximated using a
Gaussian distribution. However, the distinctions lie in the relationship
between the latent manifold, labels, and high-dimensional data. LDGD employs a
probabilistic framework, whereas VAEs utilize a deterministic framework
(neural networks). Fast LDGD occupies a middle ground, inferring $\mathbf{X}$
through a neural network while leveraging a Gaussian Process (GP) model for
label prediction. Another key difference resides in the loss functions used:
VAEs incorporate mean squared error (MSE), categorical entropy, and Kullback-
Leibler (KL) divergence, whereas LDGD relies on the ELBO. Figure 6.A shows the
probabilistic graphs for these three methods. As we can see in LDGD, we are
trying to find the posterior over $\mathbf{X}$ given $\mathbf{Y}^{r}$ and
$\mathbf{Y}^{c}$.

LDGD has multiple advantages over the VAE model, making it a more suitable
model for analyzing high-dimensional data, specifically when the size of
dataset is limited. Using LDGD, we can identify the optimal dimension of the
latent process as a part of learning. In contrast, we would have to tune the
number of dimensions as a hyperparameter using the VAE model. Inferencing and
training in LDGD is not a data-greedy process; we can build the model with
much smaller data than needed for VAE model training. Compared to a point
estimate in VAE, Bayesian inferencing helps us have higher-order statistics,
such as confidence in our prediction of a label or data point, giving us a
better sense of our predictions.

Another key advantage is the inherent robustness of LDGD against overfitting,
a common pitfall in many machine learning models. LDGD’s use of Bayesian
inferencing inherently guards against overfitting, making LDGD particularly
suited for dealing with small datasets. This framework doesn’t merely fit the
model to the data; it enables the model to correctly infer and adapt to the
data, even when the dataset is small. To assess the model’s performance in
such scenarios, we have used the synthetic-20 dataset. We used a different
amount of training data with the same amount of test data and compared VAE,
LDGD, and fast LDGD. As shown in Figure 6-B and 6-C, in small datasets, LDGD
and fast LDGD’s performance in classification is significantly better. For
large datasets, the performance is comparable. These advantages collectively
position LDGD as a mere alternative and a significant advancement over
traditional encoder-decoder models, especially in high-dimensional data
analysis.

Figure 6: Comparison of LDGD, Fast LDGD, and VAE Models. (A) shows the
graphical models for VAE, fast LDGD, and LDGD from left to right. (B, C)
illustrate the relationship between the number of training samples (x-axis)
and the performance metrics on the test set, specifically accuracy (B) and
F-measure (C), where the blue line represents VAE performance, red indicates
LDGD, and green denotes fast LDGD performance.

##  4 Discussion

In this research, we developed different modeling steps of LDGD, including its
training, inferencing, and decoding (classification). In developing the
framework, we introduced new approaches, including a shared latent variable
with doubly stochastic variational inference with trainable inducing points.
To assess different attributes of the model, we applied the framework in the
Synthetic, Iris, and Oil Flow datasets.

The modeling results in both synthetic and real-world data suggest LDGD’s
capability to correctly infer latent space and identify the latent space’s
proper dimension. The classification results using LDGD are on par with the
state-of-the-art models, and in some of the data points, it even surpasses
other approaches.

We pointed to different attributes and potential advantages of LDGD. We argued
its capabilities in identifying the optimal dimension of the latent process.
For instance, in synthetic data, we showed that the framework would suggest
2-D dimensional latent space, which matches the intrinsic dimensions of the
original data. This can be done using an ARD kernel where its coefficients are
trained during the model training. This is a remarkable feature of the model,
which is not present in other algorithms such as VAEs. We also pointed to the
LDGD’s expressiveness and ability to classify data labels precisely. We showed
that LDGD can be a perfect Iris and Oil Flow dataset classifier, surpassing
state-of-the-art approaches. We can see LDGD as a pipeline with adaptive
feature extraction and a classifier with features capable of reconstructing
the data. Our results demonstrated its capability to process larger datasets,
such as MNIST, efficiently, which was achieved by adopting batch training, the
inducing point concept, and a doubly variational strategy. Specifically, this
approach highlighted LDGD’s scalability to larger datasets, a feature
distinctly lacking in models like SGPLVM and SLLGPLVM. We also explored the
generative prowess of LDGD, that is, its ability to create data from high-
dimensional spaces. Our experiments confirmed that LDGD can adeptly generate
samples from the Iris, Oil Flow, and MNIST datasets. The generated samples
from MNIST were remarkably coherent.

We observed a couple of challenges with LDGD as well. A core challenge is
hyperparameter selection, such as determining the number of inducing points.
As with many in machine learning, the model’s performance hinges on
hyperparameter choices. In Gaussian Process-based models, the number of
inducing points is critical; too few may lead to underfitting and inadequate
data complexity representation, while too many can escalate computational
demands and risk overfitting. Finding the optimal number often requires a
dataset-specific balance, posing practical challenges. Although LDGD improves
scalability, it may struggle with extremely large or high-dimensional datasets
due to the inherent complexity of Bayesian models and Gaussian Processes,
limiting its feasibility in some real-time applications.

LDGD’s reliance on label quality for dimensionality reduction means its
performance is contingent on label accuracy. The model’s efficiency may
diminish in cases of scarce, inaccurate, or biased labels. Despite its
proficiency in uncertainty quantification, the interpretability of complex
models like LDGD can be daunting, especially for users without a machine
learning background. The contribution of various modeling steps to create the
outcome might make the pipeline hard to digest. Furthermore, while LDGD shows
promise for particular datasets, its performance across diverse data types,
such as time series, images, or text, has not been thoroughly tested, leaving
its applicability in these areas yet to be fully established.

Several improvements are suggested to overcome these limitations and extend
the applicability of LDGD. These include the development of methods for the
automatic selection of inducing points, such as Bayesian optimization and
cross-validation approaches, as well as extensions of the model to handle time
series data. Several modifications are suggested to enhance the LDGD model for
time series analysis. Firstly, incorporate temporal dynamics by integrating
autoregressive components, state space models, or using time series-specific
kernels in Gaussian Processes (Wong and Li, 2000; Fard et al., 2023; Gu and
Dao, 2023). Secondly, implementing methods to handle sequential data, possibly
using recurrent neural network structures like LSTMs (Yu et al., 2019) or GRUs
(Dey and Salem, 2017) within the LDGD framework. Thirdly, explore time-variant
inducing points that adapt to changes in time series data, improving the
model’s ability to capture temporal dynamics. Parallel computing, GPU
acceleration, and sparse variational frameworks are also recommended to
enhance algorithmic efficiency and scalability.

In general, it is essential to establish standard benchmarks and evaluation
methodologies to advance the field of dimensionality reduction and ensure the
efficacy of techniques like the LDGD model. Additionally, it’s crucial to
continuously benchmark the enhanced model against the latest state-of-the-art
models in dimensionality reduction. This ongoing comparison will ensure the
model’s competitiveness and help identify potential areas for further
improvement.

Our future work will focus on extending LDGD for specific applications, such
as datasets we record during neuroscience experiments, which are generally
small in size and show complex variability (Ziaei et al., 2023). We also aim
to enhance algorithmic efficiency and conduct robustness and generalizability
studies. In summary, the LDGD model offers promising advancements in
dimensionality reduction and classification, with wide-ranging applications
and potential for significant contributions to various fields. However,
realizing its full potential requires continuous research and development to
address its current limitations and expand its applicability in addressing
complex real-world problems.

##  5 Conclusion

In this research, we showcased the utilities of non-parametric models for
dimensionality reduction, feature discovery, and decoding. Despite the
emergence of powerful machine-learning tools, we still lack principled tools
capable of dealing with intrinsic stochasticity and variability in data or
scenarios where the size of available data is limited. The LDGD framework
addresses these needs and achieves comparable or even higher prediction
accuracy and generative capability. With LDGD, we can study data from domains
such as neuroscience, which are generally complex and limited in size. The
ideas in LDGD can be extended to more realistic data, such as time series,
thus requiring further model development, which is the core of our future
research.

##  6 Code Availability

The source code for our model and usage examples is freely available to the
research community for further exploration and development at our GitHub
repository: https://github.com/Navid-Ziaei/LDGD. This initiative promotes
transparency, reproducibility, and collaborative advancements within the
field.

Acknowledgments

This project is sponsored by the Defense Advanced Research Projects Agency
(DARPA) under cooperative agreement No. N010492402. The content of the
information does not necessarily reflect the position or the policy of the
Government, and no official endorsement should be inferred.

##  7 Appendixes

##  A Exploring the Connection between PPCA and GPLVM

The GPLVM is developed based on the probabilistic principal component analysis
(PPCA) framework. In PPCA, the observed data vector $\mathbf{y}_{i}$ for each
observation $i$ is modeled as a linear transformation of the corresponding
latent variable $\mathbf{x}_{i}$ with additive Gaussian noise
$\bm{\epsilon}_{i}$:

$\mathbf{y}_{i}=\mathbf{W}\mathbf{x}_{i}+\bm{\epsilon}_{i},$

where $\mathbf{y}_{i}\in\mathbb{R}^{D}$ denotes the observations,
$\mathbf{x}_{i}\in\mathbb{R}^{Q}$ represents the low-dimensional
representation, and $\mathbf{W}\in\mathbb{R}^{D\times Q}$ is the
transformation matrix mapping from the low-dimensional space to the high-
dimensional space. The prior distribution over $\mathbf{X}$ is assumed to be
Gaussian with a zero mean and an identity covariance matrix:

$p(\mathbf{X})=\prod_{i=1}^{N}\mathcal{N}(\mathbf{x}_{i}\mid\mathbf{0},\mathbb{I}_{Q}).$

Consequently, the likelihood of $\mathbf{Y}$ given the transformation matrix
$\mathbf{W}$ is:

$p(\mathbf{Y}\mid\mathbf{W})=\prod_{i=1}^{N}\mathcal{N}(\mathbf{y}_{i}\mid\mathbf{0},\mathbf{C}),$

with the covariance matrix $\mathbf{C}$ defined as:

$\mathbf{C}=\mathbf{W}\mathbf{W}^{T}+\sigma^{2}\mathbb{I}_{D}.$

In the dual formulation of PPCA, we consider the model from the perspective of
the observed data:

$\mathbf{y}_{:,d}=\mathbf{X}\mathbf{w}_{:,d}+\bm{\epsilon}_{d},$

where $\mathbf{X}\in\mathbb{R}^{N\times Q}$ represents the low-dimensional
representation of all samples, $\mathbf{y}_{:,d}$ is the $d$th feature in
high-dimensional space for all samples, and $\mathbf{w}_{:,d}$ is the $d$th
column of $\mathbf{W}$, which is the linear transformation matrix. The prior
distribution over $\mathbf{W}$ is assumed to be Gaussian with a zero mean and
an identity covariance matrix:

$p(\mathbf{W})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{w}_{:,d}\mid\mathbf{0},\mathbb{I}_{Q}).$

Here, the likelihood of $\mathbf{Y}$ given the latent variables $\mathbf{X}$
is:

$p(\mathbf{Y}\mid\mathbf{X})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{y}_{:,d}\mid\mathbf{0},\mathbf{K}),$

where $\mathbf{K}$ is the covariance matrix reflecting a linear relationship
to $\mathbf{X}$:

$\mathbf{K}=\mathbf{X}\mathbf{X}^{T}+\sigma^{2}\mathbb{I}_{N}.$

To accommodate the nonlinear relationships present in many real-world
datasets, we can extend the dual PPCA formulation using kernel functions
$k_{\theta}(\mathbf{x},\mathbf{x}^{\prime})$. This formulation enables us to
express the covariance matrices through kernel functions as follows:

$\mathbf{K}=\mathbf{K}_{NN}+\sigma^{2}\mathbb{I}_{N},$

where $\mathbf{K}_{NN}$ denotes the covariance matrix, with each element
$k_{ij}$ in the $i$th row and $j$th column calculated by the kernel function
$k_{\theta}(\mathbf{x}_{i},\mathbf{x}_{j})$, reflecting the covariance between
two points $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ in the input space. The
incorporation of kernel functions into the PPCA model gives rise to the
Gaussian Process Latent Variable Model (GPLVM), allowing for the capture of
nonlinear relationships within the data. This enhancement transforms the model
into a more adaptable framework, capable of accommodating the data’s inherent
complexity, as described by:

$p(\mathbf{Y}\mid\mathbf{X})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{y}_{:,d}\mid\mathbf{0},\mathbf{K}_{NN}+\sigma^{2}\mathbb{I}_{N}),$

where $\mathbf{Y}$ represents the high-dimensional data, $\mathbf{X}$ denotes
the latent low-dimensional representation, and $\mathbf{K}_{NN}$ is the
covariance matrix derived from the kernel function, capturing the
nonlinearities in the data. Thus, GPLVM extends the linear dimensionality
reduction approach of PPCA to a nonlinear setting, offering a more powerful
tool for uncovering the underlying structure in complex datasets.

##  B Evidence Lower Bound for GPLVM

In this appendix, we demonstrate that maximizing the ELBO (Equation 22) is
equivalent to minimizing the KL divergence between the true posterior
distribution and the variational posterior distribution. Equations (17) and
(21) depict LDGD’s true and approximated posteriors, respectively. Let’s
denote the true posterior as $\hat{P}\triangleq
p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X}\mid\mathbf{Y}^{r},\mathbf{Y}^{c})$
and approximated posterior as $\hat{Q}\triangleq
q(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})=q(\mathbf{F}^{c},\mathbf{U}^{c})q(\mathbf{F}^{r},\mathbf{U}^{r})q_{\phi}$.
We begin by computing the KL divergence between the exact posterior and our
approximation:

KL
$\displaystyle(\hat{Q}\parallel\hat{P})=\int{\hat{Q}\log\frac{\hat{Q}}{\hat{P}}d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}}$
(30)

The posterior distribution of LDGD, given the observed data $\mathbf{Y}^{r}$
and $\mathbf{Y}^{c}$, the latent functions $\mathbf{F}^{r}$ and
$\mathbf{F}^{c}$, the inducing variables $\mathbf{U}^{r}$ and
$\mathbf{U}^{c}$, and the latent features $\mathbf{X}$, is defined as:

$\displaystyle\hat{P}=$
$\displaystyle\frac{p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})}{p(\mathbf{Y}^{r},\mathbf{Y}^{c})}.$

We then substitute this into Equation (30):

$\displaystyle\textit{KL}(\hat{Q}\parallel\hat{P})=$
$\displaystyle\int{\hat{Q}\log\frac{p(\mathbf{Y}^{r},\mathbf{Y}^{c})\hat{Q}}{p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})}d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}}$
$\displaystyle=$ $\displaystyle\log
p(\mathbf{Y}^{r},\mathbf{Y}^{c})+\int{\hat{Q}\log\frac{\hat{Q}}{p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})}d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}}$
$\displaystyle\triangleq$ $\displaystyle\log
p(\mathbf{Y}^{r},\mathbf{Y}^{c})+\text{A}$ (31)

We focus on the integral term. We utilize Equations (16) and (21) to factorize
the joint distribution and the approximated posterior, respectively:

$\displaystyle\text{A}=$
$\displaystyle\int{\hat{Q}\log\frac{\hat{Q}}{p(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X},\mathbf{Y}^{r},\mathbf{Y}^{c})}d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}}$
$\displaystyle=$
$\displaystyle\int{\hat{Q}\log\frac{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}}$
$\displaystyle=$ $\displaystyle
E_{\hat{Q}}\left[\log\frac{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}\right]$

Given that the KL divergence is always non-negative, we have:

$\displaystyle\log p(\mathbf{Y}^{r},\mathbf{Y}^{c})\geq
E_{\hat{Q}}\left[\log\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}\right]\triangleq\textit{ELBO}$

Thus, we have identified a lower bound for the marginal likelihood, known as
the Evidence Lower Bound (ELBO). An alternative approach to achieve the same
result involves marginalizing the joint distribution, followed by multiplying
and dividing the joint distribution by the joint posterior approximation:

$\displaystyle p(\mathbf{Y}^{r},\mathbf{Y}^{c})=$ $\displaystyle\int
p(\mathbf{Y}^{r},\mathbf{Y}^{c},\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})\,d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}$
$\displaystyle=$
$\displaystyle\int\frac{\hat{Q}}{\hat{Q}}p(\mathbf{Y}^{r},\mathbf{Y}^{c},\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})\,d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}$
$\displaystyle\int\hat{Q}\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}\,d\mathbf{F}^{r}\,d\mathbf{F}^{c}\,d\mathbf{U}^{r}\,d\mathbf{U}^{c}\,d\mathbf{X}$
$\displaystyle=$ $\displaystyle
E_{\hat{Q}}\left[\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}\right]$

Using Jensen inequality we can find the evidence lower-bound (ELBO):

$\displaystyle\log p(\mathbf{Y}^{r},\mathbf{Y}^{c})\geq$ $\displaystyle
E_{q(\mathbf{F}^{r},\mathbf{F}^{c},\mathbf{U}^{r},\mathbf{U}^{c},\mathbf{X})}\left[\log\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})p(\mathbf{X})}{q_{\phi}(\mathbf{X})q_{\lambda}(\mathbf{U}^{c})q_{\gamma}(\mathbf{U}^{r})}\right]\triangleq\textit{ELBO}$

Having calculated the ELBO through two approaches, we now proceed to decompose
the ELBO into simpler terms:

$\displaystyle\textit{ELBO}=$ $\displaystyle
E_{q(\mathbf{F}^{c},\mathbf{U}^{c},\mathbf{X})}\left[\log\frac{p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})p(\mathbf{U}^{c})}{q_{\lambda}(\mathbf{U}^{c})}\right]+E_{q(\mathbf{F}^{r},\mathbf{U}^{r},\mathbf{X})}\left[\log\frac{p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})p(\mathbf{U}^{r})}{q_{\lambda}(\mathbf{U}^{c})}\right]$
$\displaystyle-\textit{KL}(q_{\phi}(\mathbf{X})\parallel p(\mathbf{X}))$
$\displaystyle=$ $\displaystyle
E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{c}\mid\mathbf{U}^{c},\mathbf{X})q_{\lambda}(\mathbf{U}^{c})}\left[\log
p(\mathbf{Y}^{c}\mid\mathbf{F}^{c})\right]\right]+$ $\displaystyle
E_{q_{\phi}(\mathbf{X})}\left[E_{p(\mathbf{F}^{r}\mid\mathbf{U}^{r},\mathbf{X})q_{\gamma}(\mathbf{U}^{r})}\left[\log
p(\mathbf{Y}^{r}\mid\mathbf{F}^{r})\right]\right]-$
$\displaystyle\textit{KL}(q_{\lambda}(\mathbf{U}^{c})\parallel
p(\mathbf{U}^{c}))-\textit{KL}(q_{\lambda}(\mathbf{U}^{r})\parallel
p(\mathbf{U}^{r}))-\textit{KL}(q_{\phi}(\mathbf{X})\parallel p(\mathbf{X}))$
$\displaystyle=$
$\displaystyle\text{ELL}^{\text{reg}}+\text{ELL}^{\text{cls}}-\text{KL}_{u}^{c}-\text{KL}_{u}^{r}-\text{KL}_{X}.$

##  C Computation of Predictive Distribution

In this appendix, we detail the calculation of the predictive distribution.
Given that the calculation for both continuous and labeled data follows the
same process, we focus on the continuous case for simplicity. The LDGD
framework employs a Gaussian prior for the inducing variables, formalized as
$p(\mathbf{u}_{:,d}^{r}\mid\mathbf{Z})=\mathcal{N}(\mathbf{u}_{:,d}^{r}\mid\mathbf{0},\mathbf{K}_{M_{r}M_{r}}^{r})$.
Moreover, the variational distribution within this model is specified as
$q_{\gamma}(\mathbf{U}^{r})=\prod_{d=1}^{D}\mathcal{N}(\mathbf{u}_{:,d}^{r}\mid\mathbf{m}_{d}^{r},\mathbf{S}_{d}^{r})$.
Utilizing the transformation matrix $A$, defined as
$A\triangleq{\mathbf{K}_{NM_{r}}^{r}}{\mathbf{K}_{M_{r}M_{r}}^{r}}^{-1}$, the
distribution of function values, conditional on the inducing points and latent
variables, is then expressed as:

$p(\mathbf{f}_{:,d}^{r}\mid\mathbf{u}_{:,d}^{r},\mathbf{X})=\mathcal{N}(\mathbf{A}\mathbf{u}_{:,d}^{r},\mathbf{K}_{NN}^{r}-\mathbf{A}\mathbf{K}_{M_{r}M_{r}}^{r}\mathbf{A}^{T}).$

This formulation leverages Gaussian distribution properties for the marginal
distribution. Subsequently, the predictive distribution, represented as
$q_{\gamma}(\mathbf{f}_{:,d}^{r}\mid\mathbf{X})$, is derived by integrating
out the inducing variables. Given that both terms within the integral are
Gaussian distributions, this integration has an analytical solution, allowing
for straightforward computation:

$\displaystyle q_{\gamma}(\mathbf{f}_{:,d}^{r}\mid\mathbf{X})$
$\displaystyle\triangleq\int
p(\mathbf{f}_{:,d}^{r}\mid\mathbf{u}_{:,d}^{r},\mathbf{X})q_{\gamma}(\mathbf{u}_{:,d}^{r})\,d\mathbf{u}_{:,d}^{r}=\mathcal{N}(\mathbf{f}_{:,d}^{r}\mid\bm{\mu}_{f}^{r},\bm{\Sigma}_{f}^{r}),$
$\displaystyle\bm{\mu}_{f}^{r}$
$\displaystyle=\mathbf{A}^{T}\mathbf{m}_{d}^{r},$
$\displaystyle\bm{\Sigma}_{f}^{r}$
$\displaystyle=\mathbf{K}_{NN}^{r}+\mathbf{A}^{T}(\mathbf{S}^{r}_{d}-{\mathbf{K}_{M_{r}M_{r}}^{r}})\mathbf{A}.$

For stabilizing the learning of variational parameters, we utilize the
whitening trick proposed by Murray and Adams (2010). Let $\mathbf{L}^{r}$
denote the Cholesky factor of the covariance matrix of the prior over inducing
variables, such that
${\mathbf{K}_{M_{r}M_{r}}^{r}}=\mathbf{L}^{r}(\mathbf{L}^{r})^{T}$. We
introduce whitened parameters
$\mathbf{m}^{r}_{d}=\mathbf{L}^{r}\hat{\mathbf{m}}^{r}_{d}$ and
$\mathbf{W}^{r}_{d}=\mathbf{L}^{r}\hat{\mathbf{W}}^{r}_{d}$, where
$\hat{\mathbf{m}}^{r}_{d}$ and $\hat{\mathbf{W}}^{r}_{d}$ are our free
parameters. Here, $\mathbf{W}^{r}_{d}$ and $\hat{\mathbf{W}}^{r}_{d}$ are
upper triangular matrices. The covariance matrix of the variational
distribution is computed as
$\mathbf{S}^{r}_{d}=\mathbf{W}^{r}_{d}(\mathbf{W}^{r}_{d})^{T}$. Given these
definitions and transformations, the mean and covariance of the predictive
distribution are expressed as
$\mathbf{\mu}^{r}_{f}=\hat{\mathbf{A}}^{T}\hat{\mathbf{m}}^{r}_{d}$ and
$\mathbf{\Sigma}^{r}_{f}=\mathbf{K}_{NN}^{r}+\hat{\mathbf{A}}^{T}\left(\hat{\mathbf{W}}^{r}_{d}(\hat{\mathbf{W}}^{r}_{d})^{T}-\mathbb{I}\right)\hat{\mathbf{A}}$
, where $\hat{\mathbf{A}}={L^{r}}^{-1}\mathbf{K}^{r}_{M_{r}N}$ . The
introduction of whitened parameterization and the transformation of the mean
vector simplify computation, ensuring numerical stability and computational
efficiency. This approach is applicable to categorical data (labels) as well:

$\displaystyle q_{\lambda}(\mathbf{f}_{:,k}^{c}\mid\mathbf{X})$
$\displaystyle\triangleq\int
p(\mathbf{f}_{:,k}^{c}\mid\mathbf{u}_{:,k}^{c},\mathbf{X})q_{\lambda}(\mathbf{u}_{:,k}^{c})\,d\mathbf{u}_{:,k}^{c}=\mathcal{N}(\mathbf{f}_{:,k}^{c}\mid\bm{\mu}_{f},\bm{\Sigma}_{f}),$
$\displaystyle\bm{\mu}^{c}_{f}$ $\displaystyle=B^{T}\mathbf{m}_{k}^{c},$
$\displaystyle\bm{\Sigma}^{c}_{f}$
$\displaystyle=\mathbf{K}_{NN}^{c}+B^{T}(\mathbf{S}^{c}_{k}-{\mathbf{K}_{M_{c}M_{c}}^{c}}),B$

where $B\triangleq{\mathbf{K}_{M_{c}M_{c}}^{c}}^{-1}\mathbf{K}_{M_{c}N}^{c}$.
Using whitened parameters, we can express it as
$\bm{\mu}_{f}^{c}=\hat{\mathbf{B}}^{T}\hat{\mathbf{m}}_{k}^{c}$ and
$\bm{\Sigma}_{f}^{c}=\mathbf{K}_{NN}^{c}+\hat{\mathbf{B}}^{T}(\hat{\mathbf{W}}_{k}^{c}(\hat{\mathbf{W}}_{k}^{c})^{T}-\mathbb{I})\hat{\mathbf{B}}$,
where $\hat{\mathbf{B}}=(\mathbf{L}^{c})^{-1}\mathbf{K}_{M_{c}N}^{c}$ and
$\mathbf{K}_{M_{c}M_{c}}^{c}=\mathbf{L}^{c}(\mathbf{L}^{c})^{T}$.

Implementation details for the predictive distribution are provided in
Algorithm 2. The calculation process is identical for both continuous and
labeled data; hence, we omit the superscripts ’r’ and ’c’ used in the main
text to differentiate between the two data types.

Algorithm 2 Calculate Predictive Distribution

1:

2:Sampled latent variables $\mathbf{X}_{b}$

3:Kernel matrices $\mathbf{K}_{NN},{\mathbf{K}_{MM}},\mathbf{K}_{MN}$

4:Whitened variational mean $\hat{\mathbf{m}}$

5:Whitened variational covariance
$\hat{\mathbf{S}}=\hat{\mathbf{W}}\hat{\mathbf{W}}^{T}$ ($\hat{\mathbf{W}}$ is
upper triangular)

6:function PredictiveDistribution($\mathbf{X}_{b}$, $\mathbf{K}_{NN}$,
${\mathbf{K}_{MM}}$, $\mathbf{K}_{MN}$, $\hat{\mathbf{m}}$,
$\hat{\mathbf{S}}$)

7: Compute kernel matrices $\mathbf{K}_{NN},{\mathbf{K}_{MM}},\mathbf{K}_{MN}$
using kernel function

8: Add jitters to ${\mathbf{K}_{MM}}$

9: Decompose ${\mathbf{K}_{MM}}$ into $LL^{T}$ using Cholesky decomposition

10: Compute interpolation term $\hat{\mathbf{A}}\leftarrow
L^{-1}\mathbf{K}_{MN}$

11: Calculate predictive mean
$\mathbf{\mu}_{f}\leftarrow\hat{\mathbf{A}}^{T}\hat{\mathbf{m}}$

12: Calculate predictive covariance
$\mathbf{\Sigma}_{f}\leftarrow\hat{\mathbf{A}}^{T}(\hat{\mathbf{W}}\hat{\mathbf{W}}^{T}-I)\hat{\mathbf{A}}$

13: $\mathbf{\Sigma}_{f}\leftarrow\mathbf{K}_{NN}^{r}+\mathbf{\Sigma}_{f}$

14: return $\mathbf{\mu}_{f}$ and $\mathbf{\Sigma}_{f}$

15:end function

##  D Expected Log-Likelihood Implementation

The implementation details for calculating $\text{ELL}^{\text{reg}}$ based on
Equation (2.3.7) and $\text{ELL}^{\text{cls}}$ based on Equation (24) are
provided in Algorithms 3 and 4, respectively.

Algorithm 3 Calculate Expected Log-Likelihood (ELL) for Regression

1:

2:Observations $\mathbf{y}$

3:Predictive distribution $p(y_{i,d}^{r}\mid
f_{d}^{r}(\mathbf{(}X)))=\mathcal{N}(f_{d}^{r}(\mathbf{(}X)),\mathbf{\sigma}_{d}^{2})$

4:Noise variance $\sigma_{d}^{2}$

5:Sampled points $\mathbf{X}^{(j)}$, $j=1,\ldots,J$

6:function ExpectedLogLikelihoodRegression($\mathbf{y}$,
$\mathbf{\mu}_{\mathbf{f}_{:,d}^{r}}$,
$\mathbf{\Sigma}_{\mathbf{f}_{:,d}^{r}}$, $\mathbf{\sigma}_{d}^{2}$,
$\mathbf{X}^{(j)}$)

7: Initialize $ELL_{i,d}^{\text{reg}}\leftarrow 0$

8: for $j=1$ to $J$ do

9: Compute $\mathbf{\mu}_{\mathbf{f}_{:,d}^{r}}(x^{(j)}_{i})$ and
$\mathbf{\Sigma}_{\mathbf{f}_{:,d}^{r}}(x^{(j)}_{i})$

10: Calculate $A(x^{(j)}_{i})\leftarrow-0.5\left[\log
2\pi+\log\mathbf{\sigma}_{d}^{2}+\frac{(y_{i,d}^{r}-\mathbf{\mu}_{\mathbf{f}_{:,d}^{r}}(x^{(j)}_{i}))^{2}+\mathbf{\Sigma}_{\mathbf{f}_{:,d}^{r}}(x^{(j)}_{i})}{\mathbf{\sigma}_{d}^{2}}\right]$

11: Update $ELL_{i,d}^{\text{reg}}\leftarrow
ELL_{i,d}^{\text{reg}}+A(x^{(j)}_{i})$

12: end for

13: Normalize $ELL_{i,d}^{\text{reg}}\leftarrow-\frac{1}{2J}\times
ELL_{i,d}^{\text{reg}}$

14: return $ELL_{i,d}^{\text{reg}}$

15:end function

Algorithm 4 Calculate Expected Log Probability using Gauss-Hermite Quadrature

1:

2:Observations $\mathbf{y}$

3:Predictive distribution
$q_{\gamma}(\mathbf{f}_{k}^{c}\mid\mathbf{X})=\mathcal{N}\left(\mathbf{f}_{k}^{c}\mid\mathbf{\mu}_{f_{k}^{c}}(\mathbf{X}),\mathbf{\Sigma}_{f_{k}^{c}}(\mathbf{X})\right)$

4:Gauss-Hermite Quadrature nodes $\mathit{Loc}$

5:Gauss-Hermite Quadrature weights $\omega$

6:function ExpectedLogLikelihoodClassification($\mathbf{y}$, $f(\mathbf{x})$,
$\mathit{Loc}$, $\mathbf{\omega}$)

7: $\mathbf{y}\leftarrow 2\mathbf{y}-1$ $\triangleright$ Normalize
$\mathbf{y}$ to $\\{-1,1\\}$

8:
$\mathit{Loc}_{shifted}\leftarrow\sqrt{2\mathbf{\Sigma}_{f_{k}^{c}}(\mathbf{X})}\times\mathit{Loc}+\mathbf{\mu}_{f_{k}^{c}}(\mathbf{X})$
$\triangleright$ Shift and scale quadrature locations

9: $log\\_probs\leftarrow\log\Phi(\mathit{Loc}_{shifted}\cdot\mathbf{y})$
$\triangleright$ Compute log probabilities using normal CDF $\Phi$

10: $\mathbb{E}_{q_{\gamma}(\mathbf{f}_{k}^{c}\mid\mathbf{X})}[\log
P(\mathbf{y}\mid
f_{k}^{c}]\approx\frac{1}{J\sqrt{\pi}}\sum_{j=1}^{J}log\\_probs^{(j)}\times\omega^{(j)}$
$\triangleright$ Approximate expected log probability

11: return Approximated expected log probability

12:end function

## References

  * Abramowitz et al. (1988) Milton Abramowitz, Irene A Stegun, and Robert H Romer.  Handbook of mathematical functions with formulas, graphs, and mathematical tables, 1988.
  * Ali et al. (2017) M Usman Ali, Shahzad Ahmed, Javed Ferzund, Atif Mehmood, and Abbas Rehman.  Using pca and factor analysis for dimensionality reduction of bio-informatics data.  _arXiv preprint arXiv:1707.07189_ , 2017.
  * Alvear-Sandoval et al. (2019) Ricardo F Alvear-Sandoval, José L Sancho-Gómez, and Aníbal R Figueiras-Vidal.  On improving cnns performance: The case of mnist.  _Information Fusion_ , 52:106–109, 2019.
  * Balakrishnama and Ganapathiraju (1998) Suresh Balakrishnama and Aravind Ganapathiraju.  Linear discriminant analysis-a brief tutorial.  _Institute for Signal and information Processing_ , 18(1998):1–8, 1998.
  * Baudat and Anouar (2000) Gaston Baudat and Fatiha Anouar.  Generalized discriminant analysis using a kernel approach.  _Neural computation_ , 12(10):2385–2404, 2000.
  * Bishop and James (1993) Christopher M Bishop and Gwilym D James.  Analysis of multiphase flows using dual-energy gamma densitometry and neural networks.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 327(2-3):580–593, 1993.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Bro and Smilde (2014) Rasmus Bro and Age K Smilde.  Principal component analysis.  _Analytical methods_ , 6(9):2812–2831, 2014.
  * Cieslak et al. (2020) Matthew C Cieslak, Ann M Castelfranco, Vittoria Roncalli, Petra H Lenz, and Daniel K Hartline.  t-distributed stochastic neighbor embedding (t-sne): A tool for eco-physiological transcriptomic analysis.  _Marine genomics_ , 51:100723, 2020.
  * Cunningham and Yu (2014) John P Cunningham and Byron M Yu.  Dimensionality reduction for large-scale neural recordings.  _Nature neuroscience_ , 17(11):1500–1509, 2014.
  * Deng (2012) Li Deng.  The mnist database of handwritten digit images for machine learning research.  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.
  * Dey and Salem (2017) Rahul Dey and Fathi M Salem.  Gate-variants of gated recurrent unit (gru) neural networks.  In _2017 IEEE 60th international midwest symposium on circuits and systems (MWSCAS)_ , pages 1597–1600. IEEE, 2017.
  * Ek and Lawrence (2009) Carl Henrik Ek and PHTND Lawrence.  _Shared Gaussian process latent variable models_.  PhD thesis, Oxford Brookes University Oxford, UK, 2009.
  * Fard et al. (2023) Reza Saadati Fard, Navid Ziaei, and Ali Yousefi.  Latent dynamical model to characterize brain network-level rhythmic dynamics.  In _2023 45th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_, pages 1–5. IEEE, 2023.
  * Fisher (1988) R. A. Fisher.  Iris.  UCI Machine Learning Repository, 1988.  DOI: https://doi.org/10.24432/C56C76.
  * Gao et al. (2010) Xinbo Gao, Xiumei Wang, Dacheng Tao, and Xuelong Li.  Supervised gaussian process latent variable model for dimensionality reduction.  _IEEE transactions on systems, man, and cybernetics, Part B (Cybernetics)_ , 41(2):425–434, 2010.
  * Gilks et al. (1995) Walter R Gilks, Sylvia Richardson, and David Spiegelhalter.  _Markov chain Monte Carlo in practice_.  CRC press, 1995.
  * Gu and Dao (2023) Albert Gu and Tri Dao.  Mamba: Linear-time sequence modeling with selective state spaces.  _arXiv preprint arXiv:2312.00752_ , 2023.
  * Heller and David (2022) Charles R Heller and Stephen V David.  Targeted dimensionality reduction enables reliable estimation of neural population coding accuracy from trial-limited data.  _PloS one_ , 17(7):e0271136, 2022.
  * Hensman et al. (2013) James Hensman, Nicolo Fusi, and Neil D Lawrence.  Gaussian processes for big data.  _arXiv preprint arXiv:1309.6835_ , 2013.
  * Hensman et al. (2015) James Hensman, Alexander Matthews, and Zoubin Ghahramani.  Scalable variational gaussian process classification.  In _Artificial Intelligence and Statistics_ , pages 351–360. PMLR, 2015.
  * Jiang et al. (2012) Xinwei Jiang, Junbin Gao, Tianjiang Wang, and Lihong Zheng.  Supervised latent linear gaussian process latent variable model for dimensionality reduction.  _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , 42(6):1620–1632, 2012.
  * Kingma et al. (2019) Diederik P Kingma, Max Welling, et al.  An introduction to variational autoencoders.  _Foundations and Trends® in Machine Learning_ , 12(4):307–392, 2019.
  * Koh et al. (2023) Tze Hui Koh, William E Bishop, Takashi Kawashima, Brian B Jeon, Ranjani Srinivasan, Yu Mu, Ziqiang Wei, Sandra J Kuhlman, Misha B Ahrens, Steven M Chase, et al.  Dimensionality reduction of calcium-imaged neuronal population activity.  _Nature Computational Science_ , 3(1):71–85, 2023.
  * Lawrence and Hyvärinen (2005) Neil Lawrence and Aapo Hyvärinen.  Probabilistic non-linear principal component analysis with gaussian process latent variable models.  _Journal of machine learning research_ , 6(11), 2005.
  * Liu et al. (2018) Haitao Liu, Jianfei Cai, and Yew-Soon Ong.  Remarks on multi-output gaussian process regression.  _Knowledge-Based Systems_ , 144:102–121, 2018.
  * Ma and Dai (2011) Shuangge Ma and Ying Dai.  Principal component analysis based methods in bioinformatics studies.  _Briefings in bioinformatics_ , 12(6):714–722, 2011.
  * McInnes et al. (2018) Leland McInnes, John Healy, and James Melville.  Umap: Uniform manifold approximation and projection for dimension reduction.  _arXiv preprint arXiv:1802.03426_ , 2018.
  * Murray and Adams (2010) Iain Murray and Ryan P Adams.  Slice sampling covariance hyperparameters of latent gaussian models.  _Advances in neural information processing systems_ , 23, 2010.
  * Salimbeni and Deisenroth (2017) Hugh Salimbeni and Marc Deisenroth.  Doubly stochastic variational inference for deep gaussian processes.  _Advances in neural information processing systems_ , 30, 2017.
  * Seeger (2004) Matthias Seeger.  Gaussian processes for machine learning.  _International journal of neural systems_ , 14(02):69–106, 2004.
  * Snelson and Ghahramani (2005) Edward Snelson and Zoubin Ghahramani.  Sparse gaussian processes using pseudo-inputs.  _Advances in neural information processing systems_ , 18, 2005.
  * Tipping and Bishop (1999) Michael E Tipping and Christopher M Bishop.  Probabilistic principal component analysis.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 61(3):611–622, 1999.
  * Titsias (2009) Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In _Artificial intelligence and statistics_ , pages 567–574. PMLR, 2009.
  * Titsias and Lawrence (2010) Michalis Titsias and Neil D Lawrence.  Bayesian gaussian process latent variable model.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 844–851. JMLR Workshop and Conference Proceedings, 2010.
  * Wang et al. (2014) Wei Wang, Yan Huang, Yizhou Wang, and Liang Wang.  Generalized autoencoder: A neural network framework for dimensionality reduction.  In _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , pages 490–497, 2014.
  * Williams and Seeger (2000) Christopher Williams and Matthias Seeger.  Using the nyström method to speed up kernel machines.  _Advances in neural information processing systems_ , 13, 2000.
  * Wong and Li (2000) Chun Shan Wong and Wai Keung Li.  On a mixture autoregressive model.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 62(1):95–115, 2000.
  * Yu et al. (2006) Shipeng Yu, Kai Yu, Volker Tresp, Hans-Peter Kriegel, and Mingrui Wu.  Supervised probabilistic principal component analysis.  In _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 464–473, 2006.
  * Yu et al. (2019) Yong Yu, Xiaosheng Si, Changhua Hu, and Jianxun Zhang.  A review of recurrent neural networks: Lstm cells and network architectures.  _Neural computation_ , 31(7):1235–1270, 2019.
  * Zebari et al. (2020) Rizgar Zebari, Adnan Abdulazeez, Diyar Zeebaree, Dilovan Zebari, and Jwan Saeed.  A comprehensive review of dimensionality reduction techniques for feature selection and feature extraction.  _Journal of Applied Science and Technology Trends_ , 1(2):56–70, 2020.
  * Zhang et al. (2017) Jiayuan Zhang, Ziqi Zhu, and Jixin Zou.  Supervised gaussian process latent variable model based on gaussian mixture model.  In _2017 International Conference on Security, Pattern Analysis, and Cybernetics (SPAC)_ , pages 124–129. IEEE, 2017.
  * Zhong and Enke (2017) Xiao Zhong and David Enke.  Forecasting daily stock market return using dimensionality reduction.  _Expert systems with applications_ , 67:126–139, 2017.
  * Ziaei et al. (2023) Navid Ziaei, Reza Saadatifard, Ali Yousefi, Behzad Nazari, Sydney S Cash, and Angelique C Paulk.  Bayesian time-series classifier for decoding simple visual stimuli from intracranial neural activity.  In _International Conference on Brain Informatics_ , pages 227–238. Springer, 2023.