\section{Introduction}\label{intro}
There is increasing interest in the problem of nonparametric regression with high-dimensional predictors. In a variety of fields, from computer science to environmental science, one often encounters high-dimensional data (e.g., ‘point cloud data’) perturbed by high-dimensional noise but cantering around some lower-dimensional manifolds. In more precise mathematical terms, manifolds are topological spaces equipped with some differential/smooth structure, the geometry of which is in general different from the usual Euclidean geometry. Naively applying traditional multivariate analysis to manifold-valued data that ignores the geometry of the space can potentially lead to highly misleading predictions and inferences. Gaussian Processes (GP) are among the powerful tools in statistics and machine learning for regression with Euclidean predictors\citep{Rasmussen2004}. However, if the predictors locate on a manifold, traditional smoothing or modelling methods like GP that do not respect the intrinsic geometry of the space and the boundary constraints, produce poor results \citep{extrinsicGP, niu2019}. In this work, we focus on the following model throughout the paper,
\begin{align} \label{eqn:reg}
y_i = f(s_i) +\epsilon_i, \ \ \  \epsilon_i \sim \mathcal{N}(0, \sigma_{noise}^2)
\end{align}
where  $f:\mathbb{M} \rightarrow \mathbb{R}$ is an unknown regression function. $y_i \in \mathbb{R}$ is an observed response variable. $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is an observed predictor on the embedded manifold in $\mathbb{R}^p$, $s_i = (s_{i}^1,\ldots,s_{i}^p)$. $\mathbb{M}$ is a complete and orientable Riemannian manifold embedded in $\mathbb{R}^p$. We would like to do inferences on how the output $y$ varies with the input $s$, including predicting $y$ values at new locations not represented in the training dataset. A GP prior can be assigned to $f$ with a covariance function $k(s_i,s_j)$. The choice of $k(s_i,s_j)$ has a fundamental impact on the results. Most common choices of covariance function such as the squared exponential kernel and Matern kernel https://www.overleaf.com/project/5ff041180b93f20680f0d6ecdepend critically on the Euclidean distance between $s_i$ and $s_j$ and ignore the intrinsic geometry of $\mathbb M$.

One of the paramount challenges in developing GP models on manifold is the difficulty in specifying the covariance structure via constructing valid and computable covariance kernels on manifolds.  \cite{niu2019} proposed to use the heat kernel to construct the intrinsic Gaussian process (In-GP) on complex constrained domain. Although the heat kernel provides a natural and canonical choice theoretically, it is analytically intractable to directly evaluate \citep{hsu1988}. \cite{Castillo2014} instead proposed to use randomly rescaled solutions of the heat equation to define a valid covariance kernel for reasonably broad classes of compact manifolds. But their proposed heat kernels are computationally intractable. Alternatively, it can be estimated as the transition density of the Brownian Motion (BM) on the manifold \citep{hsu2008}.  By using the heat kernel, In-GP can accommodate spatial domains arising as complex subsets of Euclidean space. \cite{niu2019} has demonstrated the application of In-GP in modelling the chlorophyll concentration levels in Aral Sea \citep{wood} and smooth functions on Swiss roll.  It respects the potentially complex boundary or interior conditions as well as the intrinsic geometry of the spaces. However \cite{niu2019}'s method is only applicable when the geometry of the manifold is known and the dimension of the dataset space is small such as the subspace of $\mathbb{R}^2$ or $\mathbb{R}^3$.  In this paper, we develop a computational inferential framework able to estimate the heat kernel and construct the In-GP on a unknown manifold embedded in the high dimensional point cloud.

\begin{figure}[htb]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{map.pdf}
    \caption{ \label{fig:mapping}
    \footnotesize{
    { An illustration of the parameterisation of the manifold: $\phi$ maps the chart into the embedding space }}
    }
\end{figure}
In statistics and machine learning community, we study probabilistic latent variable models (LVMs) as the parameterisation or embeddings of the manifold in the observation space. The deterministic LVM is first studied to describe curved surface By \cite{gauss1828}. He shown a low dimensional surface embedded in a high dimensional Euclidean space can be well described through a mapping function.  The low dimensional representation of the surface is known as the chart which corresponds to the latent space in probabilistic LVMs. Intuitively, the chart provides a distorted view of the surface. An illustration is shown in Fig. \ref{fig:mapping}. $x^1$ and $x^2$ coordinates in the chart represent the radius and width of the Swiss roll in $\mathbb{R}^3$.  Points in the chart can be mapped to the embedded space (Swiss roll)  through the function $\phi$. In other words, the Swiss roll can be parameterised by the chart. Measurements on the manifold can be computed in the chart locally, and integrated to provide global measures. This gives rise to the definition of a local inner product, known as a Riemannian metric tensor. A Riemannian metric $\mathbf{g}$ on a manifold $\mathbb{M}$ is symmetric and positive definite matrix which defines a smoothly varying inner product in the tangent space.
\begin{align} \label{eqn:deter_metric}
\mathbf{g} = \mathcal{J}^T \mathcal{J}, \ \ \  \mathcal{J}_{i,j} = \frac{\partial \phi_i }{ \partial  x^j}
\end{align}
where $\mathcal J$ denotes the Jacobian of the $\phi$. In \cite{niu2019} , both the mapping function $\phi$ and the metric tensor $\mathbf{g}$ are known.

In modern big data analysis, the data in the point cloud, often high dimensional, are not directly observed on the manifolds, instead they are observed in a potentially high-dimensional ambient space but concentrate around some unknown lower dimensional structure. The mapping function $\phi$ is unknown. One needs to learn this geometry before utilizing it for inference. Most recent research in \cite{dunson2020diffusion} tackled this problem by approximating the heat kernel with finitely-many eigenpairs of the Graph Laplacian. Since they are using a compact operator to approximate an unbounded operator, hence the eigenvalues are hard to recover and need many points in the point cloud. \cite{dunson2020diffusion}  has shown 3000 points are required in order to approximate the heat kernel of a simple one-dimensional unit circle or spiral from the point cloud in $\mathbb{R}^2$. In addition, they  did not consider the boundary condition. In many applications such as image processing, the dimension of the data in the point cloud  easily scales up to thousands. Neither  \cite{niu2019}   or \cite{dunson2020diffusion} work for this type of problems. There are still some critical gaps in current practice in particular in lacking of robust methods for carrying out intrinsic statistical inference, and largely missing effective models for regressions with manifold-valued data embedding in a point cloud.  In this paper, we will use the Bayesian dimension reduction method to learn the metric tensor of the underlying manifold of the point cloud and define its boundaries. The heat kernel can be estimated as the transition density of the BM on the learned manifold. We have proved the BM transition density estimation is coordinate independent. Different metric learning methods such as GPLVM (Gaussian Processes Latent Variable Model) and Bayesian GPLVM are compared. This framework allows us to build the In-GP models on the implicit manifold embedded in the point cloud. It can incorporate fully the intrinsic geometry of the implicit manifold for inference while respecting the potentially complex interior structure and boundary.

\section{Intrinsic Gaussian process on point cloud}
\label{sec:GPmani}
We propose to construct the intrinsic Gaussian processes on implicit manifold by using the heat kernel as the covariance kernel. Let $\mathbb{M}$ be a $q$-dimensional complete and orientable Riemannian manifold with the Riemannian metric $\mathbf{g}$, $\partial \mathbb{M}$ its boundary that is continuous. A heat kernel of $\mathbb{M}$ is a smooth function $K(s_0,s,t)$ on $ \mathbb{M} \times \mathbb{M} \times \mathbb{R}^{+}$ that satisfies the heat equation:
\begin{align*}
\frac{\partial}{\partial t}K_{heat}(s_0,s,t)=\frac{1}{2}\Delta_s K_{heat}(s_0,s,t),\qquad  s_0, s\in \mathbb{M},
\end{align*}
where $\Delta_s$ is the Laplacian-Beltrami operator on $\mathbb{M}$, and $t>0$ is the diffusion time. The heat kernel satisfies the initial condition $\lim_{t\rightarrow 0}K_{heat}(s_0,s,0)=\delta(s_0,s)$ with $\delta$ the Dirac delta function. The heat kernel becomes unique when we impose a suitable condition along the boundary $\partial \mathbb{M}$, such as the Neumann boundary condition. $ \partial K / \partial \mathbf{n} = 0$ along $\partial \mathbb{M}$ where $\mathbf{n}$ denotes a normal vector of $\partial \mathbb{M}$.  The Neumann boundary condition can be expressed as no heat transfer across the boundary. Here we use $K_{heat}^t$ to denote the heat kernel of arbitrary manifold. In the local coordinate (chart), the Laplacian-Beltrami operator can be written as
 \begin{align}
\label{eqn:delta_s}
\Delta_s K_{heat}^t = \frac{1}{\sqrt{G}} \frac{\partial}{\partial x^j} \left (    \sqrt{G} \mathbf{g}^{ij} \frac{\partial K_{heat}^t}{\partial x^i}  \right ),
\end{align}
where $G$ is the determinant of metric $\mathbf{g}$ and $\mathbf{g}^{ij}$ is the $(i,j)$ element of its inverse. Take the special case $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$, $\mathbf{g}$ become an identity matrix. The Laplace Beltrami operator become the Laplace operator.    The heat kernel has a closed form corresponding to a time-varying Gaussian function:
\begin{align*}
K_{heat}(\mathbf{s}_0,\mathbf{s},t)
=\frac{1}{(2\pi t)^{q/2}}\,
  \exp\left\{-\frac{||\mathbf{s}_0-\mathbf{s}||^2}{ 2t }\right\}, \; \mathbf{s}\in \mathbb R^q.
\end{align*}
The diffussion time $t$ controlling the rate of decay of the covariance. By varying the time parameter, one can vary the bumpiness of the realisations of the in-GP.

We now consider the heat kernel as the covariance function of the In-GP. Let $\mathcal{D} =  \{ (s_i, y_i), i = 1, . . . , n \} $, be the data, with $n \geq 1$ the number of labeled observations. $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is the predictor and $y_i \in \mathbb{R}$ is the corresponding response variable. Suppose we are also given an unlabeled dataset $\mathcal{V} = \{ s_i, i=n+1 . . . ,n+v \} $ where $v \geq 1$. Consider the regression model in \ref{eqn:reg}, we would like to do inferences on $f$ with the labeled dataset $\mathcal{D}$ and predict $y$ values at the unlabeled dataset $\mathcal{V} $. Under an in-GP prior for the unknown regression function as $f \sim GP(0,K_{heat}^t(.,.) )$, we have
\begin{align*}
p( \mathbf{f} | s_1, s_2, ..., s_n) \sim \mathcal{N} (0, \Sigma_{\mathbf{f} \mathbf{f}}),
\end{align*}
where $\mathbf{f}$ is the discretization of $f$ over $s_1,s_2,\cdot\cdot\cdot,s_n$ so that $f_i = f(s_i)$. $\Sigma_{\mathbf{f} \mathbf{f}} \in \mathbb{R}^{n\times n}$ is the covariance matrix induced from the heat kernel. The $(i,j)$ entries of $\Sigma_{\mathbf{f}\f}$ corresponds to $\Sigma_{\mathbf{f}\f _{i,j}} = \sigma_h^2 K_{heat}^t(s_i,s_j)$. We introduce the rescaling hyperparameter $\sigma_h^2$ to add extra flexibility of the heat kernel. Let ${\bf f}_*$ be a vector of values of $f(\cdot)$ at the unlabelled grid points $\mathcal{V}$. Under the regression model in equation \ref{eqn:reg} and the GP prior, we have the joint distribution of ${\bf y} $ and ${\bf f}_*$ is:
\begin{align}
p( {\text {\bf y}}, {\text {\bf f}}_*) = \mathcal{N} \left ( 0,
 \left [ \begin{array}{cc}
 \Sigma_{{\text {\bf f}} {\text {\bf f}}} + \sigma_{noise}^2 I_n &\Sigma_{ {\text {\bf f}}  {\text {\bf f}}_*} \\[0.3em]
 \Sigma_{ {\text {\bf f}}_* {\text {\bf f}}} &\Sigma_{ {\text {\bf f}}_* {\text {\bf f}}_*} \\[0.3em]
                \end{array} \right ]
                \right ),
\end{align}
where $\Sigma_{{\text {\bf f}}_*{\text {\bf f}}}$ is the covariance matrix for training and unlabelled grid data points. The predictive distribution is derived by marginalising out {\bf f}:
\begin{align}
p( {\text {\bf f}}_*| \bm y) = \int p( {\text {\bf f}}_*{\text {\bf f}}| \bm y) d{\text {\bf f}} = \mathcal N \left ( \Sigma_{ {\text {\bf f}}_*{\text {\bf f}}}  \left(\Sigma_{{\text {\bf f}}{\text {\bf f}}}  + \sigma_{noise}^2 I \right)^{-1}\bm y , \ \  \Sigma_{{\text {\bf f}}_*{\text {\bf f}}_*}-\left(\Sigma_{{\text {\bf f}}{\text {\bf f}}} +\sigma_{noise}^2 I \right)^{-1} \Sigma_{{\text {\bf f}}{\text {\bf f}}_*} \right)
\end{align}

The main difficulty of using In-GP is the analytic intractability of evaluating $K_{heat}^t$. The analytical solution of heat kernel only exists for very few special cases such as Euclidean space. Solving the heat equation for an arbitrary manifold is nontrivial. \cite{niu2019} estimate the heat kernel as the BM transition density by simulating BM paths on a known manifold. In the case of point cloud, both the metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$ are unknown. In this work, we will use the probabilistic latent variable models to learn the metric $\mathbf{g}$  and defining the boundary $\partial \mathbb{M}$ based on the uncertainty quantification of the mapping. The heat kernel can then be estimated by simulating BM paths on the learned manifold.

\section{Learning metric and boundary}
Let $\mathcal S = \{ s_i |  i=1,\cdot \cdot \cdot, n+v \}$,  $s_i \in \mathbb{R}^p$, include both predictors in the labeled dataset $\mathcal{D}$ and unlabelled dataset $\mathcal{V}$. We have $\mathcal{S} \in \mathbb{R}^{ (n+v) \times p }$. Suppose we perform probabilistic nonlinear dimensionality reduction by defining a latent variable model that introduces a set of latent (unobserved) variables $\mathcal{X}= \{x_i | i= 1,\cdot \cdot \cdot, n+v \}$, $x_i \in \mathbb{R}^q$ and $q < p$.  $\mathcal{X} \in \mathbb{R}^{(n+v)\times q}$. $\mathcal X$ is related to $\mathcal S$ which is observed in a higher dimensional space. The latent space is dominated by a prior distribution $p( \mathcal X )$ which induces a distribution over $\mathcal S$ under the assumption of a probabilistic mapping
\begin{align} \label{eqn:mapmodel}
s_{i}^j = \phi^j(x_i) + e_{i}^j,
\end{align}
where $x_i \in \mathbb{R}^q$ is the $i^{th}$ latent point associated with $s_i \in \mathbb{R}^p$ in the original space, $i<n$. $j$ is the index of the features (dimensions) of $s$ in the observation space and $j<p$. $e_{i}^j$ is a noise term. The noise is Gaussian distributed $e_i^{j} \sim \mathcal N(0,\beta^2)$, $\beta^2$ is the variance. An illustration is given in Fig. \ref{fig:mapping} where $q=2$ and $p=3$. If the mapping $\phi$ is linear and the prior $p( \mathcal X)$ is Gaussian, this model is known as probabilistic principal component analysis \citep{tipping1998}. In this work, we do not restrict to this linear assumption and consider some nonlinear dimensionality reduction methods such as Gaussian Process Latent Variable Model (GP-LVM) and Bayesian Gaussian Process Latent Variable Model (BGP-LVM).

When the mapping $\phi$ in Fig.\ref{fig:mapping} is differentiable, it can be interpreted as the mapping between the chart (or latent space) $\mathbb{R}^q$ and the embedding space $\mathbb{M} \subset \mathbb{R}^p$.  The natural Riemannian metric of the given model can be computed as in equation \ref{eqn:deter_metric}. In the case of probabilistic LVMs where the conditional probability over the Jacobian follow a Gaussian distribution, this naturally induces a distribution over the local metric tensor $\mathbf{g}$. We denote the distribution of the Jacobian as $\mathbf{J}$. Assuming independent rows of $\mathbf J$
\begin{align}
p(\mathbf J | \mathcal{S}, \mathbf{f}) = \prod_{j=1}^p \mathcal{N}( \mu_{\mathbf{J}}^j, \Sigma_{\mathbf J}  )
\end{align}
The expression of the mean $\mu_{\mathbf J}$ and variance $\Sigma_{\mathbf J} $ of the Jacobian are model specific and given in section \ref{gplvm} and \ref{bgplvm}. This resulting the metric $\mathbf{g}$ follow a non-central Wishart distribution \citep{anderson1946}
\begin{align} \label{eqn:metric}
\mathbf{g} \sim \mathcal{W}_q \left(  p, \Sigma_{\mathbf{J}}, \E(\mathbf{J}^T)  \E(\mathbf{J})  \right).
\end{align}
From this distribution, the expected metric tensor can be computed as
\begin{align} \label{eqn:exmetric}
\E( \mathbf{g} ) =  \mathcal{G} = E( \mathbf{J}^T) \E( \mathbf{J} ) + p \Sigma_{\mathbf{J}}.
\end{align}
We denote the expectation by $\mathcal{G}$. We also need the gradient of the expected metric to simulate BM as in section \ref{trandens}.
\begin{align}
\frac{ \partial \E[ \mathbf{g}]}{\partial x^l} = \frac{\partial \mathcal{G}}{\partial x^l} = \frac{\partial \E[\mathbf{J}^T] } { \partial x^l }  \E[\mathbf{J}] + \E[\mathbf{J}^T] \frac{ \partial  \E[\mathbf{J}]  }{\partial x^l}  + p \frac{ \partial \Sigma_\mathbf{J} }{ \partial x^l }
\end{align}
The boundary of the learned manifold can be defined by checking $Var( \phi(x) | x)$, the variance of the mapping at $x$. Any $x$ outside of the boundary has $Var( \phi(x) | x ) > \alpha$, $\alpha$ is a predefined threshold. Both $\mathcal{G}$ and $\partial \mathbb{M}$  are used to construct BM on manifold in section \ref{}.   ???
\begin{align}\label{eqn:bound}
\partial \mathbb{M} = \{  x \in \mathbb{R}^q | Var( \phi(x) | x) = \alpha   \}.
\end{align}

\subsection{GPLVM}\label{gplvm}
Gaussian Process Latent Variables Model \citep{lawrence2005} is a nonlinear dimensionality reduction method. It defines a generative mapping from the latent space $x\in \mathbb{R}^q$ to observation space $s \in \mathbb{M} \subset R^p$ that is governed by Gaussian processes. GPs priors are placed on $\phi$ in equation \ref{eqn:mapmodel}. The GPs are taken to be independent across the dimensions of the observation space. This leads to the formulation of GP-LVM. Given the construction outlined above, the joint probability over the observed data $\mathcal{S}$ and all latent variables $\mathcal{X}$ is written as follows:
\begin{align}\label{eqn:gplvmJoint}
p(\mathcal{S} \ \Phi \ \mathcal{X} | \beta) =   p(\mathcal{S} | \Phi \beta ) p(\Phi | \mathcal{X}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j, \beta ) p( \bm{\phi}_{:}^j | \mathcal{X}) p(\mathcal{X})
\end{align}
 where $\Phi = \{ \bm{\phi}_i | i=1,\cdot\cdot\cdot, n+v \} $, $\Phi \in \mathbb{R}^{(n+v)\times p}$ and $ \bm{\phi}_i^j = \phi(x_i)^j$.  $\mathbf{s}_{:}^j$ represents the $j_{th}$ dimension of every points in $\mathcal S$ and $\mathcal S \in \mathbb{R}^{ (n+v)\times p}$. The likelihood $p( \mathcal S | \mathcal{X} )$ is computed by marginalising out $\Phi$ and optimising the latent variables $\mathcal{X}$.
\begin{align}
p(\mathcal{S} | \mathcal{X} ) = \prod_{j=1}^p \mathcal N( \mathbf{s}_{:}^j, K_{\mathcal{X} \mathcal{X}} + \beta^2 I )
\end{align}
$K_{\mathcal{X} \mathcal{X}} $ is the $(n+v)\times (n+v)$ covariance matrix defined by the squared exponential kernel. Since $x \in \mathcal{X}$ is a latent variable, we can assign it a prior density given by the standard normal density. \cite{lawrence2005} estimate all latent variables $\mathcal X$ and kernel hyper-parameters of GP-LVM by following the maximum likelihood estimation. If the covariance kernel is differentiable, the derivative of a GP is again a GP \citep{Rasmussen2004}. This property allows predictions about the derivative of GP. The Jacobian $\mathbf{J}$ of the GP-LVM mapping can be computed as the partial derivative $\frac{\partial \phi_*}{ \partial x^l}$ with respect to the $l_{th}$ dimension for any point $x_*$  in the latent space.
\begin{align} \label{eqn:Jacobian}
\mathbf{J}^{T} = \frac{\partial \phi_*}{ \partial x} = \begin{bmatrix}
 \frac{ \partial \phi(x_*)^1 }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^1}   \\
  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot \\
  \frac{ \partial \phi(x_*)^1 }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^q}
\end{bmatrix}
\end{align}
where $\frac{\partial \phi_*}{ \partial x}$ is a $q \times p$ matrix. Considering the independence across the dimensions of the observation space, the joint distribution of the $j_{th}$ dimension of the mapping $\phi$ and the $j_{th}$ column of the Jacobian can be written as
\begin{align}
\begin{bmatrix}
 \phi( \mathcal{X})^j\\
 \frac{ \partial \phi(x_*)^j }{\partial x}
\end{bmatrix}, \
\sim \mathcal N \left( 0, \begin{bmatrix}
K_{\mathcal{X},\mathcal{X}} & \partial K_{\mathcal{X},*} \\
 \partial K_{\mathcal{X},*}^T & \partial^2 K_{*,*}
\end{bmatrix}  \right )
\end{align}
The expressions of $K_{\mathcal{X},\mathcal{X}}$ ,  $\partial K_{\mathcal{X},*}$ , and  $\partial^2 K_{*,*} $ are given in the appendix ??. The GP-LVM model provides an explicit mapping from the latent space to the observed space. This mapping defines the support of the observed data $\mathcal{S}$ as a $q$ dimensional manifold embedded into $\mathbb{R}^p$. The distribution of the Jacobian of the GP-LVM mapping is the product of $p$ independent Gaussian distributions (one for each dimension of the observation space). For a latent point $x_*$, the  distribution of the Jacobian takes this form
\begin{align}
p(  \mathbf{J} | \mathcal{X}, \mathcal{S}) &= \prod^p_{j=1} \mathcal N( \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \partial K^T_{\mathcal{X},*} K^{-1}_{\mathcal{X},\mathcal{X}} \mathbf{s}_{:}^j  \ , \partial^2K_{*,*} - \partial K^T_{\mathcal{X},*} K_{\mathcal{X},\mathcal{X}}^{-1} \partial K_{\mathcal{X},*}  ) \nonumber
\end{align}
From this distribution, the expected metric tensor can be computed as in equation \ref{eqn:exmetric}. The boundary $\partial M$ is defined by $Var( \phi(x_*) | x_*)$. For any $x_*$ in the latent space,
\begin{align}\label{eqn:varb}
Var(\phi(x_*) | x_*) = K_{*,*} - K_{*,\mathcal{X}}K_{\mathcal{X},\mathcal{X}}^{-1} K_{\mathcal{X},*}.
\end{align}

\subsection{Bayesian-GPLVM }\label{bgplvm}
The GP-LVM based models rely on MAP training procedures for optimising the latent inputs and the hyperparameters \citep{lawrence2005}. However, the standard GP-LVM does not mariginalise out the latent inputs $\mathcal{X}$ implies that it could be sensitive to overfitting \citep{damianou2016}. The marginal likelihood takes this form
\begin{align*}
p(\mathcal S)=\int p(\mathcal{S} \mid \mathcal{X}) p(\mathcal{X}) d \mathcal{X}
\end{align*}
The  integral is intractable as the inputs $\mathcal{X}$ of the latent variable model appears nonlinearly in the inverse of the covariance matrix. Bayesian Gaussian process latent variable model (BGP-LVM, \cite{titsias2010}) introduce a variational inference framwork for training the latent variable model and perform the nonlinear dimensional reduction. It allows variationally integrate out the input variables and compute a lower bound on the exact marginal likelihood of the nonlinear latent variable model. The maximisation of the variational lower bound provides a Bayesian training procedure that is robust to overfitting.

The key to make the variational Bayes approach tractable is the application of variational inference in an expanded probability model where the GP prior on $\phi$ is augmented to include auxiliary variables. More specifically, we expand the joint probability model in equation \ref{eqn:joint} by including $m$ extra samples (inducing points) of the GP latent mapping $\phi(x)$, so that $u_i = \phi(x_{ui}) \in \mathbb{R}^p$ is such a sample \citep{lawrence2007}. These inducing points are denoted by $\mathcal{U} = \{u_i| i=1,\cdot\cdot\cdot,m   \}$, $\mathcal{U} \in \mathbb{R}^{m\times p}$ and constitute latent function evaluations at a set of pseudo-inputs $\mathcal{X}_u = \{ x_{ui} | i=1,\ldots,m \} $, $\mathcal{X}_u \in \mathbb{R}^{m \times q}$. The inducing inputs $\mathcal{X}_u$ are not random variables. $\Phi$ is defined as in equation \ref{eqn:gplvmJoint}. The augmented joint probability and the marginal likelihood can be written as
\begin{align}\label{eqn:full-marginal}
p(\mathcal{S} \ \Phi \ \mathcal{U} \  \mathcal{X} | \beta) &=   p(\mathcal{S} | \Phi ) p(\Phi | \mathcal{U} \mathcal{X}) p(\mathcal {U}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j) p( \bm{\phi}_{:}^j |  u_{:}^j, \mathcal{X}) p(u_{:}^j) p(\mathcal{X}) \\
p(\mathcal{S}) &= \int  \int  \int  p(\mathcal{S}, \Phi, \mathcal{U}, \mathcal{X}) \ d \mathcal{U}     \  d\Phi \ d \mathcal{X}
\end{align}
Using this variational inference, $log\left(p(\mathcal{S} )\right)$ can be expressed as the Jensen’s lower bound. We can analytically maximize the lower bound. The full technical details of the lower bound are given in the appendix ??.  In contrast to the standard GPLVM, Bayesian GPLVM variationally integrate the latent variables $\mathcal{X}$ out and obtain extra regularisation. The Jacobian have the same shape as in equation \ref{eqn:Jacobian}. For a latent point $x_*$ the distribution the distribution of the Jacobian for the Bayesian GP-LVM mapping takes the form:
\begin{align}
p(  \mathbf{J} | \mathcal{X},\mathcal{U},\mathcal{S}) &= \prod^p_{j=1} \mathcal N( \  \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \  \partial K^T_{\mathcal{X}_u,*} K^{-1}_{\mathcal{X}_u,\mathcal{X}_u} u_{:}^j , \  \partial^2K_{*,*} - \partial K^T_{\mathcal{X}_u,*} K_{\mathcal{X}_u,\mathcal{X}_u}^{-1} \partial K_{\mathcal{X}_u,*}  ) \nonumber
\end{align}
The expressions of $ \partial K_{\mathcal{X}_u,*} $ and $\partial^2K_{*,*} $ are given in the appendix??. From this distribution, the expected metric tensor can be computed as in equation \ref{eqn:exmetric}.

\section{Estimate heat kernel as BM transition density}
The construction of the heat kernel associated with the Laplace-Beltrami operator is not a trivial task and belongs to the fields of partial differential equations and differential geometry. It is carried out in great detail in \cite{chavel1984}. The Laplace Beltrami operator in equation \ref{eqn:delta_s} is also the infinitesimal generator of Brownian Motion on the manifold \citep{hsu2008}.  Let
$\mathbf{S}_t$ be a Brownian motion with transition density $K_{heat}^t$. This transition density satisfies the heat equation. In this section, we will estimate the heat kernel as the BM transition density to avoid solving the heat equation on manifold.

\subsection{ Simulating Brownian motion on manifolds }
\label{trandens}
From last section, we have learned the parameterisation of the unknown manifold $\mathbb{M}$ as the latent space $\mathcal{X}$, the distribution of the associated metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$. In order to estimate the BM transition density on $\mathbb{M}$, we first need to simulate BM trajectories. Simulating a sample path of BM on $\mathbb{M} \subset \mathbb{R}^p$ is equivalent to simulating a stochastic process in the latent space in $\mathbb{R}^q$, $q<p$. An illustration example is shown in Fig. \ref{fig:BM}.

\begin{figure}[htb]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{bmlatent1.pdf}
    \caption{ \label{fig:BM}
    \footnotesize{
    { BM on $\mathbb{M}$ and its equivalent stochastic process in chart in $\mathbb{R}^2$.  $\phi: \mathbb{R}^2 \rightarrow \mathbb{M}$ is a parametrisation of $\mathbb{M}$.  }}
    }
\end{figure}

The BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations in the Ito form \citep{hsu1988,hsu2008}. We will use the expectation $\mathcal{G} = \E[\mathbf{g}]$ to construct this SDE.
\begin{align}
\label{eqnBM}
dx^i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x^j} \left(  \mathcal{G}^{ij}G^{1/2} \right) dt + \left( \mathcal{G}^{-1/2} dB(t)\right)_i,
\end{align}
where $x^i$ represent the $i_th$ dimension of the latent space (chart). $\mathcal{G}$ is the metric tensor of $\mathbb{M}$, $G$ is the determinant of $\mathcal{G}$ and $B(t)$ represents an independent BM in the Euclidean space. In order to do forward simulation, the discrete form of equation (\ref{eqnBM}) is derived in equation \eqref{disBM}.
\begin{align}
\label{disBM}
x^i(t) &= x^i(t-1) + \frac{1}{2} \sum^{q}_{j=1} \left(  -\mathcal{G}^{-1} \frac{\partial \mathcal{G}}{\partial x^j} \mathcal{G}^{-1} \right)_{ij} \Delta t + \frac{1}{4} \sum_{j=1}^q(\mathcal{G}^{-1})_{ij}tr(\mathcal{G}^{-1}\frac{\partial \mathcal{G}}{\partial x^j}) \Delta t + \left( \mathcal{G}^{-1/2} dB(t)\right)_i \nonumber \\
         &= \mu( x^i(t-1),\Delta t)_i + \left( \sqrt{\Delta t} \mathcal{G}^{-1/2} z^q\right)_i,
\end{align}
where $\Delta t$ is the diffusion time of the BM simulation step and $z^q$ is a $q$-dimensional standard normal random variable. The discrete form of the SDE also defines the proposal mechanism of the BM:
\begin{align}
\label{eqn:disProp}
q\left(x(t)|x(t-1)\right) = \mathbb N \left( x(t) |  \mu \left( x\left(t-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right ).
\end{align}

\begin{theorem} \label{theorem:sde}
The stochastic process defined in equation \ref{eqnBM} is coordinate independent. With a given $\partial t$, simulating any choice of local coordinates as above is equivalent as a step in $\mathbb{M}$.
\end{theorem}

The proof of Theorem \ref{theorem:sde} is given in appendix \ref{ax:independent}. This implies that the BM sample paths simulated from different choices of metric $\mathcal{G}$ should have the same properties. The proposal in equation \ref{eqn:disProp} guides BM steps according to the metric tensor. The boundary $\partial \mathbb{M}$ of the manifold is also quantified by the uncertainty of the mapping. We apply the Neumann boundary condition as in section \ref{sec:GPmani}. It implies the simulated sample paths only exist within the boundary. ???  An example of BM trajectories on Swiss roll is given in Fig. \ref{fig:BM}. The latent space and the associated metric is learned from Bayesian GPLVM. The grey colour in the latent space indicate the square root of the determinant of the metric $\sqrt{det(\mathcal{G}) }$. In machine learning community, $\sqrt{det(\mathcal{G}) }$ is also called the magnification factor. The geometry interpretation of the magnification factor is how much a small piece of the latent space in $\mathbb{R}^q$ will be stretched or compressed when it is mapped to the manifold embedded in $\mathbb{R}^p$. The horizontal axis of the latent space can be interpreted as the scaled radius of the Swiss roll, vertical axis is the width of the Swiss roll. When the radius is bigger, the corresponding area is darker and the magnification factor is larger. The trajectory of the stochastic process in the latent space is denser in the darker area and more spread in the bright area. This implies the stochastic process move with smaller step when the magnification factor is large and travel with bigger step when  the magnification factor is small. When we map the latent space back to the manifold in $\mathbb{R}^3$, the darker region will be stretched as the area of the outer layer of the Swill roll is bigger. The BM trajectory in the Swill roll still move with similar step size.

\subsection{Estimate the heat kernel with the transition density of Brownian motion} \label{BMinGP}
Considering the BM $\{ \mathbf{S}(t) | t>0 \}$ on $\mathbb{M} \subset \mathbb{R}^p$. The BM start from $\mathbf{S}(0)= s_0$ at time $0$. We simulate $N_{BM}$ sample paths. Given a point $s \in \mathbb{M}$, we define a small neighbourhood of $s$ as $\mathbf{A}_s \in \mathbb{M}$. For any $t > 0$, the probability of $\mathbf{S}(t)$ reach $\mathbf{A}_s$ at time t, $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$, can be approximated by counting how many paths reach $\mathbf{A}_s$.
\begin{align}
p( \mathbf{S}(t)\in \mathbf{A}_s\,| \ \mathbf{S}(0)=s_0)  \approx  \frac{N_{A_s}}{N_{BM}}
\end{align}
where $N_{A_s}$ is the number of sample paths reaching $\mathbf{A}_s$ at time $t$. This BM transition probability is defined as the integral of the BM transition density over $\mathbf{A}_s$ which is also the heat kernel $K_{heat}^t$.
\begin{align}
\label{BMprob}
p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0) = \int_{A_s} K_{heat}^t(s_0,s) ds
\end{align}
The integral is defined with respect to the volume form of $\mathbb{M}$. Since we do not have the analytical expression of the transition probability,  the transition density cannot be derived by taking the derivative of $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$. Instead, $K_{heat}^t$ can be numerically approximated as $ \hat{K}_{heat}^t $:
\begin{align}
\label{k:heat}
K_{heat}^t(s_0,s) \approx  \frac{p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0 )}{V(\mathbf{A}_s)} \approx \frac{1}{V(\mathbf{A}_s)}\cdot \frac{N_{A_s}}{N_{BM}} = \hat{K}_{heat}^t
\end{align}
where  $V(\mathbf{A}_s)$ is the Riemannian volume of $\mathbf{A}_s$, which is parameterised with the radius of $\mathbf{A}_s$. The Neumann boundary condition corresponds to BM reflecting at the boundary. This can be approximated by pausing time and resampling the next step until it falls into the interior of $\mathbb{M}$. This estimator is asymptotically unbiased and consistent \citep{niu2019}.  Note that $t$ is the BM diffusion time.  If $t$ is large, the BM paths have higher probability to reach the neighbourhood of the target point and leads to higher covariance and vice versa. The transition density can be estimated using algorithm \ref{alg:bm}.

\begin{corollary} \label{corollary:metric}
The transition density estimation in equation \ref{k:heat} is coordinate independent. The transition density estimation is equal for different choices of coordinate system.
\end{corollary}
By Theorem \ref{theorem:sde}, it is straight forward to have corollary \ref{corollary:metric}. The transition density estimation in equation \ref{k:heat} of the stochastic process in equation \ref{eqnBM} is also coordinate independent. Based on corollary \ref{corollary:metric}, we can compare the heat kernel estimates from different parameterisation of the same manifold. Here we take the Swiss roll as an example. Numerical accuracy of estimates of heat kernel (or BM transition density) on Swiss roll are shown in Fig. \ref{fig:metriccompare}. The heat kernel is estimated using equation \ref{k:heat} at twenty nine points with equally spaced radius $\in (3,9)$ of Swiss roll (see appendix ?? for details). The diffusion time is fixed as 50. The BM trajectories are simulated by using three different metric tensors. The red solid line in Fig. \ref{fig:metriccompare} represents the heat kernel estimated based on the analytical metric, assuming the geometry of the Swiss roll is known. The derivation of the analytical metric is shown appendix ??. When the geometry of the Swiss roll is unknown, the green dash line represents the heat kernel estimated based on the metric learned by Bayesian GPLVM and the blue dot line for GPLVM. The numerical comparison of GPLVM and Bayesian GPLVM is also shown inTable \ref{tb:metriccompare}. It is clear that the Bayesian GPLVM approach achieved better performance than GPLVM and give similar estimates as the analytical metric.

\begin{figure}[htb]
 \centering
\includegraphics[width=0.5\textwidth,height=0.4\textwidth]{metricCompare.pdf}
    \caption{ \label{fig:metriccompare}
       {\footnotesize
     Comparison of estimated BM transition density using analytical metric, metric from GPLVM and metric from Bayesian GPLVM.  The red solid line represents the heat kernel estimated based on the analytical metric. The green dash lines represent the heat kernel estimated with the Bayesian GPLVM metric. The blue dot line represent the heat kernel estimated with the metric learned from GPLVM. The black dot lines are the red curve plus and minus two times standard deviation of the estimated transition density with the analytical metric. The green line locates between the two black dotted lines.
     } }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Comparison of estimates of heat kernel using metric from GPLVM and B-GPLVM.  The table shows the root mean squared differences between the BM transition density estimated using analytical metric and the BM transition density using GPLVM metric and Bayesian GPLVM. }
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em Metric learning &\em Root mean squared differences\\
\hline
GPLVM & 1.2e-3\\
B-GPLVM & 7.4e-5 \\
\end{tabular}}
\end{table}

 \begin{algorithm}[!h]
 Learn the metric $\mathcal{G}$ from the point cloud $\mathcal{S} = \{  s_i | i = 1,\ldots,n+v\ \} $.        \Comment{ \footnotesize use eqn \ref{eqn:exmetric}  }
\caption{ { Simulating BM sample paths for estimating $K_{heat}^t$ } }
\label{alg:bm}
\begin{algorithmic}
\State 1.1 Generate Brownian motion trajectories
\For{$i = 1,\ldots,n$}    \Comment{ {\small $n$ is the size of data points}  }
\For{$j = 1,\ldots,N_{BM}$}    \Comment{ {\small$N_{BM}$ is No. of trajectories} }
\For{$ l = 1,\ldots,T$}    \Comment{ {\small T steps Brownian motion, $T*\Delta t$  $\rightarrow$ max diffusion time} }
\State {\bf do}    \Comment{ {\small keep proposing x until it falls in the boundary } }
\State $q\left(x_{i,j}(l)|x_{i,j}(l-1)\right) \leftarrow \mathbb N \left( x_{i,j}(l) |  \mu \left( x_{i,j}\left(l-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right )$ \Comment{ {\small use eqn \ref{eqn:disProp} } }
\State {\bf While} $x_{i,j}$ is not in the boundary $\partial \mathbb{M}$, $Var( \phi(x_{i,j})  ) > \alpha$.   \Comment{ \footnotesize use equation \ref{eqn:varb}  }
\EndFor
\EndFor
\EndFor
\Return $\bm x$
\State 1.2 Given a discrete choice of the diffusion time $t\in \{ \Delta t, 2*\Delta t,\cdot \cdot \cdot, T*\Delta t \}$, the covariance matrix $\Sigma^{t}$ is estimated based on the BM simulation from 1.1.
\For{$i = 1,\ldots,n$}
\For{$j = 1,\ldots,n$}
\State $N_{\mathbf{A}_j}$ = which( x(t) $\in$ $\mathbf{A}_j$ )  \Comment{ {\small counting how many BM paths reach $\mathbf{A}_j$}  }
\State $K_{heat}^t(s_i,s_j)= \frac{ N_{\mathbf{A}_j}}{N_{BM}*V( \mathbf{A}_j)}$  \Comment{ {\small use eqn \ref{k:heat} } }
\State $\Sigma_{ij}^t = \sigma_h^2 K_{heat}^t(s_i,s_j)$
\EndFor
\EndFor
\Return $\Sigma^t$
\end{algorithmic}
\end{algorithm}

\subsection{Optimising the kernel hyper parameters and make predication} \label{sec:optimpar}
Given a diffusion time $t$, we can generate the covariance matrix $\Sigma_{ {\text {\bf ff}} }^{t}$ for the training data using algorithm \ref{alg:bm}. $\Sigma_{ \text{\bf ff} }^t$ can be obtained as follows:  for the $i{th}$ row of $\Sigma_{ \text{\bf ff} }^t$, $N_{BM}$ trajectories are simulated, with the starting point the $i{th}$ data point corresponding to the index of row. For each element of $\Sigma_{ \text{\bf ff} }$, $\hat{K}^t(s_i,s_j)$ is then estimated using \eqref{k:heat}. The hyperparameters can be obtained by maximising the log of the marginal likelihood. The maximum of the BM diffusion time is set as $N_t \times \Delta t$, and $\Delta t$ is the BM simulation time step as defined in \eqref{eqn:disProp}.  $N_t$ covariance matrices $\Sigma_{ {\text {\bf ff}} }^{1 \ldots N_t}$ can be generated based on the BM simulations. Optimisation of diffusion time $t$ can be done by selecting the corresponding $\Sigma_{ {\text {\bf ff}} }^{t}$ that maximises the log marginal likelihood. Estimation of $\sigma_h$ follows the standard optimisation routines, such as quasi-Newton. The log marginal likelihood function (over $f$) is given by \citep{Rasmussen2004}:
\begin{align}
\label{loglike}
p( \bm{y}| s) =  \int p( \bm{y}| {\text {\bf f} }) p({\text {\bf f} }|s) d{\text {\bf f} }   = -\frac{1}{2} \bm y^T (\Sigma_{ {\text {\bf ff}} }^{t} + \sigma_{noise}^2 I)^{-1} \bm y - \frac{1}{2} \log|\Sigma_{ {\text {\bf ff}} }^t +\sigma_{noise}^2I | - \frac{n}{2}\log2\pi.
\end{align}

\section{Simulation studies on Swiss roll} \label{simdata}
 In this section, we carry out simulation studies for a regression model with true regression functions defined on a 2-dimensional Swiss Roll embedded in $\mathbb R^3$. Examples of BM sample paths on the Swiss Roll is shown in Figure \ref{fig:BM}.

\begin{figure}[htp]
    \centering
   \subfigure[ Swiss roll in $\mathbb{R}^3$ ]{\label{u:path}  \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swiss1.png} }
       \subfigure[ Latent space and variance of mapping]{ \label{fig:swissBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swissVar.png} }
    \subfigure[Latent space and metric]{ \label{fig:torusBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swissMag.png} }
    \caption{ \label{fig:BMexamples}
    \footnotesize{
The latent space constructed from B-GPLVM is plotted with the variance of the mapping as the background in (b), it is also plotted with the magnification factor as the background in (c).
      }
  }
\end{figure}

 \begin{figure}[htp]
    \centering
    \subfigure[True function and data points]{\label{u:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPtrue.pdf}}
     \subfigure[ \small{GP in $\mathbb{R}^3$ prediction }]{\label{u16:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGP3.pdf}}
    \subfigure[In-GP on manifold prediction]{\label{u:gp}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissmanifoldGP.pdf}}
    \subfigure[GP in $\mathbb{R}^2$ prediction]{ \label{u16:gp} \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGP2-2.pdf}}
    \caption{\label{fig:swissGP}
   { \footnotesize
     Comparison of in-GP  and Euclidean GP in the Swiss roll example. The true function and data points are shown in (a). The predictions of the GP which use the euclidean distance of data points in $\mathbb{R}^3$ are in (b). The in-GP on manifold prediction are plotted in (c). If we ignore the metric tensor and boundary, and construct a GP which use the $\mathbb{R}^2$ Euclidean distance in the latent space, the results are shown in (d).
    }
 }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Root mean squared error of regression on Swiss Roll }
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em &  \em $\mathbb{R}^3 GP$ &\em $\mathbb{R}^2 GP$ &\em  In-GP on manifold  \\
\hline
RMSE & 0.471  & 0.281 & 0.158 \\
\end{tabular}}
\end{table}

 The performance of in-GP is compared to that of an Euclidean  GP (the standard GP as in \cite{Rasmussen2004}) in $\mathbb{R}^3$.

\section{Application to image data} \label{simdata}
We consider images from the COIL data set ?? which consist of images from a fixed camera depicting objects on a motorised turntable against a black background. Each image is acquired after a 5 degree rotation. Here we consider the 66 images of object 14 (a toy cat). The camera angles are from 15 degree to 340 degree. Each image is treated as a 1024 dimensional vector.

\begin{figure}[H]
    \centering
   \subfigure[ COIL images ]{\label{fig:coilimage}  \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilimage.png} }
       \subfigure[ Latent space and variance of mapping]{ \label{fig:coilvar} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilLatentVar.png} }
    \subfigure[Latent space and metric]{ \label{fig:torusBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilLatentMag.png} }
    \caption{ \label{fig:coilMag}
    \footnotesize{
The latent space constructed from B-GPLVM is plotted with the variance of the mapping as the background in (b). The latent space is also plotted with the magnification factor as the background in (c).
      }
  }
\end{figure}

\begin{figure}[htb]
 \centering
\includegraphics[width=0.5\textwidth,height=0.4\textwidth]{coilBM.png}
    \caption{ \label{fig:coilBM}
       {\footnotesize
BM on the implicit manifold of COIL data.
     } }
\end{figure}

 \begin{figure}[H]
    \centering
    \subfigure[True function and data points]{\label{coil:true}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilTrue.png}}
     \subfigure[ \small{GP in $\mathbb{R}^{1024}$ prediction }]{\label{coil:Hd}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilHd.png}}
    \subfigure[In-GP on manifold prediction]{\label{coil:in-GP}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilIn-GP.png}}
    \subfigure[GP in $\mathbb{R}^2$ prediction]{ \label{coil:hd} \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coil2D.png}}
    \caption{\label{fig:coilGP}
   { \footnotesize
     Comparison of in-GP  and Euclidean GP in COIL images example. The true function is plotted with colour in the latent space in (a). The predictions of the GP which use the euclidean distance of data points in $\mathbb{R}^{1024}$ are in (b). The in-GP on manifold prediction are plotted in (c). If we ignore the metric tensor and boundary, and construct a GP which use the $\mathbb{R}^2$ Euclidean distance in the latent space, the results are shown in (d).
    }
 }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Root mean squared error of regression on COIL images}
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em &  \em $\mathbb{R}^{1024} GP$ &\em $\mathbb{R}^2 GP$ &\em  In-GP on manifold  \\
\hline
RMSE & 0.098  & 0.116 & 0.043 \\
\end{tabular}}
\end{table}

\section{Discussion} \label{discuss}
Our work proposes a novel class of intrinsic Gaussian processes on manifolds and complex constrained domains employing the equivalence relationship between heat kernels and the transition density of Brownian motions on manifolds. One of the key features of the in-GP is to fully incorporate the intrinsic geometry of the spaces for inference while respecting the potentially complex boundary or interior constraints. To reduce the computational cost of simulating BM sample paths when the sample size is large, sparse in-GPs are developed leveraging ideas from the literature on fast computation in GPs in Euclidean spaces. The results in section \ref{simdata} and \ref{real} indicate that in-GP achieves significant improvement over usual GPs. Although we did not conduct a formal asymptotic study for in-GPs, with insights gained from the Euclidean GP with squared exponential kernel in the Euclidean space (see, e.g., \cite{vandervaart}), we expect in-GPs to yield posterior consistency with respect to appropriate neighbourhoods of the true regression function $f_0$. The focus of this article has been on developing in-GPs on manifolds with known metric tensors.  There has been abundant interest in learning of unknown lower-dimensional  manifold structure in high-dimensional data.  in-GPs can be combined with these approaches for performing supervised learning on lower-dimensional latent manifolds.

\section{Introduction}\label{intro}
There is increasing interest in the problem of nonparametric regression with high-dimensional predictors. In a variety of fields, from computer science to environmental science, one often encounters high-dimensional data (e.g., ‘point cloud data’) perturbed by high-dimensional noise but cantering around some lower-dimensional manifolds. In more precise mathematical terms, manifolds are topological spaces equipped with some differential/smooth structure, the geometry of which is in general different from the usual Euclidean geometry. Naively applying traditional multivariate analysis to manifold-valued data that ignores the geometry of the space can potentially lead to highly misleading predictions and inferences. Gaussian Processes (GP) are among the powerful tools in statistics and machine learning for regression with Euclidean predictors\citep{Rasmussen2004}. However, if the predictors locate on a manifold, traditional smoothing or modelling methods like GP that do not respect the intrinsic geometry of the space and the boundary constraints, produce poor results \citep{extrinsicGP, niu2019}. In this work, we focus on the following model throughout the paper,
\begin{align} \label{eqn:reg}
y_i = f(s_i) +\epsilon_i, \ \ \  \epsilon_i \sim \mathcal{N}(0, \sigma_{noise}^2)
\end{align}
where  $f:\mathbb{M} \rightarrow \mathbb{R}$ is an unknown regression function. $y_i \in \mathbb{R}$ is an observed response variable. $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is an observed predictor on the embedded manifold in $\mathbb{R}^p$, $s_i = (s_{i}^1,\ldots,s_{i}^p)$. $\mathbb{M}$ is a complete and orientable Riemannian manifold embedded in $\mathbb{R}^p$. We would like to do inferences on how the output $y$ varies with the input $s$, including predicting $y$ values at new locations not represented in the training dataset. A GP prior can be assigned to $f$ with a covariance function $k(s_i,s_j)$. The choice of $k(s_i,s_j)$ has a fundamental impact on the results. Most common choices of covariance function such as the squared exponential kernel and Matern kernel https://www.overleaf.com/project/5ff041180b93f20680f0d6ecdepend critically on the Euclidean distance between $s_i$ and $s_j$ and ignore the intrinsic geometry of $\mathbb M$.

One of the paramount challenges in developing GP models on manifold is the difficulty in specifying the covariance structure via constructing valid and computable covariance kernels on manifolds.  \cite{niu2019} proposed to use the heat kernel to construct the intrinsic Gaussian process (In-GP) on complex constrained domain. Although the heat kernel provides a natural and canonical choice theoretically, it is analytically intractable to directly evaluate \citep{hsu1988}. \cite{Castillo2014} instead proposed to use randomly rescaled solutions of the heat equation to define a valid covariance kernel for reasonably broad classes of compact manifolds. But their proposed heat kernels are computationally intractable. Alternatively, it can be estimated as the transition density of the Brownian Motion (BM) on the manifold \citep{hsu2008}.  By using the heat kernel, In-GP can accommodate spatial domains arising as complex subsets of Euclidean space. \cite{niu2019} has demonstrated the application of In-GP in modelling the chlorophyll concentration levels in Aral Sea \citep{wood} and smooth functions on Swiss roll.  It respects the potentially complex boundary or interior conditions as well as the intrinsic geometry of the spaces. However \cite{niu2019}'s method is only applicable when the geometry of the manifold is known and the dimension of the dataset space is small such as the subspace of $\mathbb{R}^2$ or $\mathbb{R}^3$.  In this paper, we develop a computational inferential framework able to estimate the heat kernel and construct the In-GP on a unknown manifold embedded in the high dimensional point cloud.

\begin{figure}[htb]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{map.pdf}
    \caption{ \label{fig:mapping}
    \footnotesize{
    { An illustration of the parameterisation of the manifold: $\phi$ maps the chart into the embedding space }}
    }
\end{figure}
In statistics and machine learning community, we study probabilistic latent variable models (LVMs) as the parameterisation or embeddings of the manifold in the observation space. The deterministic LVM is first studied to describe curved surface By \cite{gauss1828}. He shown a low dimensional surface embedded in a high dimensional Euclidean space can be well described through a mapping function.  The low dimensional representation of the surface is known as the chart which corresponds to the latent space in probabilistic LVMs. Intuitively, the chart provides a distorted view of the surface. An illustration is shown in Fig. \ref{fig:mapping}. $x^1$ and $x^2$ coordinates in the chart represent the radius and width of the Swiss roll in $\mathbb{R}^3$.  Points in the chart can be mapped to the embedded space (Swiss roll)  through the function $\phi$. In other words, the Swiss roll can be parameterised by the chart. Measurements on the manifold can be computed in the chart locally, and integrated to provide global measures. This gives rise to the definition of a local inner product, known as a Riemannian metric tensor. A Riemannian metric $\mathbf{g}$ on a manifold $\mathbb{M}$ is symmetric and positive definite matrix which defines a smoothly varying inner product in the tangent space.
\begin{align} \label{eqn:deter_metric}
\mathbf{g} = \mathcal{J}^T \mathcal{J}, \ \ \  \mathcal{J}_{i,j} = \frac{\partial \phi_i }{ \partial  x^j}
\end{align}
where $\mathcal J$ denotes the Jacobian of the $\phi$. In \cite{niu2019} , both the mapping function $\phi$ and the metric tensor $\mathbf{g}$ are known.

In modern big data analysis, the data in the point cloud, often high dimensional, are not directly observed on the manifolds, instead they are observed in a potentially high-dimensional ambient space but concentrate around some unknown lower dimensional structure. The mapping function $\phi$ is unknown. One needs to learn this geometry before utilizing it for inference. Most recent research in \cite{dunson2020diffusion} tackled this problem by approximating the heat kernel with finitely-many eigenpairs of the Graph Laplacian. Since they are using a compact operator to approximate an unbounded operator, hence the eigenvalues are hard to recover and need many points in the point cloud. \cite{dunson2020diffusion}  has shown 3000 points are required in order to approximate the heat kernel of a simple one-dimensional unit circle or spiral from the point cloud in $\mathbb{R}^2$. In addition, they  did not consider the boundary condition. In many applications such as image processing, the dimension of the data in the point cloud  easily scales up to thousands. Neither  \cite{niu2019}   or \cite{dunson2020diffusion} work for this type of problems. There are still some critical gaps in current practice in particular in lacking of robust methods for carrying out intrinsic statistical inference, and largely missing effective models for regressions with manifold-valued data embedding in a point cloud.  In this paper, we will use the Bayesian dimension reduction method to learn the metric tensor of the underlying manifold of the point cloud and define its boundaries. The heat kernel can be estimated as the transition density of the BM on the learned manifold. We have proved the BM transition density estimation is coordinate independent. Different metric learning methods such as GPLVM (Gaussian Processes Latent Variable Model) and Bayesian GPLVM are compared. This framework allows us to build the In-GP models on the implicit manifold embedded in the point cloud. It can incorporate fully the intrinsic geometry of the implicit manifold for inference while respecting the potentially complex interior structure and boundary.

\section{Intrinsic Gaussian process on point cloud}
\label{sec:GPmani}
We propose to construct the intrinsic Gaussian processes on implicit manifold by using the heat kernel as the covariance kernel. Let $\mathbb{M}$ be a $q$-dimensional complete and orientable Riemannian manifold with the Riemannian metric $\mathbf{g}$, $\partial \mathbb{M}$ its boundary that is continuous. A heat kernel of $\mathbb{M}$ is a smooth function $K(s_0,s,t)$ on $ \mathbb{M} \times \mathbb{M} \times \mathbb{R}^{+}$ that satisfies the heat equation:
\begin{align*}
\frac{\partial}{\partial t}K_{heat}(s_0,s,t)=\frac{1}{2}\Delta_s K_{heat}(s_0,s,t),\qquad  s_0, s\in \mathbb{M},
\end{align*}
where $\Delta_s$ is the Laplacian-Beltrami operator on $\mathbb{M}$, and $t>0$ is the diffusion time. The heat kernel satisfies the initial condition $\lim_{t\rightarrow 0}K_{heat}(s_0,s,0)=\delta(s_0,s)$ with $\delta$ the Dirac delta function. The heat kernel becomes unique when we impose a suitable condition along the boundary $\partial \mathbb{M}$, such as the Neumann boundary condition. $ \partial K / \partial \mathbf{n} = 0$ along $\partial \mathbb{M}$ where $\mathbf{n}$ denotes a normal vector of $\partial \mathbb{M}$.  The Neumann boundary condition can be expressed as no heat transfer across the boundary. Here we use $K_{heat}^t$ to denote the heat kernel of arbitrary manifold. In the local coordinate (chart), the Laplacian-Beltrami operator can be written as
 \begin{align}
\label{eqn:delta_s}
\Delta_s K_{heat}^t = \frac{1}{\sqrt{G}} \frac{\partial}{\partial x^j} \left (    \sqrt{G} \mathbf{g}^{ij} \frac{\partial K_{heat}^t}{\partial x^i}  \right ),
\end{align}
where $G$ is the determinant of metric $\mathbf{g}$ and $\mathbf{g}^{ij}$ is the $(i,j)$ element of its inverse. Take the special case $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$, $\mathbf{g}$ become an identity matrix. The Laplace Beltrami operator become the Laplace operator.    The heat kernel has a closed form corresponding to a time-varying Gaussian function:
\begin{align*}
K_{heat}(\mathbf{s}_0,\mathbf{s},t)
=\frac{1}{(2\pi t)^{q/2}}\,
  \exp\left\{-\frac{||\mathbf{s}_0-\mathbf{s}||^2}{ 2t }\right\}, \; \mathbf{s}\in \mathbb R^q.
\end{align*}
The diffussion time $t$ controlling the rate of decay of the covariance. By varying the time parameter, one can vary the bumpiness of the realisations of the in-GP.

We now consider the heat kernel as the covariance function of the In-GP. Let $\mathcal{D} =  \{ (s_i, y_i), i = 1, . . . , n \} $, be the data, with $n \geq 1$ the number of labeled observations. $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is the predictor and $y_i \in \mathbb{R}$ is the corresponding response variable. Suppose we are also given an unlabeled dataset $\mathcal{V} = \{ s_i, i=n+1 . . . ,n+v \} $ where $v \geq 1$. Consider the regression model in \ref{eqn:reg}, we would like to do inferences on $f$ with the labeled dataset $\mathcal{D}$ and predict $y$ values at the unlabeled dataset $\mathcal{V} $. Under an in-GP prior for the unknown regression function as $f \sim GP(0,K_{heat}^t(.,.) )$, we have
\begin{align*}
p( \mathbf{f} | s_1, s_2, ..., s_n) \sim \mathcal{N} (0, \Sigma_{\mathbf{f} \mathbf{f}}),
\end{align*}
where $\mathbf{f}$ is the discretization of $f$ over $s_1,s_2,\cdot\cdot\cdot,s_n$ so that $f_i = f(s_i)$. $\Sigma_{\mathbf{f} \mathbf{f}} \in \mathbb{R}^{n\times n}$ is the covariance matrix induced from the heat kernel. The $(i,j)$ entries of $\Sigma_{\mathbf{f}\f}$ corresponds to $\Sigma_{\mathbf{f}\f _{i,j}} = \sigma_h^2 K_{heat}^t(s_i,s_j)$. We introduce the rescaling hyperparameter $\sigma_h^2$ to add extra flexibility of the heat kernel. Let ${\bf f}_*$ be a vector of values of $f(\cdot)$ at the unlabelled grid points $\mathcal{V}$. Under the regression model in equation \ref{eqn:reg} and the GP prior, we have the joint distribution of ${\bf y} $ and ${\bf f}_*$ is:
\begin{align}
p( {\text {\bf y}}, {\text {\bf f}}_*) = \mathcal{N} \left ( 0,
 \left [ \begin{array}{cc}
 \Sigma_{{\text {\bf f}} {\text {\bf f}}} + \sigma_{noise}^2 I_n &\Sigma_{ {\text {\bf f}}  {\text {\bf f}}_*} \\[0.3em]
 \Sigma_{ {\text {\bf f}}_* {\text {\bf f}}} &\Sigma_{ {\text {\bf f}}_* {\text {\bf f}}_*} \\[0.3em]
                \end{array} \right ]
                \right ),
\end{align}
where $\Sigma_{{\text {\bf f}}_*{\text {\bf f}}}$ is the covariance matrix for training and unlabelled grid data points. The predictive distribution is derived by marginalising out {\bf f}:
\begin{align}
p( {\text {\bf f}}_*| \bm y) = \int p( {\text {\bf f}}_*{\text {\bf f}}| \bm y) d{\text {\bf f}} = \mathcal N \left ( \Sigma_{ {\text {\bf f}}_*{\text {\bf f}}}  \left(\Sigma_{{\text {\bf f}}{\text {\bf f}}}  + \sigma_{noise}^2 I \right)^{-1}\bm y , \ \  \Sigma_{{\text {\bf f}}_*{\text {\bf f}}_*}-\left(\Sigma_{{\text {\bf f}}{\text {\bf f}}} +\sigma_{noise}^2 I \right)^{-1} \Sigma_{{\text {\bf f}}{\text {\bf f}}_*} \right)
\end{align}

The main difficulty of using In-GP is the analytic intractability of evaluating $K_{heat}^t$. The analytical solution of heat kernel only exists for very few special cases such as Euclidean space. Solving the heat equation for an arbitrary manifold is nontrivial. \cite{niu2019} estimate the heat kernel as the BM transition density by simulating BM paths on a known manifold. In the case of point cloud, both the metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$ are unknown. In this work, we will use the probabilistic latent variable models to learn the metric $\mathbf{g}$  and defining the boundary $\partial \mathbb{M}$ based on the uncertainty quantification of the mapping. The heat kernel can then be estimated by simulating BM paths on the learned manifold.

\section{Learning metric and boundary}
Let $\mathcal S = \{ s_i |  i=1,\cdot \cdot \cdot, n+v \}$,  $s_i \in \mathbb{R}^p$, include both predictors in the labeled dataset $\mathcal{D}$ and unlabelled dataset $\mathcal{V}$. We have $\mathcal{S} \in \mathbb{R}^{ (n+v) \times p }$. Suppose we perform probabilistic nonlinear dimensionality reduction by defining a latent variable model that introduces a set of latent (unobserved) variables $\mathcal{X}= \{x_i | i= 1,\cdot \cdot \cdot, n+v \}$, $x_i \in \mathbb{R}^q$ and $q < p$.  $\mathcal{X} \in \mathbb{R}^{(n+v)\times q}$. $\mathcal X$ is related to $\mathcal S$ which is observed in a higher dimensional space. The latent space is dominated by a prior distribution $p( \mathcal X )$ which induces a distribution over $\mathcal S$ under the assumption of a probabilistic mapping
\begin{align} \label{eqn:mapmodel}
s_{i}^j = \phi^j(x_i) + e_{i}^j,
\end{align}
where $x_i \in \mathbb{R}^q$ is the $i^{th}$ latent point associated with $s_i \in \mathbb{R}^p$ in the original space, $i<n$. $j$ is the index of the features (dimensions) of $s$ in the observation space and $j<p$. $e_{i}^j$ is a noise term. The noise is Gaussian distributed $e_i^{j} \sim \mathcal N(0,\beta^2)$, $\beta^2$ is the variance. An illustration is given in Fig. \ref{fig:mapping} where $q=2$ and $p=3$. If the mapping $\phi$ is linear and the prior $p( \mathcal X)$ is Gaussian, this model is known as probabilistic principal component analysis \citep{tipping1998}. In this work, we do not restrict to this linear assumption and consider some nonlinear dimensionality reduction methods such as Gaussian Process Latent Variable Model (GP-LVM) and Bayesian Gaussian Process Latent Variable Model (BGP-LVM).

When the mapping $\phi$ in Fig.\ref{fig:mapping} is differentiable, it can be interpreted as the mapping between the chart (or latent space) $\mathbb{R}^q$ and the embedding space $\mathbb{M} \subset \mathbb{R}^p$.  The natural Riemannian metric of the given model can be computed as in equation \ref{eqn:deter_metric}. In the case of probabilistic LVMs where the conditional probability over the Jacobian follow a Gaussian distribution, this naturally induces a distribution over the local metric tensor $\mathbf{g}$. We denote the distribution of the Jacobian as $\mathbf{J}$. Assuming independent rows of $\mathbf J$
\begin{align}
p(\mathbf J | \mathcal{S}, \mathbf{f}) = \prod_{j=1}^p \mathcal{N}( \mu_{\mathbf{J}}^j, \Sigma_{\mathbf J}  )
\end{align}
The expression of the mean $\mu_{\mathbf J}$ and variance $\Sigma_{\mathbf J} $ of the Jacobian are model specific and given in section \ref{gplvm} and \ref{bgplvm}. This resulting the metric $\mathbf{g}$ follow a non-central Wishart distribution \citep{anderson1946}
\begin{align} \label{eqn:metric}
\mathbf{g} \sim \mathcal{W}_q \left(  p, \Sigma_{\mathbf{J}}, \E(\mathbf{J}^T)  \E(\mathbf{J})  \right).
\end{align}
From this distribution, the expected metric tensor can be computed as
\begin{align} \label{eqn:exmetric}
\E( \mathbf{g} ) =  \mathcal{G} = E( \mathbf{J}^T) \E( \mathbf{J} ) + p \Sigma_{\mathbf{J}}.
\end{align}
We denote the expectation by $\mathcal{G}$. We also need the gradient of the expected metric to simulate BM as in section \ref{trandens}.
\begin{align}
\frac{ \partial \E[ \mathbf{g}]}{\partial x^l} = \frac{\partial \mathcal{G}}{\partial x^l} = \frac{\partial \E[\mathbf{J}^T] } { \partial x^l }  \E[\mathbf{J}] + \E[\mathbf{J}^T] \frac{ \partial  \E[\mathbf{J}]  }{\partial x^l}  + p \frac{ \partial \Sigma_\mathbf{J} }{ \partial x^l }
\end{align}
The boundary of the learned manifold can be defined by checking $Var( \phi(x) | x)$, the variance of the mapping at $x$. Any $x$ outside of the boundary has $Var( \phi(x) | x ) > \alpha$, $\alpha$ is a predefined threshold. Both $\mathcal{G}$ and $\partial \mathbb{M}$  are used to construct BM on manifold in section \ref{}.   ???
\begin{align}\label{eqn:bound}
\partial \mathbb{M} = \{  x \in \mathbb{R}^q | Var( \phi(x) | x) = \alpha   \}.
\end{align}

\subsection{GPLVM}\label{gplvm}
Gaussian Process Latent Variables Model \citep{lawrence2005} is a nonlinear dimensionality reduction method. It defines a generative mapping from the latent space $x\in \mathbb{R}^q$ to observation space $s \in \mathbb{M} \subset R^p$ that is governed by Gaussian processes. GPs priors are placed on $\phi$ in equation \ref{eqn:mapmodel}. The GPs are taken to be independent across the dimensions of the observation space. This leads to the formulation of GP-LVM. Given the construction outlined above, the joint probability over the observed data $\mathcal{S}$ and all latent variables $\mathcal{X}$ is written as follows:
\begin{align}\label{eqn:gplvmJoint}
p(\mathcal{S} \ \Phi \ \mathcal{X} | \beta) =   p(\mathcal{S} | \Phi \beta ) p(\Phi | \mathcal{X}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j, \beta ) p( \bm{\phi}_{:}^j | \mathcal{X}) p(\mathcal{X})
\end{align}
 where $\Phi = \{ \bm{\phi}_i | i=1,\cdot\cdot\cdot, n+v \} $, $\Phi \in \mathbb{R}^{(n+v)\times p}$ and $ \bm{\phi}_i^j = \phi(x_i)^j$.  $\mathbf{s}_{:}^j$ represents the $j_{th}$ dimension of every points in $\mathcal S$ and $\mathcal S \in \mathbb{R}^{ (n+v)\times p}$. The likelihood $p( \mathcal S | \mathcal{X} )$ is computed by marginalising out $\Phi$ and optimising the latent variables $\mathcal{X}$.
\begin{align}
p(\mathcal{S} | \mathcal{X} ) = \prod_{j=1}^p \mathcal N( \mathbf{s}_{:}^j, K_{\mathcal{X} \mathcal{X}} + \beta^2 I )
\end{align}
$K_{\mathcal{X} \mathcal{X}} $ is the $(n+v)\times (n+v)$ covariance matrix defined by the squared exponential kernel. Since $x \in \mathcal{X}$ is a latent variable, we can assign it a prior density given by the standard normal density. \cite{lawrence2005} estimate all latent variables $\mathcal X$ and kernel hyper-parameters of GP-LVM by following the maximum likelihood estimation. If the covariance kernel is differentiable, the derivative of a GP is again a GP \citep{Rasmussen2004}. This property allows predictions about the derivative of GP. The Jacobian $\mathbf{J}$ of the GP-LVM mapping can be computed as the partial derivative $\frac{\partial \phi_*}{ \partial x^l}$ with respect to the $l_{th}$ dimension for any point $x_*$  in the latent space.
\begin{align} \label{eqn:Jacobian}
\mathbf{J}^{T} = \frac{\partial \phi_*}{ \partial x} = \begin{bmatrix}
 \frac{ \partial \phi(x_*)^1 }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^1}   \\
  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot \\
  \frac{ \partial \phi(x_*)^1 }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^q}
\end{bmatrix}
\end{align}
where $\frac{\partial \phi_*}{ \partial x}$ is a $q \times p$ matrix. Considering the independence across the dimensions of the observation space, the joint distribution of the $j_{th}$ dimension of the mapping $\phi$ and the $j_{th}$ column of the Jacobian can be written as
\begin{align}
\begin{bmatrix}
 \phi( \mathcal{X})^j\\
 \frac{ \partial \phi(x_*)^j }{\partial x}
\end{bmatrix}, \
\sim \mathcal N \left( 0, \begin{bmatrix}
K_{\mathcal{X},\mathcal{X}} & \partial K_{\mathcal{X},*} \\
 \partial K_{\mathcal{X},*}^T & \partial^2 K_{*,*}
\end{bmatrix}  \right )
\end{align}
The expressions of $K_{\mathcal{X},\mathcal{X}}$ ,  $\partial K_{\mathcal{X},*}$ , and  $\partial^2 K_{*,*} $ are given in the appendix ??. The GP-LVM model provides an explicit mapping from the latent space to the observed space. This mapping defines the support of the observed data $\mathcal{S}$ as a $q$ dimensional manifold embedded into $\mathbb{R}^p$. The distribution of the Jacobian of the GP-LVM mapping is the product of $p$ independent Gaussian distributions (one for each dimension of the observation space). For a latent point $x_*$, the  distribution of the Jacobian takes this form
\begin{align}
p(  \mathbf{J} | \mathcal{X}, \mathcal{S}) &= \prod^p_{j=1} \mathcal N( \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \partial K^T_{\mathcal{X},*} K^{-1}_{\mathcal{X},\mathcal{X}} \mathbf{s}_{:}^j  \ , \partial^2K_{*,*} - \partial K^T_{\mathcal{X},*} K_{\mathcal{X},\mathcal{X}}^{-1} \partial K_{\mathcal{X},*}  ) \nonumber
\end{align}
From this distribution, the expected metric tensor can be computed as in equation \ref{eqn:exmetric}. The boundary $\partial M$ is defined by $Var( \phi(x_*) | x_*)$. For any $x_*$ in the latent space,
\begin{align}\label{eqn:varb}
Var(\phi(x_*) | x_*) = K_{*,*} - K_{*,\mathcal{X}}K_{\mathcal{X},\mathcal{X}}^{-1} K_{\mathcal{X},*}.
\end{align}

\subsection{Bayesian-GPLVM }\label{bgplvm}
The GP-LVM based models rely on MAP training procedures for optimising the latent inputs and the hyperparameters \citep{lawrence2005}. However, the standard GP-LVM does not mariginalise out the latent inputs $\mathcal{X}$ implies that it could be sensitive to overfitting \citep{damianou2016}. The marginal likelihood takes this form
\begin{align*}
p(\mathcal S)=\int p(\mathcal{S} \mid \mathcal{X}) p(\mathcal{X}) d \mathcal{X}
\end{align*}
The  integral is intractable as the inputs $\mathcal{X}$ of the latent variable model appears nonlinearly in the inverse of the covariance matrix. Bayesian Gaussian process latent variable model (BGP-LVM, \cite{titsias2010}) introduce a variational inference framwork for training the latent variable model and perform the nonlinear dimensional reduction. It allows variationally integrate out the input variables and compute a lower bound on the exact marginal likelihood of the nonlinear latent variable model. The maximisation of the variational lower bound provides a Bayesian training procedure that is robust to overfitting.

The key to make the variational Bayes approach tractable is the application of variational inference in an expanded probability model where the GP prior on $\phi$ is augmented to include auxiliary variables. More specifically, we expand the joint probability model in equation \ref{eqn:joint} by including $m$ extra samples (inducing points) of the GP latent mapping $\phi(x)$, so that $u_i = \phi(x_{ui}) \in \mathbb{R}^p$ is such a sample \citep{lawrence2007}. These inducing points are denoted by $\mathcal{U} = \{u_i| i=1,\cdot\cdot\cdot,m   \}$, $\mathcal{U} \in \mathbb{R}^{m\times p}$ and constitute latent function evaluations at a set of pseudo-inputs $\mathcal{X}_u = \{ x_{ui} | i=1,\ldots,m \} $, $\mathcal{X}_u \in \mathbb{R}^{m \times q}$. The inducing inputs $\mathcal{X}_u$ are not random variables. $\Phi$ is defined as in equation \ref{eqn:gplvmJoint}. The augmented joint probability and the marginal likelihood can be written as
\begin{align}\label{eqn:full-marginal}
p(\mathcal{S} \ \Phi \ \mathcal{U} \  \mathcal{X} | \beta) &=   p(\mathcal{S} | \Phi ) p(\Phi | \mathcal{U} \mathcal{X}) p(\mathcal {U}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j) p( \bm{\phi}_{:}^j |  u_{:}^j, \mathcal{X}) p(u_{:}^j) p(\mathcal{X}) \\
p(\mathcal{S}) &= \int  \int  \int  p(\mathcal{S}, \Phi, \mathcal{U}, \mathcal{X}) \ d \mathcal{U}     \  d\Phi \ d \mathcal{X}
\end{align}
Using this variational inference, $log\left(p(\mathcal{S} )\right)$ can be expressed as the Jensen’s lower bound. We can analytically maximize the lower bound. The full technical details of the lower bound are given in the appendix ??.  In contrast to the standard GPLVM, Bayesian GPLVM variationally integrate the latent variables $\mathcal{X}$ out and obtain extra regularisation. The Jacobian have the same shape as in equation \ref{eqn:Jacobian}. For a latent point $x_*$ the distribution the distribution of the Jacobian for the Bayesian GP-LVM mapping takes the form:
\begin{align}
p(  \mathbf{J} | \mathcal{X},\mathcal{U},\mathcal{S}) &= \prod^p_{j=1} \mathcal N( \  \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \  \partial K^T_{\mathcal{X}_u,*} K^{-1}_{\mathcal{X}_u,\mathcal{X}_u} u_{:}^j , \  \partial^2K_{*,*} - \partial K^T_{\mathcal{X}_u,*} K_{\mathcal{X}_u,\mathcal{X}_u}^{-1} \partial K_{\mathcal{X}_u,*}  ) \nonumber
\end{align}
The expressions of $ \partial K_{\mathcal{X}_u,*} $ and $\partial^2K_{*,*} $ are given in the appendix??. From this distribution, the expected metric tensor can be computed as in equation \ref{eqn:exmetric}.

\section{Estimate heat kernel as BM transition density}
The construction of the heat kernel associated with the Laplace-Beltrami operator is not a trivial task and belongs to the fields of partial differential equations and differential geometry. It is carried out in great detail in \cite{chavel1984}. The Laplace Beltrami operator in equation \ref{eqn:delta_s} is also the infinitesimal generator of Brownian Motion on the manifold \citep{hsu2008}.  Let
$\mathbf{S}_t$ be a Brownian motion with transition density $K_{heat}^t$. This transition density satisfies the heat equation. In this section, we will estimate the heat kernel as the BM transition density to avoid solving the heat equation on manifold.

\subsection{ Simulating Brownian motion on manifolds }
\label{trandens}
From last section, we have learned the parameterisation of the unknown manifold $\mathbb{M}$ as the latent space $\mathcal{X}$, the distribution of the associated metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$. In order to estimate the BM transition density on $\mathbb{M}$, we first need to simulate BM trajectories. Simulating a sample path of BM on $\mathbb{M} \subset \mathbb{R}^p$ is equivalent to simulating a stochastic process in the latent space in $\mathbb{R}^q$, $q<p$. An illustration example is shown in Fig. \ref{fig:BM}.

\begin{figure}[htb]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{bmlatent1.pdf}
    \caption{ \label{fig:BM}
    \footnotesize{
    { BM on $\mathbb{M}$ and its equivalent stochastic process in chart in $\mathbb{R}^2$.  $\phi: \mathbb{R}^2 \rightarrow \mathbb{M}$ is a parametrisation of $\mathbb{M}$.  }}
    }
\end{figure}

The BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations in the Ito form \citep{hsu1988,hsu2008}. We will use the expectation $\mathcal{G} = \E[\mathbf{g}]$ to construct this SDE.
\begin{align}
\label{eqnBM}
dx^i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x^j} \left(  \mathcal{G}^{ij}G^{1/2} \right) dt + \left( \mathcal{G}^{-1/2} dB(t)\right)_i,
\end{align}
where $x^i$ represent the $i_th$ dimension of the latent space (chart). $\mathcal{G}$ is the metric tensor of $\mathbb{M}$, $G$ is the determinant of $\mathcal{G}$ and $B(t)$ represents an independent BM in the Euclidean space. In order to do forward simulation, the discrete form of equation (\ref{eqnBM}) is derived in equation \eqref{disBM}.
\begin{align}
\label{disBM}
x^i(t) &= x^i(t-1) + \frac{1}{2} \sum^{q}_{j=1} \left(  -\mathcal{G}^{-1} \frac{\partial \mathcal{G}}{\partial x^j} \mathcal{G}^{-1} \right)_{ij} \Delta t + \frac{1}{4} \sum_{j=1}^q(\mathcal{G}^{-1})_{ij}tr(\mathcal{G}^{-1}\frac{\partial \mathcal{G}}{\partial x^j}) \Delta t + \left( \mathcal{G}^{-1/2} dB(t)\right)_i \nonumber \\
         &= \mu( x^i(t-1),\Delta t)_i + \left( \sqrt{\Delta t} \mathcal{G}^{-1/2} z^q\right)_i,
\end{align}
where $\Delta t$ is the diffusion time of the BM simulation step and $z^q$ is a $q$-dimensional standard normal random variable. The discrete form of the SDE also defines the proposal mechanism of the BM:
\begin{align}
\label{eqn:disProp}
q\left(x(t)|x(t-1)\right) = \mathbb N \left( x(t) |  \mu \left( x\left(t-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right ).
\end{align}

\begin{theorem} \label{theorem:sde}
The stochastic process defined in equation \ref{eqnBM} is coordinate independent. With a given $\partial t$, simulating any choice of local coordinates as above is equivalent as a step in $\mathbb{M}$.
\end{theorem}

The proof of Theorem \ref{theorem:sde} is given in appendix \ref{ax:independent}. This implies that the BM sample paths simulated from different choices of metric $\mathcal{G}$ should have the same properties. The proposal in equation \ref{eqn:disProp} guides BM steps according to the metric tensor. The boundary $\partial \mathbb{M}$ of the manifold is also quantified by the uncertainty of the mapping. We apply the Neumann boundary condition as in section \ref{sec:GPmani}. It implies the simulated sample paths only exist within the boundary. ???  An example of BM trajectories on Swiss roll is given in Fig. \ref{fig:BM}. The latent space and the associated metric is learned from Bayesian GPLVM. The grey colour in the latent space indicate the square root of the determinant of the metric $\sqrt{det(\mathcal{G}) }$. In machine learning community, $\sqrt{det(\mathcal{G}) }$ is also called the magnification factor. The geometry interpretation of the magnification factor is how much a small piece of the latent space in $\mathbb{R}^q$ will be stretched or compressed when it is mapped to the manifold embedded in $\mathbb{R}^p$. The horizontal axis of the latent space can be interpreted as the scaled radius of the Swiss roll, vertical axis is the width of the Swiss roll. When the radius is bigger, the corresponding area is darker and the magnification factor is larger. The trajectory of the stochastic process in the latent space is denser in the darker area and more spread in the bright area. This implies the stochastic process move with smaller step when the magnification factor is large and travel with bigger step when  the magnification factor is small. When we map the latent space back to the manifold in $\mathbb{R}^3$, the darker region will be stretched as the area of the outer layer of the Swill roll is bigger. The BM trajectory in the Swill roll still move with similar step size.

\subsection{Estimate the heat kernel with the transition density of Brownian motion} \label{BMinGP}
Considering the BM $\{ \mathbf{S}(t) | t>0 \}$ on $\mathbb{M} \subset \mathbb{R}^p$. The BM start from $\mathbf{S}(0)= s_0$ at time $0$. We simulate $N_{BM}$ sample paths. Given a point $s \in \mathbb{M}$, we define a small neighbourhood of $s$ as $\mathbf{A}_s \in \mathbb{M}$. For any $t > 0$, the probability of $\mathbf{S}(t)$ reach $\mathbf{A}_s$ at time t, $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$, can be approximated by counting how many paths reach $\mathbf{A}_s$.
\begin{align}
p( \mathbf{S}(t)\in \mathbf{A}_s\,| \ \mathbf{S}(0)=s_0)  \approx  \frac{N_{A_s}}{N_{BM}}
\end{align}
where $N_{A_s}$ is the number of sample paths reaching $\mathbf{A}_s$ at time $t$. This BM transition probability is defined as the integral of the BM transition density over $\mathbf{A}_s$ which is also the heat kernel $K_{heat}^t$.
\begin{align}
\label{BMprob}
p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0) = \int_{A_s} K_{heat}^t(s_0,s) ds
\end{align}
The integral is defined with respect to the volume form of $\mathbb{M}$. Since we do not have the analytical expression of the transition probability,  the transition density cannot be derived by taking the derivative of $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$. Instead, $K_{heat}^t$ can be numerically approximated as $ \hat{K}_{heat}^t $:
\begin{align}
\label{k:heat}
K_{heat}^t(s_0,s) \approx  \frac{p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0 )}{V(\mathbf{A}_s)} \approx \frac{1}{V(\mathbf{A}_s)}\cdot \frac{N_{A_s}}{N_{BM}} = \hat{K}_{heat}^t
\end{align}
where  $V(\mathbf{A}_s)$ is the Riemannian volume of $\mathbf{A}_s$, which is parameterised with the radius of $\mathbf{A}_s$. The Neumann boundary condition corresponds to BM reflecting at the boundary. This can be approximated by pausing time and resampling the next step until it falls into the interior of $\mathbb{M}$. This estimator is asymptotically unbiased and consistent \citep{niu2019}.  Note that $t$ is the BM diffusion time.  If $t$ is large, the BM paths have higher probability to reach the neighbourhood of the target point and leads to higher covariance and vice versa. The transition density can be estimated using algorithm \ref{alg:bm}.

\begin{corollary} \label{corollary:metric}
The transition density estimation in equation \ref{k:heat} is coordinate independent. The transition density estimation is equal for different choices of coordinate system.
\end{corollary}
By Theorem \ref{theorem:sde}, it is straight forward to have corollary \ref{corollary:metric}. The transition density estimation in equation \ref{k:heat} of the stochastic process in equation \ref{eqnBM} is also coordinate independent. Based on corollary \ref{corollary:metric}, we can compare the heat kernel estimates from different parameterisation of the same manifold. Here we take the Swiss roll as an example. Numerical accuracy of estimates of heat kernel (or BM transition density) on Swiss roll are shown in Fig. \ref{fig:metriccompare}. The heat kernel is estimated using equation \ref{k:heat} at twenty nine points with equally spaced radius $\in (3,9)$ of Swiss roll (see appendix ?? for details). The diffusion time is fixed as 50. The BM trajectories are simulated by using three different metric tensors. The red solid line in Fig. \ref{fig:metriccompare} represents the heat kernel estimated based on the analytical metric, assuming the geometry of the Swiss roll is known. The derivation of the analytical metric is shown appendix ??. When the geometry of the Swiss roll is unknown, the green dash line represents the heat kernel estimated based on the metric learned by Bayesian GPLVM and the blue dot line for GPLVM. The numerical comparison of GPLVM and Bayesian GPLVM is also shown inTable \ref{tb:metriccompare}. It is clear that the Bayesian GPLVM approach achieved better performance than GPLVM and give similar estimates as the analytical metric.

\begin{figure}[htb]
 \centering
\includegraphics[width=0.5\textwidth,height=0.4\textwidth]{metricCompare.pdf}
    \caption{ \label{fig:metriccompare}
       {\footnotesize
     Comparison of estimated BM transition density using analytical metric, metric from GPLVM and metric from Bayesian GPLVM.  The red solid line represents the heat kernel estimated based on the analytical metric. The green dash lines represent the heat kernel estimated with the Bayesian GPLVM metric. The blue dot line represent the heat kernel estimated with the metric learned from GPLVM. The black dot lines are the red curve plus and minus two times standard deviation of the estimated transition density with the analytical metric. The green line locates between the two black dotted lines.
     } }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Comparison of estimates of heat kernel using metric from GPLVM and B-GPLVM.  The table shows the root mean squared differences between the BM transition density estimated using analytical metric and the BM transition density using GPLVM metric and Bayesian GPLVM. }
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em Metric learning &\em Root mean squared differences\\
\hline
GPLVM & 1.2e-3\\
B-GPLVM & 7.4e-5 \\
\end{tabular}}
\end{table}

 \begin{algorithm}[!h]
 Learn the metric $\mathcal{G}$ from the point cloud $\mathcal{S} = \{  s_i | i = 1,\ldots,n+v\ \} $.        \Comment{ \footnotesize use eqn \ref{eqn:exmetric}  }
\caption{ { Simulating BM sample paths for estimating $K_{heat}^t$ } }
\label{alg:bm}
\begin{algorithmic}
\State 1.1 Generate Brownian motion trajectories
\For{$i = 1,\ldots,n$}    \Comment{ {\small $n$ is the size of data points}  }
\For{$j = 1,\ldots,N_{BM}$}    \Comment{ {\small$N_{BM}$ is No. of trajectories} }
\For{$ l = 1,\ldots,T$}    \Comment{ {\small T steps Brownian motion, $T*\Delta t$  $\rightarrow$ max diffusion time} }
\State {\bf do}    \Comment{ {\small keep proposing x until it falls in the boundary } }
\State $q\left(x_{i,j}(l)|x_{i,j}(l-1)\right) \leftarrow \mathbb N \left( x_{i,j}(l) |  \mu \left( x_{i,j}\left(l-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right )$ \Comment{ {\small use eqn \ref{eqn:disProp} } }
\State {\bf While} $x_{i,j}$ is not in the boundary $\partial \mathbb{M}$, $Var( \phi(x_{i,j})  ) > \alpha$.   \Comment{ \footnotesize use equation \ref{eqn:varb}  }
\EndFor
\EndFor
\EndFor
\Return $\bm x$
\State 1.2 Given a discrete choice of the diffusion time $t\in \{ \Delta t, 2*\Delta t,\cdot \cdot \cdot, T*\Delta t \}$, the covariance matrix $\Sigma^{t}$ is estimated based on the BM simulation from 1.1.
\For{$i = 1,\ldots,n$}
\For{$j = 1,\ldots,n$}
\State $N_{\mathbf{A}_j}$ = which( x(t) $\in$ $\mathbf{A}_j$ )  \Comment{ {\small counting how many BM paths reach $\mathbf{A}_j$}  }
\State $K_{heat}^t(s_i,s_j)= \frac{ N_{\mathbf{A}_j}}{N_{BM}*V( \mathbf{A}_j)}$  \Comment{ {\small use eqn \ref{k:heat} } }
\State $\Sigma_{ij}^t = \sigma_h^2 K_{heat}^t(s_i,s_j)$
\EndFor
\EndFor
\Return $\Sigma^t$
\end{algorithmic}
\end{algorithm}

\subsection{Optimising the kernel hyper parameters and make predication} \label{sec:optimpar}
Given a diffusion time $t$, we can generate the covariance matrix $\Sigma_{ {\text {\bf ff}} }^{t}$ for the training data using algorithm \ref{alg:bm}. $\Sigma_{ \text{\bf ff} }^t$ can be obtained as follows:  for the $i{th}$ row of $\Sigma_{ \text{\bf ff} }^t$, $N_{BM}$ trajectories are simulated, with the starting point the $i{th}$ data point corresponding to the index of row. For each element of $\Sigma_{ \text{\bf ff} }$, $\hat{K}^t(s_i,s_j)$ is then estimated using \eqref{k:heat}. The hyperparameters can be obtained by maximising the log of the marginal likelihood. The maximum of the BM diffusion time is set as $N_t \times \Delta t$, and $\Delta t$ is the BM simulation time step as defined in \eqref{eqn:disProp}.  $N_t$ covariance matrices $\Sigma_{ {\text {\bf ff}} }^{1 \ldots N_t}$ can be generated based on the BM simulations. Optimisation of diffusion time $t$ can be done by selecting the corresponding $\Sigma_{ {\text {\bf ff}} }^{t}$ that maximises the log marginal likelihood. Estimation of $\sigma_h$ follows the standard optimisation routines, such as quasi-Newton. The log marginal likelihood function (over $f$) is given by \citep{Rasmussen2004}:
\begin{align}
\label{loglike}
p( \bm{y}| s) =  \int p( \bm{y}| {\text {\bf f} }) p({\text {\bf f} }|s) d{\text {\bf f} }   = -\frac{1}{2} \bm y^T (\Sigma_{ {\text {\bf ff}} }^{t} + \sigma_{noise}^2 I)^{-1} \bm y - \frac{1}{2} \log|\Sigma_{ {\text {\bf ff}} }^t +\sigma_{noise}^2I | - \frac{n}{2}\log2\pi.
\end{align}

\section{Simulation studies on Swiss roll} \label{simdata}
 In this section, we carry out simulation studies for a regression model with true regression functions defined on a 2-dimensional Swiss Roll embedded in $\mathbb R^3$. Examples of BM sample paths on the Swiss Roll is shown in Figure \ref{fig:BM}.

\begin{figure}[htp]
    \centering
   \subfigure[ Swiss roll in $\mathbb{R}^3$ ]{\label{u:path}  \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swiss1.png} }
       \subfigure[ Latent space and variance of mapping]{ \label{fig:swissBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swissVar.png} }
    \subfigure[Latent space and metric]{ \label{fig:torusBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{swissMag.png} }
    \caption{ \label{fig:BMexamples}
    \footnotesize{
The latent space constructed from B-GPLVM is plotted with the variance of the mapping as the background in (b), it is also plotted with the magnification factor as the background in (c).
      }
  }
\end{figure}

 \begin{figure}[htp]
    \centering
    \subfigure[True function and data points]{\label{u:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPtrue.pdf}}
     \subfigure[ \small{GP in $\mathbb{R}^3$ prediction }]{\label{u16:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGP3.pdf}}
    \subfigure[In-GP on manifold prediction]{\label{u:gp}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissmanifoldGP.pdf}}
    \subfigure[GP in $\mathbb{R}^2$ prediction]{ \label{u16:gp} \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGP2-2.pdf}}
    \caption{\label{fig:swissGP}
   { \footnotesize
     Comparison of in-GP  and Euclidean GP in the Swiss roll example. The true function and data points are shown in (a). The predictions of the GP which use the euclidean distance of data points in $\mathbb{R}^3$ are in (b). The in-GP on manifold prediction are plotted in (c). If we ignore the metric tensor and boundary, and construct a GP which use the $\mathbb{R}^2$ Euclidean distance in the latent space, the results are shown in (d).
    }
 }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Root mean squared error of regression on Swiss Roll }
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em &  \em $\mathbb{R}^3 GP$ &\em $\mathbb{R}^2 GP$ &\em  In-GP on manifold  \\
\hline
RMSE & 0.471  & 0.281 & 0.158 \\
\end{tabular}}
\end{table}

 The performance of in-GP is compared to that of an Euclidean  GP (the standard GP as in \cite{Rasmussen2004}) in $\mathbb{R}^3$.

\section{Application to image data} \label{simdata}
We consider images from the COIL data set ?? which consist of images from a fixed camera depicting objects on a motorised turntable against a black background. Each image is acquired after a 5 degree rotation. Here we consider the 66 images of object 14 (a toy cat). The camera angles are from 15 degree to 340 degree. Each image is treated as a 1024 dimensional vector.

\begin{figure}[H]
    \centering
   \subfigure[ COIL images ]{\label{fig:coilimage}  \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilimage.png} }
       \subfigure[ Latent space and variance of mapping]{ \label{fig:coilvar} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilLatentVar.png} }
    \subfigure[Latent space and metric]{ \label{fig:torusBM} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{coilLatentMag.png} }
    \caption{ \label{fig:coilMag}
    \footnotesize{
The latent space constructed from B-GPLVM is plotted with the variance of the mapping as the background in (b). The latent space is also plotted with the magnification factor as the background in (c).
      }
  }
\end{figure}

\begin{figure}[htb]
 \centering
\includegraphics[width=0.5\textwidth,height=0.4\textwidth]{coilBM.png}
    \caption{ \label{fig:coilBM}
       {\footnotesize
BM on the implicit manifold of COIL data.
     } }
\end{figure}

 \begin{figure}[H]
    \centering
    \subfigure[True function and data points]{\label{coil:true}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilTrue.png}}
     \subfigure[ \small{GP in $\mathbb{R}^{1024}$ prediction }]{\label{coil:Hd}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilHd.png}}
    \subfigure[In-GP on manifold prediction]{\label{coil:in-GP}  \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coilIn-GP.png}}
    \subfigure[GP in $\mathbb{R}^2$ prediction]{ \label{coil:hd} \includegraphics[width=0.35\textwidth,height=0.3\textwidth]{coil2D.png}}
    \caption{\label{fig:coilGP}
   { \footnotesize
     Comparison of in-GP  and Euclidean GP in COIL images example. The true function is plotted with colour in the latent space in (a). The predictions of the GP which use the euclidean distance of data points in $\mathbb{R}^{1024}$ are in (b). The in-GP on manifold prediction are plotted in (c). If we ignore the metric tensor and boundary, and construct a GP which use the $\mathbb{R}^2$ Euclidean distance in the latent space, the results are shown in (d).
    }
 }
\end{figure}

\begin{table}
\caption{\label{tb:metriccompare}Root mean squared error of regression on COIL images}
\centering
\fbox{%
\begin{tabular}{*{10}{c}}
\em &  \em $\mathbb{R}^{1024} GP$ &\em $\mathbb{R}^2 GP$ &\em  In-GP on manifold  \\
\hline
RMSE & 0.098  & 0.116 & 0.043 \\
\end{tabular}}
\end{table}

\section{Discussion} \label{discuss}
Our work proposes a novel class of intrinsic Gaussian processes on manifolds and complex constrained domains employing the equivalence relationship between heat kernels and the transition density of Brownian motions on manifolds. One of the key features of the in-GP is to fully incorporate the intrinsic geometry of the spaces for inference while respecting the potentially complex boundary or interior constraints. To reduce the computational cost of simulating BM sample paths when the sample size is large, sparse in-GPs are developed leveraging ideas from the literature on fast computation in GPs in Euclidean spaces. The results in section \ref{simdata} and \ref{real} indicate that in-GP achieves significant improvement over usual GPs. Although we did not conduct a formal asymptotic study for in-GPs, with insights gained from the Euclidean GP with squared exponential kernel in the Euclidean space (see, e.g., \cite{vandervaart}), we expect in-GPs to yield posterior consistency with respect to appropriate neighbourhoods of the true regression function $f_0$. The focus of this article has been on developing in-GPs on manifolds with known metric tensors.  There has been abundant interest in learning of unknown lower-dimensional  manifold structure in high-dimensional data.  in-GPs can be combined with these approaches for performing supervised learning on lower-dimensional latent manifolds.

\section{Introduction}

Gaussian Processes (GPs) have been widely used as data efficient modelling approaches that produce good uncertainty estimates. They also power many decision making approaches such as Bayesian optimisation, multi-armed bandits and experiment design. Characteristics of the function prior such as differentiability, periodicity and symmetry can be easily controlled by constructing a specific covariance function. Most widely used GP covariance functions are defined over Euclidean space. However, in real world applications, data often lies on a manifold within the original space. Predictions are only valid on the manifold and should only be extrapolated from observations along the manifold. For example, traffic flows can only be measured over networks of roads, surface tension is only measured on the surface of a specific object and the joints of a robot arm can only be moved safely within a manifold of the joint space of all the joints. The traffic connectivity can also easily differ from road network due to road maintenance or traffic congestion. The inference of the manifold based on data brings the accuracy and flexibility.
Directly applying a GP with a covariance function defined over Euclidean space would not be ideal as the extrapolation would not respect the geometry of the manifold.

Previously, several GP on manifold methods \citep{niu2019,lin2019,borovitskiy2020} have been proposed under the assumption that the geometry of the manifold is known. \cite{niu2019} defined a heat kernel on a manifold and constructed an intrinsic Gaussian Process on the manifold. The intrinsic GP refers to a GP that employs the intrinsic Riemannian geometry of the manifold, including the boundary features and interior conditions. \cite{lin2019} proposed extrinsic framework for GP modeling on manifolds, which relies on embedding of the manifold into a Euclidean space and then constructing extrinsic kernels for GPs on their images. \cite{dunson2020diffusion} , \cite{borovitskiy2021} and \cite{bolin2022} focused on developing GPs on graphs and metric graphs formed by data observed on the manifold.

In this work, we study GP modelling on manifolds of which the geometry is unknown. This is a more realistic setting for real world applications because measuring the exact geometry of data manifolds is often impossible or overly expensive. In particular, we focus on the scenario where only a sparse set of data points can be collected. As shown in our experiments, Graph Laplacian based methods perform poorly in this scenario due to the poor graph approximation. In contrast, we estimate the probabilistic parameterisation of the unknown manifolds using probabilistic latent variable models and propose a practical and general intrinsic GP on unknown manifolds (GPUM) methodology. This is the major novel contribution of the paper. Specifically, we investigate the geometrical structure of the unknown manifolds using Riemannian geometry. The distributions of the metric tensor and the boundaries of resulting manifolds are estimated using the Bayesian Gaussian process latent variable models (B-GPLVM). Brownian Motion(BM) sample paths on the unknown manifold are simulated using the probabilistic metric and respecting the boundaries.
The covariance kernel on the unknown manifold is estimated by employing the equivalence relationship between the heat kernel and the transition density of BM on the unknown manifold. We prove this estimator is coordinate independent.
 Our method can incorporate the intrinsic geometry and the uncertainty of the unknown manifold for inference and respect the interior constraints. With numerical experiments on synthetic and real world datasets, we compared our method against Graph Laplacian based methods and GP regressions without manifold estimations. We showed that our method is significantly better than these other methods on all the datasets.

In the following sections, concepts of Riemannian geometry are introduced in Section \ref{sec:concepts}. The metric learning algorithms are explained in Section \ref{sec:metric}. We prove the BM paths on manifolds simulated using different metrics have the same transition density in Section \ref{sec:heatBM}. The heat kernel estimates derived from the analytical metric and different metric learning methods such as Gaussian Processes Latent Variable Model (GPLVM) and B-GPLVM are compared in Section \ref{sec:density}. Applications of GPUM  on a synthetic dataset on Swiss roll, high dimensional real datasets of WiFi signals \citep{ferris2007} and COIL images \citep{nene1996} are illustrated and compared to Graph Laplacian GPs \citep{dunson2020diffusion}, Graph Mat\'{e}rn GPs \citep{borovitskiy2021} and Euclidean GPs in Section \ref{sec:simdata}, \ref{sec:wifi} and \ref{sec:real}.

\section{Concepts of Riemannian Geometry and Theoretical Background}
\label{sec:concepts}

One way of representing a high dimensional dataset is to relate it to a lower dimensional set of latent variables through a set of mapping functions (potentially nonlinear). \cite{tosi2014metric,arvanitidis2019} investigated the geometrical structure of probabilistic generative dimensionality reduction models (or latent variable models) using Riemannian metrics and computed the geodesic distances on the manifolds learned from data. A manifold embedded in a high dimensional Euclidean space can be well described by a probabilistic mapping function and the corresponding latent space. If the dimension of the latent space is the same as the intrinsic dimension of the manifold, the latent space can be interpreted as the chart of the learned manifold. Intuitively, the chart provides a distorted view of the manifold.
An illustration is shown in Fig.\ref{fig:mapping}. The $x^1$ and $x^2$ coordinates in the chart represent the radius and width of the Swiss roll in $\mathbb{R}^3$. The blue triangles in the chart can be mapped to the black points in the embedded space (Swiss roll in $\mathbb{R}^3$)  through $\phi$. In other words, the Swiss roll can be parameterised by the chart. Measurements on the manifold can be computed in the chart locally, and integrated to provide global measures. This gives rise to the definition of a local inner product, known as a Riemannian metric tensor.

Let $\mathbb{M}$ be a $q$-dimensional complete and compact  Riemannian manifold with the Riemannian metric $\mathbf{g}$, and $\partial \mathbb{M}$ its boundary.  The Riemannian metric $\mathbf{g}$ can be represented as a symmetric, positive definite matrix-valued function, which defines a smoothly varying inner product in the tangent space of $\mathbb{M}$. Let $\mathcal J$ denote the Jacobian of $\phi$. We have
\begin{align} \label{eqn:deter_metric}
\mathbf{g} = \mathcal{J}^T \mathcal{J}, \ \ \  \mathcal{J}_{i,j} = \frac{\partial \phi^i }{ \partial  x^j}
\end{align}
The superscript indicates the $j_{th}$ dimension of the chart and the $i_{th}$ dimension of the observation space.

\begin{figure}[t]

    \subfigure[Chart of $\mathbb{M}$(Swiss roll)]{\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{map.pdf} \label{fig:mapping} }
    \subfigure[Regression on $\mathbb{M}$]{\includegraphics[width=0.2\textwidth,height=0.2\textwidth]{3dregSw.png} \label{fig:regexample} }
    \caption{
    \footnotesize{
    { An illustration of a parameterisation of the manifold $\mathbb{M}$(Swiss roll) estimated from a data cloud in (a). $\phi$ maps the chart into $\mathbb{M} \subset \mathbb{R}^3$. The blue triangles in the chart are mapped to the black dots on $\mathbb{M} \subset \mathbb{R}^3$. The Riemannian geometry can help to learn the regression function in (b) (the color indicates function values)
    .}}
    }
\end{figure}

 Moreover, based on its metric tensor, $\mathbb{M}$ has an associated Laplace-Beltrami operator, which is an intrinsically defined differential operator denoted $\Delta_s$. In local coordinates, the Laplacian-Beltrami operator i
\begin{align}
\label{eqn:delta_s}
\Delta_s f = \frac{1}{\sqrt{G}} \frac{\partial}{\partial x^j} \left (    \sqrt{G} \mathbf{g}^{ij} \frac{\partial f}{\partial x^i}  \right ),
\end{align}
where $G$ is the determinant of the metric $\mathbf{g}$, $\mathbf{g}^{ij}$ is the $(i,j)$ element of its inverse and $f$ is a smooth function on $\mathbb{M}$. Take the special case where $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$, $\mathbf{g}$ becomes an identity matrix. The Laplace-Beltrami operator $\Delta_s$ becomes the Laplace operator $\Delta$ (the sum of second partial derivatives).

Consider the heat equation on $\mathbb{M}$, given by

\begin{align*}
\frac{\partial}{\partial t}K_{heat}(s_0,s,t)=\frac{1}{2}\Delta_s K_{heat}(s_0,s,t),\qquad  s_0, s\in \mathbb{M},
\end{align*}
where $s \in \mathbb{M}$, $\Delta_s$ is the Laplacian-Beltrami operator on $\mathbb{M}$, and $t \in \mathbb{R}^{+}$ is the diffusion time. A heat kernel of $\mathbb{M}$ is a smooth function $K(s_0,s,t)$ on $ \mathbb{M} \times \mathbb{M} \times \mathbb{R}^{+}$ that satisfies the heat equation. It can be interpreted as the amount of heat that is transferred from $s_0$ to $s$ in time $t$ via diffusion.

The heat kernel satisfies the initial condition $\lim_{t\rightarrow 0}K_{heat}(s_0,s,t)=\delta(s_0,s)$ with $\delta$ the Dirac delta function. The heat kernel becomes unique when we impose a suitable condition along the boundary $\partial \mathbb{M}$, such as the Neumann boundary condition: $ \partial K / \partial \mathbf{n} = 0$ along $\partial \mathbb{M}$ where $\mathbf{n}$ denotes a normal vector of $\partial \mathbb{M}$. The Neumann boundary condition can be understood as allowing no heat transfer across the boundary. If $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$,  The heat kernel has a closed form corresponding to a time-varying Gaussian function:
\begin{align*}
K_{heat}(\mathbf{s}_0,\mathbf{s},t)
=\frac{1}{(2\pi t)^{q/2}}\,
  \exp\left\{-\frac{||\mathbf{s}_0-\mathbf{s}||^2}{ 2t }\right\}, \; \mathbf{s}\in \mathbb R^q.
\end{align*}
The diffussion time $t$ controls the rate of decay of the covariance. In the following, we will also write $K_{heat}(\mathbf{s}_0,\mathbf{s},t)$ as $K_{heat}^t(\mathbf{s}_0,\mathbf{s})$.

For arbitrary Riemannian manifold, the construction of the heat kernel associated with the Laplace-Beltrami operator is not a trivial task and belongs to the fields of partial differential equations and differential geometry \citep{chavel1984}. To circumvent solving the heat equation on manifolds, \cite{niu2019} estimates the heat kernel as the BM transition density by simulating BM paths on a known manifold. In the case of point clouds, both the metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$ are unknown. In this work, we will use probabilistic latent variable models to learn the distribution of the metric $\mathbf{g}$ and define the boundary $\partial \mathbb{M}$ based on the uncertainty quantification of the mapping. The Laplace-Beltrami operator in \eqref{eqn:delta_s} is the infinitesimal generator of BM on manifolds \citep{hsu1988}. The BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations in the Ito form \citep{hsu1988,hsu2008}:
\begin{align}
\label{eqnBMhsu}
dx_i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x_j} \left(  \mathbf{g}^{ij}G^{1/2} \right) dt + \left( \mathbf{g}^{-1/2} dB(t)\right)_i,
\end{align}
where $\mathbf{g}$ is the metric tensor of $\mathbb{M}$, $G$ is the determinant of $\mathbf{g}$ and $B(t)$ represents an independent BM in the Euclidean space.

\section{Intrinsic Gaussian Processes on unknown manifolds}
\label{sec:GPmani}

In this work, we focus on the following model,
\begin{align} \label{eqn:reg}
y_i &= f(s_i) +\epsilon_i, \ \ \  \epsilon_i \sim \mathcal{N}(0, \sigma_{\epsilon}^2)
\end{align}
where  $f:\mathbb{M} \rightarrow \mathbb{R}$ is an unknown regression function. $y_i \in \mathbb{R}$ is a response variable. $s_i= (s_{i}^1,\ldots,s_{i}^p) \in \mathbb{M} \subset \mathbb{R}^p$ is an observed predictor on a complete and compact Riemannian manifold $\mathbb{M}$ embedded in $\mathbb{R}^p$.
$\mathbb{M}$ can be parameterised by a $q$ dimensional local coordinate system (chart), and $q<p$. For example, let $\mathbb{M}$ be a unit sphere and we have a point $i$ on $\mathbb{M}$. $s_i \in \mathbb{R}^3$ is a three dimensional vector representing the Cartesian coordinates of $i$. The two dimensional chart of $\mathbb{M}$ is the longitude and latitude coordinate system. The point $i$ can also be represented as a two dimensional vector $x_i$ in the chart. There exists a mapping function $\phi$ from the Cartesian coordinates to the spherical coordinates. However in the case of point clouds where $s_i$ is observed in $\mathbb{R}^p$, the parameterisation $\phi$ of $\mathbb{M}$ is unknown and $x_i$ becomes a latent variable.
We propose to infer how the output $y$ varies with the input $s$, including predicting $y$ values at new locations not represented in the training set.

A GP prior can be assigned to $f$ with a covariance function.
The choice of
covariance kernel has a fundamental impact on the results. Most common choices of covariance kernels such as the squared exponential kernel and Mat\'{e}rn kernel depend critically on the Euclidean distance between $s_i$ and $s_j$ and ignore the intrinsic geometry of $\mathbb M$. By contrast, the heat kernel depends only on the Riemmanian metric and the intrinsic geometry of $\mathbb{M}$. It provides a natural generalisation of the RBF kernel on manifolds. The heat kernel represents the diffusion of the heat on a Riemannian manifold. However it is analytically intractable to directly evaluate \citep{hsu1988}. %
\cite{niu2019} proposed a computational framework to estimate the heat kernel on the manifold in which the analytical parameterisation is known.
However, data represented as point clouds is often high dimensional and concentrated around some unknown lower dimensional structures, \cite{niu2019} is not applicable due to the lack of analytical manifold parameterisation.
In this paper, we address this problem by using the probabilistic generative dimension reduction models to learn the geometry of the implicit manifold and define the boundary. We propose to construct GPUM  by using the heat kernel of the implicit manifold as the covariance kernel. With the help of Riemannian geometry, the Brownian motion sample paths can be simulated on $\mathbb{M}$ and the heat kernel can be estimated as the transition density of the BM. With this construction we can learn the regression function on $\mathbb{M}$ as in Fig.\ref{fig:regexample}.

Let $\mathcal{D} =  \{ (s_i, y_i), i = 1, . . . , n \} $ be the data, with $n \geq 1$ the number of labeled observations, $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is the predictor and $y_i \in \mathbb{R}$ is the corresponding response. Suppose we are also given an unlabeled dataset $\mathcal{V} = \{ s_i, i=n+1 . . . ,n+v \} $ where $v \geq 1$. Consider the regression model in \eqref{eqn:reg}, we would like to make inferences about $f$ with the labeled dataset $\mathcal{D}$ and predict $y$ values for the unlabeled dataset $\mathcal{V} $. Under an GPUM  prior for the unknown regression function as $f \sim GP(0,K_{heat}^t(.,.) )$, we have
\begin{align*}
p( \mathbf{f} | s_1, s_2, ..., s_n) \sim \mathcal{N} (0, \Sigma_{\mathbf{f} \mathbf{f}}),
\end{align*}
where $\mathbf{f}$ is the discretization of $f$ over $s_1,s_2,\cdot\cdot\cdot,s_n$ so that $f_i = f(s_i)$. $\Sigma_{\mathbf{f} \mathbf{f}} \in \mathbb{R}^{n\times n}$ is the covariance matrix induced from the heat kernel. The $(i,j)$ entry of $\Sigma_{\mathbf{f}\f}$ corresponds to $\Sigma_{\mathbf{f}\f _{i,j}} = \sigma_h^2 K_{heat}^t(s_i,s_j)$. We introduce the rescaling hyperparameter $\sigma_h^2$ to add extra flexibility to the heat kernel.

\section{Related work}
One of the paramount challenges in developing GP models on manifolds is the difficulty in specifying the covariance structure via constructing valid covariance kernels on manifolds. One might hope to achieve this by replacing Euclidean norms in the squared exponential kernel with geodesic distances. However \cite{feragen2015} proved this is not generally a well-defined kernel. \cite{extrinsicGP} proposed extrinsic Gaussian Processes on manifolds by embedding the manifolds onto a higher dimensional Euclidean space. The squared exponential kernel is applied on the images of the manifold after embedding. However, such embeddings are not easy to obtain and only available for certain manifolds when the geometry is known.
\cite{yang2016} proposed a model bypassing the need to estimate the manifold, and can be implemented using standard algorithms for posterior computation in GPs.
They show that by imposing a GP prior on the regression function with a covariance kernel defined directly on the ambient space (the embedding of the manifold in a high dimensional Euclidean space), the posterior distribution yields a posterior contraction rate depending on the intrinsic dimension of the manifold. They assume that the unknown lower dimensional space where the predictors center around are a class of submanifolds of Euclidean space.  They focus on compact manifolds without boundary.

\cite{niu2019} proposed to use heat kernels to construct the intrinsic Gaussian Processes on complex constrained domains. Heat kernel can be interpreted as the transition density of Brownian Motion on the manifold \citep{hsu1988}. It is estimated by simulating BM paths on manifolds in which the analytical parameterisation is known in \cite{niu2019}.

Alternatively if the eigen-paires of the Laplacian-Beltrami operator of the manifold are available, \cite{borovitskiy2020} approximated the heat kernel with the sum of finite-many eigen-paires of the Laplacian-Beltrami. Both \cite{borovitskiy2020} and \cite{niu2019} are only applicable when the geometry of the manifold is known and the dimensions of the observation(or embedding) space are low such as $\mathbb{M} \subset \mathbb{R}^2$ or $\mathbb{M} \subset \mathbb{R}^3$ ($p= 2$ or $3$).

Most recent research in \cite{dunson2020diffusion} tackled this problem by approximating the heat kernel kernel of a compact Riemannian manifold with finitely-many eigenpairs of the Graph Laplacian using the labeled and unlabeled predictor values. We refer the Gaussian processes constructed by this approximation as the Graph Laplacian Gaussian processes (GL-GPs). Let $\Delta$ be the Laplace-Beltrami operator of a manifold and $ \{ \lambda_i \}_{i=0}^{\infty}$ be the spectrum or eigenvalues of $-\Delta$. Denote $\varphi_i$ the corresponding eigenfunction, for each $i \in N$, we have $\Delta \varphi_i = - \lambda_i \varphi_i$. If the eigen-decomposition of $\Delta$ is known, the corresponding heat kernel of the manifold has the following expression: $K_{heat}(s,s',t) = \sum_{i=0}^{\infty} e^{-\lambda_i t} \varphi_i(s) \varphi_i(s')  $. If the geometry of the manifold and the corresponding $\Delta$ are unknown, $\Delta$ is approximated  by the Graph Laplacian matrix $\bm{L}$ in \cite{dunson2020diffusion}.  $\bm{L}$ is constructed from a point cloud whose adjacency matrix is computed by using a Gaussian function with pair-wised Euclidean distance. The heat kernel is approximated as the summation of finite eigenpairs of $\bm{L}$,  $\sum_{i=0}^{n_G} e^{-\mu_i t} v_i v_i^T $, where $\mu_i$ and $v_i$ are the $i_{th}$ normalised eigenvalue and eigenvector of $\bm{L}$. The implementation of the GL-GP is provided in  Appendix \ref{ax:GL_kern}.

\cite{borovitskiy2021} leveraged the stochastic partial differential equation characterization of Mat\'{e}rn kernel to study their analog for undirected graphs and developed the Graph Mat\'{e}rn Gaussian processes (GM-GPs). The Graph Mat\'{e}rn kernel is constructed with the sum of finitely-many eigenpairs of the normalised Graph Laplacian $\bm{L}$.
$\bm{L}$ is computed from the adjacency matrix of a predefined graph. We implemented the GM-GP by following the instructions in the Github repository in \cite{borovitskiy2021}. The key component of GL-GPs and GM-GPs is the Graph Laplacian. If the true graph connections are not known, the graph constructed based on local distances such as Delaunay triangulation can be error-prone when observations are sparse \citep{hjelle2006}. The graph based methods often result in poor approximation of the manifold when the number of the observations are low. GL-GPs and GM-GPs are applied to the simulation studies and real datasets and compared to GPUM in the later sections.

There are still some critical gaps in current practices. In particular, the lack of robust methods for carrying out intrinsic statistical inference and  effective models for regressions with manifold-valued data embedded in a point cloud.
In this work we focus on learning the manifold structure using probabilistic dimension reduction methods such as Bayesian GPLVM and constructing the GPUM by using the heat kernel of the learned manifold. Other related methods such as Auto-encoders \citep{kramer1991} and Variational Auto-encoders (VAE)\citep{kingma2019}, can also be considered in this two stage approach. Auto-encoders provide a neural network based framework for learning deep latent variable models.  However, the resulting mapping function is deterministic and the model does not have a built-in quantification of its uncertainty.
There is lack of uncertainty quantification in the learned manifold. Alternatively, VAE
address this concern directly through an explicit likelihood model and a variational approximation of the representation posterior. It learns a generative model by specifying a likelihood of observations conditioned on latent variables and a prior over the latent variables.  Both the likelihood and the variational distributions have parameters predicted by
neural networks that act similarly to the encoder–decoder pair of the classic autoencoder.

\section{Learning metrics and boundaries}\label{sec:metric}
Let $\mathcal S = \{ s_i |  i=1,\cdot \cdot \cdot, n+v \}$,  $s_i \in \mathbb{R}^p$, include all predictors in the labeled dataset $\mathcal{D}$ and unlabeled dataset $\mathcal{V}$. We have $\mathcal{S} \in \mathbb{R}^{ (n+v) \times p }$. Suppose we perform probabilistic nonlinear dimensionality reduction by defining a latent variable model that introduces a set of latent (unobserved) variables $\mathcal{X}= \{x_i | i= 1,\cdot \cdot \cdot, n+v \}$, $x_i \in \mathbb{R}^q$, $\mathcal{X} \in \mathbb{R}^{(n+v)\times q}$ and $q < p$. $\mathcal X$ is related to $\mathcal S$ which is observed in a higher dimensional space. A prior distribution is placed on the latent space which induces a distribution over $\mathcal S$ under the assumption of the probabilistic mapping
\begin{align} \label{eqn:mapmodel}
s_{i}^j = \phi^j(x_i) + e_{i}^j,
\end{align}
where $x_i \in \mathbb{R}^q$ is the latent point associated with the $i_{th}$ observation $s_i \in \mathbb{R}^p$ in the original observation space and $i\leq n$. $j$ is the index of the features (dimensions) of $s$ in the observation space and $j\leq p$. $e_{i}^j$ is a Gaussian distributed noise term,   $e_i^{j} \sim \mathcal N(0,\beta^2)$.
An illustration is given in Fig. \ref{fig:mapping} where $q=2$ and $p=3$. If the mapping $\phi$ is linear and the prior $p( \mathcal X)$ is Gaussian, this model is known as probabilistic principal component analysis \citep{tipping1998}. In this work, we do not restrict to this linear assumption and consider some nonlinear dimensionality reduction methods such as GPLVM and B-GPLVM.

When the $\phi$ function in Fig.\ref{fig:mapping} is differentiable, it can be interpreted as the mapping between the latent space and the manifold $\mathbb{M}$. If the dimensions of $\mathbb{M}$ are known, by letting $q$ equal the dimensions of $\mathbb{M}$, the latent space can be interpreted as the chart of $\mathbb{M}$. If the dimensions of $\mathbb{M}$ are unknown, $q$ is estimated by using the so called Automatic Relevance Determination (ARD) in GPLVM and B-GPLVM \citep{zwiessele2017}. ARD allows assigning scaling parameters for each dimension. These scaling parameters can be incorporated to kernels such as the RBF kernel as the inverse of the squared lengthscales.

In general
a manifold may need more than one chart to be parameterised. In this work, we focus on examples of single chart. But our method can be extended to learn multiple charts from multiple datasets.
The Riemannian metric of the given model can be computed as in \eqref{eqn:deter_metric}. In the case of probabilistic LVMs, we place a Gaussian prior over the mapping function $\phi(x)|x$. The conditional probability over the Jacobian also follows a Gaussian distribution, this naturally induces a distribution over the metric tensor $\mathbf{g}$. We denote the distribution of the Jacobian as $\mathbf{J}$ given the set $\mathcal{S}$ and the mapping $\phi$ in \eqref{eqn:mapmodel}. Assuming independent rows of $\mathbf J$ \citep{lawrence2005,titsias2010},
\begin{align}
p(\mathbf J | \mathcal{S}, \Phi) = \prod_{j=1}^p \mathcal{N}( \mu_{\mathbf{J}}^j, \Sigma_{\mathbf J}  ).
\end{align}
This independent row assumption is for the dimensions in the original
observational space. The dimensions in the learned latent space are
not independent. This assumption can be relaxed by using a multi-output GP which is more compuationaly expensive \citep{alvarez2011}. The expressions of the mean $\mu_{\mathbf J}$ and variance $\Sigma_{\mathbf J} $ of the Jacobian are model specific and given in section \ref{gplvm} and \ref{bgplvm}. The resulting metric $\mathbf{g}$ follows a non-central Wishart distribution \citep{anderson1946}
\begin{align} \label{eqn:metric}
\mathbf{g} \sim \mathcal{W}_q \left(  p, \Sigma_{\mathbf{J}}, \E(\mathbf{J}^T)  \E(\mathbf{J})  \right).
\end{align}
From this distribution, the expected metric tensor can be computed as
\begin{align} \label{eqn:exmetric}
\mathcal{G} = \E( \mathbf{g} )  = \E( \mathbf{J}^T) \E( \mathbf{J} ) + p \Sigma_{\mathbf{J}}.
\end{align}
We denote the expectation by $\mathcal{G}$. Note that the variance term $\Sigma_{\mathbf{J}}$ is included in $\mathcal{G}$. It implies that the metric tensor expands as the uncertainty over the mapping increases.
Hence the BM simulation steps in the SDE in (\ref{eqnBM}) will travel `slower' in the region of the latent space where the uncertainty is high.

We also need the gradient of the expected metric to simulate BM as in section \ref{trandens}.
\begin{align}
\frac{\partial \mathcal{G}}{\partial x^l} = \frac{ \partial \E[ \mathbf{g}]}{\partial x^l} = \frac{\partial \E[\mathbf{J}^T] } { \partial x^l }  \E[\mathbf{J}] + \E[\mathbf{J}^T] \frac{ \partial  \E[\mathbf{J}]  }{\partial x^l}  + p \frac{ \partial \Sigma_\mathbf{J} }{ \partial x^l }
\end{align}
The estimates of the mapping  $\phi(x)$ become highly unreliable when $x$ is far from the data points. The corresponding metric tensor estimates and BM simulations can also be ill defined and violate the manifold geometry. To avoid these poorly estimated region due to lack of data, the boundary of the learned manifold can be defined by $\mathbb{V} ar( \phi(x) | x)$, the variance of the mapping at $x$. $\mathbb{V} ar( \phi(x) | x)$ is also a smooth function in B-GPLVM. In most cases, the $\partial \mathbb{M}$ defined here is also a $q-1$ smooth manifold. Any $x$ outside of the boundary has $\mathbb{V} ar( \phi(x) | x ) > \alpha$.
\begin{align}\label{eqn:bound}
\partial \mathbb{M} = \{  x \in \mathbb{R}^q \  \vert \ \  \mathbb{V} ar( \phi(x) | x) = \alpha   \}.
\end{align}
where $\alpha$ is defined as the maximum variance of the mapping at the shifted latent points, $\max \left( \mathbb{V} ar\left( \phi(h) | h \right) \right)$. The shifted latent points are randomly sampled from the set $\mathbb{H}=\{ h |  \| h -x_i\| = \delta_{\mathcal{X}}, x_i\in \mathcal{X}  \}$.  $\delta_{\mathcal{X}}$ can be chosen according to the maximum distance of two neighbouring latent points. In practice, we create $h$ by moving all data points with $\delta_{\mathcal{X}}$ in random directions and let $\alpha$ equal to the maximum variance of the mapping at these relocated points. The samples are chosen to make sure that $\alpha > \max \left( \mathbb{V} ar\left( \phi(x_i) | x_i \right) \right)$, for all $x_i\in \mathcal{X}$.

\subsection{ GPLVM metric}\label{gplvm}
Gaussian Process Latent Variables Model \citep[GPLVM; ][]{lawrence2005} is a nonlinear probabilistic generative model. In this section, we will derive the distribution of the metric from GPLVM.
A sample from GPLVM defines a generative mapping from $x\in \mathbb{R}^q$ in the latent space  to $s \in \mathbb{M} \subset \mathbb{R}^p$ in the observation space.
GPs define a prior over the mapping $\phi$ in \eqref{eqn:mapmodel}. A zero mean prior is used as a default choice. If the domain
knowledge of where the prior should be centred is available, such as the embedding of $\mathbb{R}^q$ into $\mathbb{R}^p$, it can also be encoded into
the mean function. The individual dimensions of the $p$ dimensional observation space are modeled independently in the GP prior sharing the same hyperparameters. Given the construction outlined above, the probability of the observed data $\mathcal{S} =\{ s_i^j| i\in \{1,\cdots,n+v \},j\in \{1,\cdots,p\}  \}$ conditioned on all latent variables $\mathcal{X} = \{ x_i^j| i\in \{1,\cdots,n+v \},j\in \{1,\cdots,q\}   \}$ is written as follows:
\begin{align}\label{eqn:gplvmJoint}
p(\mathcal{S}, \Phi |\mathcal{X} , \beta) =   p(\mathcal{S} | \Phi ,\beta ) p(\Phi | \mathcal{X})   = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j, \beta ) p( \bm{\phi}_{:}^j | \mathcal{X}),
\end{align}
 where $\Phi = \{ \bm{\phi}_i^j | i\in \{ 1,\cdot\cdot\cdot, n+v \} , j\in \{ 1,\cdot\cdot\cdot, p \} \} $, $\Phi \in \mathbb{R}^{(n+v)\times p}$ and $ \bm{\phi}_i^j = \phi(x_i)^j$.  $\mathbf{s}_{:}^j$ represents the $j_{th}$ dimension of all points in $\mathcal S$ and $\mathcal S \in \mathbb{R}^{ (n+v)\times p}$. The likelihood $p( \mathcal S | \mathcal{X} )$ is computed by marginalising out $\Phi$ and optimising the latent variables $\mathcal{X}$
\begin{align}
p(\mathcal{S} | \mathcal{X} ) = \prod_{j=1}^p \mathcal N( \mathbf{s}_{:}^j, K_{\mathcal{X} \mathcal{X}} + \beta^2 I ).
\end{align}
$K_{\mathcal{X} \mathcal{X}} $ is the $(n+v)\times (n+v)$ covariance matrix defined by the squared exponential kernel.
\cite{lawrence2005} estimated all the latent variables $\mathcal X$ and the kernel hyper-parameters of GPLVM with the maximum likelihood estimate. If the covariance kernel is differentiable, the derivative of a GP is again a GP \citep{Rasmussen2004}. This property allows us to compute the derivative of GP. The Jacobian $\mathbf{J}$ of the GPLVM mapping can be computed as the partial derivative $\frac{\partial \phi_*}{ \partial x^l}$ with respect to the $l_{th}$ dimension for any point $x_*$  in the latent space,
\begin{align} \label{eqn:Jacobian}
\mathbf{J}^{T} = \frac{\partial \phi_*}{ \partial x} = \begin{bmatrix}
 \frac{ \partial \phi(x_*)^1 }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^1}   \\
  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot \\
  \frac{ \partial \phi(x_*)^1 }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^q}
\end{bmatrix},
\end{align}
where $\frac{\partial \phi_*}{ \partial x}$ is a $q \times p$ matrix. Considering the independence across the dimensions of the observation space, the joint distribution of the $j_{th}$ dimension of the mapping $\phi$ and the $j_{th}$ column of the Jacobian can be written as
\begin{align}
\begin{bmatrix}
 \phi( \mathcal{X})^j\\
 \frac{ \partial \phi(x_*)^j }{\partial x}
\end{bmatrix}, \
\sim \mathcal N \left( 0, \begin{bmatrix}
K_{\mathcal{X},\mathcal{X}} & \partial K_{\mathcal{X},*} \\
 \partial K_{\mathcal{X},*}^T & \partial^2 K_{*,*}
\end{bmatrix}  \right ).
\end{align}
The expressions of $K_{\mathcal{X},\mathcal{X}}$ ,  $\partial K_{\mathcal{X},*}$ , and  $\partial^2 K_{*,*} $ are given in Appendix \ref{ax:GPLVM}. GPLVM provides an explicit mapping from the latent space to the observation space. This mapping defines the support of the observed data $\mathcal{S}$ as a $q$-dimensional manifold embedded in $\mathbb{R}^p$. The distribution of the Jacobian of GPLVM is the product of $p$ independent Gaussian distributions (one for each dimension of the observation space). For a point $x_*$ in the latent space, the  distribution of the Jacobian takes the form
\begin{align}
p(  \mathbf{J} | \mathcal{X}, \mathcal{S}) &= \prod^p_{j=1} \mathcal N( \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \partial K^T_{\mathcal{X},*} K^{-1}_{\mathcal{X},\mathcal{X}} \mathbf{s}_{:}^j  \ , \partial^2K_{*,*} - \partial K^T_{\mathcal{X},*} K_{\mathcal{X},\mathcal{X}}^{-1} \partial K_{\mathcal{X},*}  ). \nonumber
\end{align}
From this distribution, the expected metric tensor can be computed as in \eqref{eqn:exmetric}. The boundary $\partial \mathbb{M}$ is defined by $\mathbb{V} ar( \phi(x_*) | x_*)$, the variance of mapping in \eqref{eqn:bound}. For any $x_*$ in the latent space
\begin{align}\label{eqn:varb}
\mathbb{V} ar(\phi(x_*) | x_*) = K_{*,*} - K_{*,\mathcal{X}}K_{\mathcal{X},\mathcal{X}}^{-1} K_{\mathcal{X},*}.
\end{align}

\subsection{Bayesian GPLVM metric }\label{bgplvm}
The maximum likelihood estimation of the latent inputs  $\mathcal{X}$ in GPLVM often leads to overfitting due to its high dimensionality \citep{damianou2016}.
This overfitting can be avoided by applying a Bayesian treatment to the latent inputs. By introducing a prior distribution to the latent inputs, the marginal likelihood takes the form
\begin{align*}
p(\mathcal S)=\int p(\mathcal{S} \mid \mathcal{X}) p(\mathcal{X}) d \mathcal{X}.
\end{align*}
The  integral is intractable as the inputs $\mathcal{X}$ to the latent variable model go through a non-linear calculation in the inverse of the covariance matrix. Bayesian Gaussian Process Latent Variable Model \citep[B-GPLVM;][]{titsias2010} introduces a variational inference framework for training the latent variable model. It variationally integrates out the input variables and computes a lower bound on the exact marginal likelihood of the nonlinear latent variable model. The maximization of the variational lower bound provides a Bayesian training procedure that is robust to overfitting.

The key to the tractable variational Bayes approach is the application of variational inference to an augmented GP formulation, known as sparse GP, where the GP prior on $\phi$ is augmented to include auxiliary variables.
More specifically, we expand the conditional probabilistic model in \eqref{eqn:gplvmJoint} by including $m$ extra samples (inducing points) to the GP latent mapping
e.g. $u_i = \phi(x_{ui}) \in \mathbb{R}^p$ is such a sample \citep{lawrence2007}.
These inducing points are denoted by $\mathcal{U} = \{u_i^j |  i \in \{ 1,\cdot\cdot\cdot,m \},  j \in \{  1,\cdot\cdot\cdot,p \}    \}$, $\mathcal{U} \in \mathbb{R}^{m\times p}$ and constitute latent function evaluations at a set of pseudo-inputs $\mathcal{X}_u = \{ x_{ui}^j | i \in \{1,\ldots,m \} , j \in \{ 1,\ldots,q\} \} $, $\mathcal{X}_u \in \mathbb{R}^{m \times q}$.
The inducing inputs $\mathcal{X}_u$ are variational parameters.
$\Phi$ is defined as in \eqref{eqn:gplvmJoint}. The augmented joint probability and the marginal likelihood can be written as \eqref{eqn:full-marginal}
\begin{align}\label{eqn:full-marginal}
p(\mathcal{S} , \Phi , \mathcal{U} ,  \mathcal{X} | \beta) &=   p(\mathcal{S} | \Phi ) p(\Phi | \mathcal{U} ,\mathcal{X}) p(\mathcal {U}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j) p( \bm{\phi}_{:}^j | \mathbf{u}_{:}^j, \mathcal{X}) p(\mathbf{u}_{:}^j) p(\mathcal{X}), \\
p(\mathcal{S}) &= \int  \int  \int  p(\mathcal{S}, \Phi, \mathcal{U}, \mathcal{X}) \ d \mathcal{U}     \  d\Phi \ d \mathcal{X}.
\end{align}
With variational inference, $\log \left(p(\mathcal{S} )\right)$ can be lower bounded by applying Jensen’s inequality \citep{damianou2016}.
The resulting lower bound can be computed analytically. The full technical details of the lower bound are given in \eqref{apxeqn:vlowbou} in Appendix \ref{ax:BGPLVM}. The Jacobian has the same shape as in \eqref{eqn:Jacobian}. For a point $x_*$ in the latent space, the distribution of the Jacobian for B-GPLVM takes the form:
\begin{align}
p(  \mathbf{J} | \mathcal{X},\mathcal{U},\mathcal{S}) &= \prod^p_{j=1} \mathcal N( \  \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \  \partial K^T_{\mathcal{X}_u,*} K^{-1}_{\mathcal{X}_u,\mathcal{X}_u} \mu_{qu}^j , \  \partial^2K_{*,*} - \partial K^T_{\mathcal{X}_u,*} \Lambda \partial K_{\mathcal{X}_u,*}  )  \nonumber   \\
 \Lambda &= K_{\mathcal{X}_u\mathcal{X}_u}^{-1} - K_{\mathcal{X}_u\mathcal{X}_u}^{-1} \Sigma_{qu} K_{\mathcal{X}_u\mathcal{X}_u}^{-1} \nonumber  \\
\mathbb{V} ar(\phi(x_*) | x_*) &= K_{*,*} - K_{*,\mathcal{X}_u} \Lambda K_{\mathcal{X}_u,*}\label{eqn:varb}
\end{align}
where $\mu_{qu}^j$ and $\Sigma_{qu}$ are the mean and variance of the variational distribution of the inducing points $\mathcal{U}$. The expressions of $\mu_{qu}$ and $\Sigma_{qu} $ are given in \eqref{apxeqn:vquqsig} in Appendix \ref{ax:BGPLVM}. The expectation of the metric tensor can be computed as in \eqref{eqn:exmetric}. The boundary $\partial \mathbb{M}$ of the implicit manifold can be defined by computing the variance of the mapping as in \eqref{eqn:varb}.

\section{Estimate heat kernel as BM transition density} \label{sec:heatBM}

 In this section, we will estimate the BM transition density by simulating BM sample paths on the implicit manifold $\mathbb{M}$.

\subsection{ Simulating Brownian motion on implicit manifolds }
\label{trandens}
From section \ref{sec:metric}, we learned the probabilistic parameterisation $\phi$ of $\mathbb{M}$, the distribution of the associated metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$. In order to estimate the BM transition density on $\mathbb{M}$, we first need to simulate BM trajectories. Simulating the sample paths of BM on $\mathbb{M} \subset \mathbb{R}^p$ is equivalent to simulating the stochastic processes in the latent space (or chart) in $\mathbb{R}^q$, $q<p$. An illustrated example is shown in Fig. \ref{fig:BM}.

\begin{figure}[t]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{bmlatent1.pdf}
    \caption{ \label{fig:BM}
    \footnotesize{
    {A BM sample path (blue line, right panel) on $\mathbb{M}$ (Swiss roll in $\mathbb{R}^3$) and its equivalent stochastic process (purple line, left panel) in the chart in $\mathbb{R}^2$. $\phi: \mathbb{R}^2 \rightarrow \mathbb{M} \subset \mathbb{R}^3$ is a parametrization of $\mathbb{M}$. The red dot is the starting location of the BM trajectory. The horizontal axis of the latent space represents the radius of Swiss roll. The vertical axis is for the width. The gray color in the latent space indicates the magnification factor. When the latent space is mapped to $\mathbb{M} \subset \mathbb{R}^3$, the darker region will be stretched more.}}
    }
\end{figure}
BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations (SDE) in Ito form \citep{hsu1988}. We use the expected metric $\mathcal{G} = \E[\mathbf{g}]$ to construct the SDEs
\begin{align}
\label{eqnBM}
dx^i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x^j} \left(  \mathcal{G}^{ij}G^{1/2} \right) dt + \left( \mathcal{G}^{-1/2} dB(t)\right)_i
\end{align}
where $x^i$ represents the $i_{th}$ dimension of the latent space (chart). $\mathcal{G}$ is defined in eqn (\ref{eqn:exmetric}),
$G$ is the determinant of $\mathcal{G}$ and $B(t)$ represents an independent BM in Euclidean space. The discrete form of \eqref{eqnBM} is derived in \eqref{disBM}.
\begin{align}
\label{disBM}
x^i(t) &= x^i(t-1) + \frac{1}{2} \sum^{q}_{j=1} \left(  -\mathcal{G}^{-1} \frac{\partial \mathcal{G}}{\partial x^j} \mathcal{G}^{-1} \right)_{ij} \Delta t + \frac{1}{4} \sum_{j=1}^q(\mathcal{G}^{-1})_{ij} tr\left(\mathcal{G}^{-1}\frac{\partial \mathcal{G}}{\partial x^j}\right) \Delta t + \left( \mathcal{G}^{-1/2} dB(t)\right)_i \nonumber \\
         &= \mu( x^i(t-1),\Delta t)_i + \left( \sqrt{\Delta t} \mathcal{G}^{-1/2} z^q\right)_i
\end{align}
where $\Delta t$ is the diffusion time of the BM simulation step and $z^q$ is a $q$-dimensional standard normal random variable. The discrete form of the SDE also defines the proposal mechanism of the BM
\begin{align}
\label{eqn:disProp}
q\left(x(t)|x(t-1)\right) = \mathbb N \left( x(t) |  \mu \left( x\left(t-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right ).
\end{align}

\begin{theorem} \label{theorem:sde}
The stochastic process defined in \eqref{eqnBM} is coordinate independent. With a given $\delta t$, simulations in any choice of local coordinates (or metric) as above are equivalent to the same step in $\mathbb{M}$.
\end{theorem}

The proof of Theorem \ref{theorem:sde} is given in Appendix \ref{ax:independent}. This implies that the BM sample paths simulated from different choices of metric $\mathcal{G}$ should have the same properties. The BM steps are sampled from the proposal distribution in \eqref{eqn:disProp}, which is defined by the metric tensor. The boundary $\partial \mathbb{M}$ of the manifold is also quantified by the uncertainty of the mapping. We apply the Neumann boundary condition as in Section \ref{sec:GPmani}. As a result the simulated sample paths stay within the boundary. An example of BM trajectory on Swiss roll is given in Fig. \ref{fig:BM}. The latent space and the associated metrics are learned from B-GPLVM. The gray color in the latent space indicates the square root of the determinant of the metric. It is also known as the magnification factor \citep{bishop1997,zwiessele2017,tosi2014}.
\begin{align*}
\mathcal{MF} = \sqrt{det(\mathcal{G}) }
\end{align*}
The geometric interpretation of the magnification factor is how much a small piece of the latent space in $\mathbb{R}^q$ will be stretched or compressed when it is mapped to $\mathbb{M} \subset \mathbb{R}^p$. For example in Fig.\ref{fig:BM} left panel, the horizontal axis of the latent space can be interpreted as the scaled radius of the Swiss roll, the vertical axis as the width of the Swiss roll. When the radius is bigger, the corresponding area in the latent space is darker and the magnification factor is larger. When the latent space is mapped back to the manifold, the darker region will be stretched more. The purple trajectory of the stochastic process in the latent space (Fig. \ref{fig:BM} left panel) is denser in the darker area and more spread out in the bright area. The stochastic process travels with smaller steps when the magnification factor is large and vice versa.
Note that, due to the BM being coordinate independent as in Theorem \ref{theorem:sde}, the BM trajectory is evenly spread out in the manifold of the Swiss roll in $\mathbb{R}^3$ (see Fig. \ref{fig:BM} right panel).

\subsection{Estimate the transition density of BM} \label{sec:density}
Considering the BM $\{ \mathbf{S}(t) | t>0 \}$ on $\mathbb{M} \subset \mathbb{R}^p$. The BM starts from $\mathbf{S}(0)= s_0$ at time $0$. We simulate $N_{BM}$ sample paths. Given a point $s \in \mathbb{M}$, we define a small neighbourhood of $s$ as $\mathbf{A}_s \in \mathbb{M}$. For any $t > 0$, the probability of $\mathbf{S}(t)$ reaching $\mathbf{A}_s$ at time t, $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$, can be approximated by
\begin{align}
\label{eqn:BMdensity}
p( \mathbf{S}(t)\in \mathbf{A}_s\,| \ \mathbf{S}(0)=s_0)  \approx  \frac{N_{A_s}}{N_{BM}}
\end{align}
where $N_{A_s}$ is the number of sample paths reaching $\mathbf{A}_s$ at time $t$. An illustrative diagram is shown in Fig. \ref{fig:BM-density}. This BM transition probability is defined as the integral of the BM transition density over $\mathbf{A}_s$ which is also the heat kernel $K_{heat}^t$.
Since we do not have the analytical expression of the transition probability,  the transition density cannot be derived by taking the derivative of $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$. Instead, $K_{heat}^t$ can be numerically approximated as $ \hat{K}_{heat}^t $
\begin{align}
\label{k:heat}
K_{heat}^t(s_0,s) \approx  \frac{p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0 )}{V(\mathbf{A}_s)} \approx \frac{1}{V(\mathbf{A}_s)}\cdot \frac{N_{A_s}}{N_{BM}} = \hat{K}_{heat}^t
\end{align}
where  $V(\mathbf{A}_s)$ is the Riemannian volume of $\mathbf{A}_s$. $V(\mathbf{A}_s)$ is parameterised by its radius $\omega$. \cite{niu2019} has provided a indication of the optimal order of magnitude of $\omega$ by minimizing the error of the kernel estimator. When $V(A_s)$ is large, the error of estimating the transition probability becomes smaller. But the error of approximating the transition density become bigger. The former is called Monte Carlo error and the later is called numerical error in Niu et. al (2019). The optimal order of magnitude of $\omega$ can be derived by balancing these two errors.

\begin{figure}[t]
 \centering
\includegraphics[width=0.45\textwidth,height=0.4\textwidth]{BM-density-example.png}
    \caption{ \label{fig:BM-density}
       {\footnotesize
    Three BM paths (red,blue,green) are simulated on the Swiss Roll, with the starting point at $s_0$ (black ball). $s$ (red ball) is the target point. $\bm{A}$ is the neighbourhood of $S$. Only the red path reach $\bm{A}$ at time $t$. The transition probability $p\{S(t)\in \bm{A} | S(0)=s_0 \}$ is 1/3.
     } }
\end{figure}

The Neumann boundary condition corresponds to BM reflecting at the boundary. This can be approximated by pausing time and resampling the next step until it falls into the interior of $\mathbb{M}$. This estimator is asymptotically unbiased and consistent \citep{niu2019}.  Note that $t$ is the BM diffusion time.  If $t$ is large, the BM paths have a higher probability to reach the neighbourhood of the target point leading to higher covariance and vice versa. The transition density can be estimated using Algorithm \ref{alg:bm}.

\begin{corollary} \label{corollary:metric}
The transition density estimate in \eqref{k:heat} is coordinate independent.
\end{corollary}
By Theorem \ref{theorem:sde}, it is straightforward to have Corollary \ref{corollary:metric}. The transition density estimate in \eqref{k:heat} of the stochastic process defined in \eqref{eqnBM} is also coordinate independent.
Based on Corollary \ref{corollary:metric},  we can evaluate the heat kernel of the implicit manifold estimated by different LVMs, independent of the specific parameterisations of their latent spaces.
If the estimated manifold is close to the true manifold, we expect the resulting BM transition density estimates to be similar to the ones estimated using the analytical parameterisation.

Here we take the Swiss roll as an example. Assuming the geometry of the Swiss roll is known, we can follow \cite{niu2019}'s approach and evaluate the heat kernel $K_{heat}^t(s_0,s)$ by simulating BM paths with the analytical metric tensor. The derivations of the analytical metric and parameterisation of the Swiss roll are shown in Appendix \ref{ax:swisspar}. Let $s_0$ with $radius=6$ and $width=3$ be the starting point of the BM. $N_{BM}=20000$ BM paths are simulated. The BM transition density is evaluated using \eqref{k:heat} at twenty nine target points $\{ s_j \in \mathbb{M} \subset \mathbb{R}^3| j \in \{1, \cdots,29 \} \}$ in the observation space. These target points are centred on $s_0$ and equally spaced. The diffusion time is fixed at 50. The results are plotted as the red solid line in Fig. \ref{fig:swiss-gl-kern}. The horizontal axis is the radius of the Swiss roll and the vertical axis is the transition density. The red density plot is asymmetric. When the radius is large, the transition density estimate decreases more quickly.

If the geometry of the Swiss roll is unknown,  Bayesian GPLVM and GPLVM can be applied to learn the metrics from 250 grid points on the Swiss roll. The derivations of B-GPLVM metrics and GPLVM metrics are given in section \ref{bgplvm} and \ref{gplvm}. The BM trajectories are simulated using the estimated metric tensors. The heat kernel estimates using B-GPLVM metrics are plotted as the green dashed line in Fig.\ref{fig:swiss-gl-kern}. It is clear that the B-GPLVM results are very close to the red solid line.
The estimates using GPLVM metrics are plotted in brown dashed line. It is also close to the solid red line, but not as good as the B-GPLVM estimates. GPLVM performs point estimate of $x$ in the latent space while B-GPLVM estimates a Bayesian posterior of $x$. As a result, B-GPLVM is more robust in terms of estimating the latent variables \citep{damianou2016}.

The heat kernel estimates from the Graph Laplacian (GL) approach (\cite{dunson2020diffusion}) in different data regime are also provided in Fig. \ref{fig:swiss-gl-kern}. When the number of points on the Swiss roll is 250, the GL kernel estimates are plotted as the blue dashed line. It is far from the solid red line and do not match the overall pattern of the analytical kernel estimates. When the number of grid points is increased to 1000 and 10000, the GL kernel estimates are plotted as the purple dashed line and black dash dotted line. The estimates are closer to the solid red line as the number of grid points increases. Comparing to the GL approach, the kernel estimates using B-GPLVM metrics achieve the best performance with much fewer points on the Swiss roll. B-GPLVM is used in the simulation study and real data applications in later sections.

\begin{figure}[h]
 \centering
\includegraphics[width=0.6\textwidth,height=0.5\textwidth]{GLkernel-c.pdf}
    \caption{ \label{fig:swiss-gl-kern}
       {\footnotesize
Comparison of heat kernel estimates using the analytical metric, B-GPLVM metric, GPLVM metric and Graph Laplacian.  The red solid line represents heat kernel estimates using the analytical metric. The green dashed line represents estimates using B-GPLVM metric. The brown dash dotted line represents estimates using GPLVM metric. The heat kernel estimates from the Graph Laplacian approach are plotted as the blue dashed line. We increase the number of the grid points on the Swiss roll from 250 to 1000 and 10000. The GL estimates are plotted as the purple dotted line and black dashed line.
     } }
\end{figure}

\begin{algorithm}[tb]
   \caption{ Simulating BM sample paths on $\mathbb{M}$ for estimating $K_{heat}^t$ }
   \label{alg:bm}
\begin{algorithmic}
  \STATE Learn the metric $\mathcal{G}$ from the point cloud $\mathcal{S} = \{  s_i | i = 1,\ldots,n+v\ \} $. \COMMENT {\footnotesize use eqn \eqref{eqn:exmetric}  }
   \STATE  1.1 Generate BM trajectories on implicit manifold.
   \FOR{$i = 1,\ldots,n$   \COMMENT{ \footnotesize $n$ is the size of data points } }
     \FOR{$j = 1,\ldots,N_{BM}$   \COMMENT{\footnotesize $N_{BM}$ is No. of trajectories }  }
    \FOR{$ l = 1,\ldots,N_t$   \COMMENT{\footnotesize $N_t$ steps Brownian motion, $N_t\times \Delta t$  $\rightarrow$ max diffusion time} }

 \STATE {\bf do}   \COMMENT{\footnotesize keep proposing x until it falls inside of the boundary }
   \STATE $q\left(x_{i,j}(l)|x_{i,j}(l-1)\right) = \mathbb N \left( x_{i,j}(l) |  \mu \left( x_{i,j}\left(l-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right )$  \COMMENT{ \footnotesize  use eqn \eqref{eqn:disProp} }
\STATE {\bf While} $\mathbb{V} ar( \phi(x_{i,j})  ) > \alpha$, $x_{i,j}$ is outside of the boundary $\partial \mathbb{M}$.   \COMMENT{ \footnotesize use eqn \eqref{eqn:varb}  }

    \ENDFOR
     \ENDFOR
    \ENDFOR
    \RETURN $\bm x$
\STATE 1.2 Given a discrete choice of the diffusion time $t\in \{ \Delta t, 2\Delta t,\cdot \cdot \cdot, N_t\Delta t \}$, the covariance matrix $\Sigma^{t}$ is estimated based on the BM simulation from Algorithm 1.1.
\FOR{$i = 1,\ldots,n$}
\FOR{$j = 1,\ldots,n$}
\STATE $N_{\mathbf{A}_j}$ = which( x(t) $\in$ $\mathbf{A}_j$ )  \COMMENT{ {\small counting how many BM paths reach $\mathbf{A}_j$}  }
\STATE $K_{heat}^t(s_i,s_j)= \frac{ N_{\mathbf{A}_j}}{N_{BM}*V( \mathbf{A}_j)}$  \COMMENT{ {\small use eqn \eqref{k:heat} } }
\STATE $\Sigma_{ij}^t = \sigma_h^2 K_{heat}^t(s_i,s_j)$
\ENDFOR
\ENDFOR
\RETURN $\Sigma^t$
\end{algorithmic}
\end{algorithm}

From  \eqref{k:heat} we can see the construction of GPUM and the heat kernel requires simulating BM sample paths at each data point. Although the BM
simulations are trivially parallelizable, the computational cost can be high when the
number of data points is large. \cite{niu2019} proposed the sparse intrinsic GP on known manifold by introducing some inducing points. The number of inducing
points is much smaller than the number of data points. BM paths only need to be
simulated starting at the inducing points instead of every data point. The inducing point approximation summarizes the training data into a small set of inducing points, so that
inference could be done more efficiently \citep{QuioneroCandelaRasmussen2005}. Similar approach can be applied in the
GPUM when the manifold is unknown. Small number of inducing
points can be introduced in the learned latent space. The focus of this paper has been on developing the GPUM. The development of the sparse GPUM is for the future research.

\subsection{Optimising the kernel hyperparameters} \label{sec:optimpar}
Given a diffusion time $t$, we can generate the covariance matrix $\Sigma_{ {\text {\bf ff}} }^{t}$ for training data $\mathcal{D}$ using Algorithm \ref{alg:bm}. $\Sigma_{ \text{\bf ff} }^t$ can be obtained as follows: with the $i_{th}$ data point as the starting point, $N_{BM}$ trajectories are simulated to generate the $i_{th}$ row of $\Sigma_{ \text{\bf ff} }^t$. For each element of $\Sigma_{ \text{\bf ff}}^t$, $\hat{K}^t(s_i,s_j)$ is then estimated using  \eqref{k:heat}. The hyperparameters can be obtained by maximizing the log of the marginal likelihood  (over $f$) in \eqref{loglike}.
The maximum BM diffusion time is set as $N_t \times \Delta t$. $\Delta t$ is the BM simulation time step as defined in \eqref{eqn:disProp}. $N_t$ is the number of simulation steps. $N_t$ covariance matrices $\Sigma_{ {\text {\bf ff}} }^{1 \ldots N_t}$ can be generated based on the BM simulations. Optimisation of the diffusion time $t$ can be done by selecting the corresponding $\Sigma_{ {\text {\bf ff}} }^{t}$ which maximizes the log marginal likelihood. Estimation of $\sigma_h$ follows standard optimisation routines, such as quasi-Newton.
\begin{align}
\label{loglike}
\log p( \bm{y}| s) &= \log \int p( \bm{y}| {\text {\bf f} }) p({\text {\bf f} }|s) d{\text {\bf f} }   \nonumber \\
&= -\frac{1}{2} \bm y^T (\Sigma_{ {\text {\bf ff}} }^{t} + \sigma_{noise}^2 I)^{-1} \bm y - \frac{1}{2} \log|\Sigma_{ {\text {\bf ff}} }^t +\sigma_{noise}^2I | - \frac{n}{2}\log2\pi.
\end{align}

Let ${\bf f}_*$ be a vector of values of $f(\cdot)$ at the unlabeled points in $\mathcal{V}$. Under the regression model in \eqref{eqn:reg} and the GPUM  prior, we have the joint distribution of ${\bf y} $ and ${\bf f}_*$:
\begin{align}
p( {\text {\bf y}}, {\text {\bf f}}_*) = \mathcal{N} \left ( 0,
 \left [ \begin{array}{cc}
 \Sigma_{{\text {\bf f}} {\text {\bf f}}} + \sigma_{noise}^2 I_n &\Sigma_{ {\text {\bf f}}  {\text {\bf f}}_*} \\[0.3em]
 \Sigma_{ {\text {\bf f}}_* {\text {\bf f}}} &\Sigma_{ {\text {\bf f}}_* {\text {\bf f}}_*} \\[0.3em]
                \end{array} \right ]
                \right ),
\end{align}
where $\Sigma_{{\text {\bf f}}_*{\text {\bf f}}}$ is the covariance matrix for training(labeled) and unlabeled data points. The predictive distribution is derived by marginalising out {\bf f}:
\begin{align}
p( {\text {\bf f}}_*| \bm y) = \int p( {\text {\bf f}}_*{\text {\bf f}}| \bm y) d{\text {\bf f}} = \mathcal N \left ( \Sigma_{ {\text {\bf f}}_*{\text {\bf f}}}  \left(\Sigma_{{\text {\bf f}}{\text {\bf f}}}  + \sigma_{noise}^2 I \right)^{-1}\bm y , \ \  \Sigma_{{\text {\bf f}}_*{\text {\bf f}}_*}-\left(\Sigma_{{\text {\bf f}}{\text {\bf f}}} +\sigma_{noise}^2 I \right)^{-1} \Sigma_{{\text {\bf f}}{\text {\bf f}}_*} \right)
\end{align}

\section{Simulation study on Swiss roll} \label{sec:simdata}

 In this section, we carry out a simulation study for a regression model with synthetic data on the Swiss roll which is a two dimensional manifold depicted by a point cloud in $\mathbb{R}^3$. The point cloud is plotted in Fig. \ref{swiss:point3d}. It is comprised of the set of labeled points $n=24$ and the set of unlabeled points $v=450$. Both labeled and unlabeled observed points are used in B-GPLVM to learn the latent space. The point cloud is unfolded into a flat surface as in Fig. \ref{swiss:var}. The unlabeled points are plotted as blue triangles and the labeled points are in red in the latent space. The variance of the mapping is plotted as the gray background of the latent space. It is clear from Fig. \ref{swiss:var} that the background color gets darker in the region further away from the blue triangles. This indicates the variance of the mapping
gets bigger in the region which is far from the observations.

 Based on the definition in \eqref{eqn:bound}, the boundary $\partial \mathbb{M}$ is plotted in Fig. \ref{swiss:boud}. The black regions on the left and right sides of Fig. \ref{swiss:boud} are outside of $\partial \mathbb{M}$. The white area in the middle is within $\partial \mathbb{M}$. The magnification factor is plotted as the background color in Fig. \ref{swiss:mag}. Similar to Fig. \ref{fig:BM}, the horizontal axis can be interpreted as the scaled radius of the Swiss roll.
 When the radius is bigger, the corresponding magnification factor is larger (the color is darker).  However, as there are fewer data points available at the tail of the Swiss roll, the estimation of the implicit manifold from B-GPLVM becomes less accurate. This can be observed on the right end of Fig. \ref{swiss:mag} where the last column of the observed latent points is further away from the rest. Once the metric $\mathcal{G}$ and the boundary $\partial \mathbb{M}$ are learned from B-GPLVM, we can estimate the heat kernel by simulating BM paths on the implicit manifold. An example of a BM sample path on Swiss roll is shown in Fig. \ref{fig:BM}.

For the labeled points, the response variables are
 \begin{align*}
y_i = f(s_i^1,s_i^2,s_i^3) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^3
 \end{align*}

where $f$ is the unknown regression function, $s_i$ is the coordinate of the observed point in $\mathbb{M} \subset \mathbb{R}^3$. For better visualisation, the true function is plotted in the unfolded Swiss roll in Fig.\ref{swiss:true} using a two dimensional analytical parameterisation. The coordinates are radius and width. 24 labeled observations are marked as black crosses. The true function values are indicated by the color codes and contours. The regression function varies slowly when the radius is small and rapidly when the radius is big.

We first apply the Euclidean $\mathbb{R}^3$ GP (the standard GP as in \cite{Rasmussen2006} Chapter 4) constructed by the squared exponential kernel with $\mathbb{R}^3$ Euclidean distance. With this kernel setting, the $\mathbb{R}^3$ GP completely ignores the interior structure of the manifold and lets the inner layer and outer layer of the Swiss roll interact. The predictive means on the unlabeled points are shown in Fig. \ref{swiss:3d}. Compared to the true function in Fig.\ref{swiss:true}, the overall shape of $\mathbb{R}^3$ GP prediction contours is more wiggly. The color coding of Fig.\ref{swiss:3d} is also very different in the regions where the radius (horizontal axis) is 6, 8 and 10. In the second model, the $\mathbb{R}^2$ Euclidean distance in the latent space is used in the squared exponential kernel to construct the $\mathbb{R}^2$ GP. The geometric properties such as the metric and magnification factors are ignored. The color coding of the predictive means for $\mathbb{R}^2$ GP is shown in Fig. \ref{swiss:r2} in Appendix \ref{ax:swisspar}. Since the regression function is nonstationary in the latent space, the $\mathbb{R}^2$ GP is underfitting. The GPUM  predictive mean is shown in Fig. \ref{swiss:gp}. The overall pattern of the  GPUM  prediction is similar to the true function. The shape of the contours and the color coding of Fig.\ref{swiss:gp} are consistent with Fig.\ref{swiss:true}. The prediction of the GL-GP is shown Fig.\ref{swiss:gl} in Appendix \ref{ax:swisspar}. Since the GL kernel estimates in Fig. \ref{fig:swiss-gl-kern} are far from the analytical kernel, the GL-GP prediction is also poor. A one dimensional comparison is generated by plotting the predictions of all methods at $z=4$ and varying the radius from 2 to 12.5 in Fig.\ref{swiss:compare1d}. It is clear the $\mathbb{R}^2$ prediction(brown dotted line) is under fitting. The $\mathbb{R}^3$ prediction(green dashed line) is oscillating in the opposite direction of the ground truth (blue solid line). The GL-GP prediction (black dashed line) is also oscillating around the blue solid line. The GPUM  prediction (red dashed line) achieves the best performance and follows the overall pattern of the ground truth.

\begin{figure}[H]
    \centering
\subfigure[Swiss roll point cloud in $\mathbb{R}^3$ ]{\label{swiss:point3d}  \includegraphics[width=0.44\textwidth,height=0.3\textwidth]{swiss1.png} }
\subfigure[ Latent space and variance of mapping]{ \label{swiss:var} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{swissVar1.png} }
\subfigure[Latent space and magnification factor]{ \label{swiss:mag} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{swissMag1.png} }
\subfigure[Boundary]{ \label{swiss:boud} \includegraphics[width=0.49\textwidth,height=0.35\textwidth]{swissBound1.png} }
    \caption{ \label{fig:swissDR}
    \footnotesize{ (a) The Swiss roll represented as a point cloud in $\mathbb{R}^3$. (b) The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and red crosses are labeled points. The background color represents the variance of the mapping, a dark background represents high uncertainty. (c) The same latent space visualization with the background color representing the magnification factor, a dark background represents high magnification factor. (d) The same latent space visualization highlighting the boundary. The dark region is outside of the boundary and the white region is inside of the boundary.}
  }
\end{figure}

In order to evaluate the performance of different methods in different data regimes, we consider three scenarios with different numbers of unlabeled grid points on the Swiss roll from $v=250$ to $v=450$ and $v=800$. These unlabeled grid points are used to estimate the manifold for GPUM , GM-GP, GL-GP and are also used as the test set of the regression accuracy. We constructed twenty sets of training points by randomly selecting $n=23$ labeled points from the labeled set. For each training set, GPUM , GM-GP, GL-GP, $\mathbb{R}^2$GP and $\mathbb{R}^3$GP have been applied to make predictions at the points in the test set. The root mean square errors (RMSE) are calculated between the true function values and the predictive means of all five models. The mean and standard deviation of the root mean square errors are reported in Table \ref{tb:swiss}. GPUM’s results are consistent
across all three scenarios and significantly better than all other methods. The GL-GP and GM-GP have similar performance. Both models perform better when the number of the grid points is large. They are significantly better than the Euclidean GPs. Comparing to the graph based approaches, GPUM achieves the minimum mean RMSE with fewer points on the manifold.

 \begin{figure}[H]
    \centering
    \subfigure[True function and data points]{\label{swiss:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPtrue.pdf}}
     \subfigure[ \small{$\mathbb{R}^3$ GP  prediction }]{\label{swiss:3d}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPr3.pdf}}
    \subfigure[GPUM  prediction]{\label{swiss:gp}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissmanifoldGP.pdf}}
    \subfigure[ Comparison of predictions when fixing z ]{ \label{swiss:compare1d} \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{fig/SwissReg-z4new.pdf}
    \caption{\label{fig:swiss2d}
   { \footnotesize
(a) The true function in the unfolded Swiss roll. The labeled points are marked with black crosses. (b) The prediction of $\mathbb{R}^3$ GP. $\mathbb{R}^3$ GP is constructed using the Euclidean distance in $\mathbb{R}^3$. (c) The prediction of GPUM . (d) Comparing the predictions of all methods on a spiral in the swiss roll by fixing the z coordinate at 4 and changing the radius from 2 to 12.5. The vertical axis represent the value of the prediction. The horizontal axis represents the radius.
    }
 }
\end{figure}

\begin{table}[h]
\caption{\label{tb:swiss}Comparison of the root mean squared errors of five methods on Swiss roll. Values in parentheses show the standard deviation.}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
    &   $\mathbb{R}^3 GP$ & $\mathbb{R}^2 GP$ &  GPUM  & GL-GP & GM-GP\\
\hline
mean RMSE $v=250$ &0.284(0.006) & 0.293(0.005)&0.163(0.020)& 0.243(0.003)&0.231(0.001) \\
mean RMSE $v=450$ & 0.298(0.007)  & 0.290(0.005) & 0.162(0.003) &0.220(0.002) & 0.207(0.002)\\
mean RMSE $v=800$ &0.287(0.006)&0.282(0.005)& 0.164(0.002)&0.216(0.001) & 0.206(0.001)\\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\begin{figure}[H]
    \centering
\subfigure[ Latent space and uncertainty of mapping ]{\label{wifi:latent}  \includegraphics[width=0.44\textwidth,height=0.3\textwidth]{wifi_Latent.png} }
\subfigure[ Latent space and boundary]{ \label{wifi:bound} \includegraphics[width=0.51\textwidth,height=0.33\textwidth]{wifi_bound.png} }
\subfigure[Latent space and magnification factor]{ \label{wifi:mag} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{wifi_mag.png} }
\subfigure[A BM trajectory]{ \label{wifi:bm} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{wifi_BMmag.png} }
    \caption{ \label{fig:wifi}
    \footnotesize{ (a)
The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and the red crosses are labeled points. The background color represents the variance of the mapping, a dark background represents high uncertainty. (b) The same latent space visualization highlighting the boundary. The dark region is outside of the boundary and the white region is inside of the boundary. (c) The same latent space visualization with the background color representing the magnification factor, a dark background represents a high magnification factor. (d) A BM trajectory is shown in purple. The red star is the starting location.}
  }
\end{figure}

\section{Location estimation from WiFi signal} \label{sec:wifi}

Indoor location estimation has tremendous value but standard location estimation techniques such as Global Positioning System (GPS) do not work indoors. Instead, indoor wireless signals from devices such as WiFi access points can be exploited for location estimation. In this section, we consider the problem of indoor 2D location estimation from WiFi access point signal strengths. We use the WiFi data collected by \cite{ferris2007}, in which a series of WiFi signal strength traces are collected by a mobile device which travels in a one floor university building. The 2D location coordinates of the mobile device are also recorded by a click-to-map based annotation program. The total number of WiFi access points in this dataset is 30. As it is often expensive to collect labeled data, we mimic a low data indoor location estimation scenario by assuming that the true locations of only three points are known and we aim at predicting the indoor locations of the remaining points based on the WiFi signals.

This location estimation problem can be treated as a regression problem, in which the location coordinate $y$ of the mobile device can be modelled as a function of the high dimensional WiFi signal
\begin{align*}
y_i = f(s_i) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^{30}
\end{align*}
where $f$ is the unknown regression function and WiFi signal $s_i$ is represented by a 30 dimensional vector. Here we consider the WiFi signal measurements at $n+v=36$ locations. Only $n=3$ of the locations are labeled with one dimensional location coordinates of the mobile device. To avoid selecting the three points clustered together, we randomly pick the points from three different regions respectively, where the union of the regions cover the whole dataset. Different methods are applied to estimate the coordinates of the mobile device in the testing sets (the unlabeled $v=33$ locations). Twenty training and testing sets are generated from this random selection.

Unlike the simulation study of the Swiss roll, we cannot plot the high dimensional point cloud of the WiFi signals. The implicit manifold is also unknown. We first estimate a $q=2$ dimensional latent space (the chart of the underlying manifold) using B-GPLVM. The value of $q$ is determined by the ARD contributions which measure how much each dimension is contributing to the latent space. The input dimensions are sorted based on the relevance assigned by the scaling in B-GPLVM. The plot of ARD contributions is shown in Appendix \ref{ax:wifi}. The latent space is plotted with the variance of the mapping as the gray background in Fig. \ref{wifi:latent}. The 36  WiFi signal strength measurements are represented by the blue triangles. The training set (the labeled points) is marked by the red crosses. The dark color value is for high uncertainty. Since the mobile device moves in a loop closure, the latent points forms a closed loop in the latent space. The boundary of the implicit manifold $\partial M$ is shown in Fig. \ref{wifi:bound} which is defined by \eqref{eqn:bound}. The magnification factor is  plotted as the gray background in Fig.\ref{wifi:mag}. The dark color value is for high magnification factor. A sample path of BM in the implicit manifold is plotted in Fig. \ref{wifi:bm}. With the Neumann boundary condition the BM path can only exist within the boundary.

The coordinates of the mobile device are estimated by five different methods. The ground truth values are marked by different colors in the latent space in Fig. \ref{wifi:true}. If we ignore the interior structure of the manifold, a $\mathbb{R}^{30}$ GP using the squared exponential kernel with $\mathbb{R}^{30}$ Euclidean distance of the WiFi signals is applied. The predictive means of the $\mathbb{R}^{30}$ GP are plotted in Fig.\ref{wifi:hd}. It is clear the $\mathbb{R}^{30}$ GP prediction is very poor. The color coding of the prediction is different from the ground truth. The GPUM   predictive mean is shown in Fig. \ref{wifi:in-GP}. The overall pattern of the GPUM  prediction is similar to the ground truth. In the third case, the $\mathbb{R}^2$ Euclidean distance in the latent space is used with the squared exponential kernel to construct the $\mathbb{R}^2$ GP. The prediction results are shown in Fig.\ref{fig:wifi2d} in Appendix \ref{ax:wifi}. The Graph Laplacian based approaches such as GL-GP and GM-GP have also been applied. Since the number of the observations are relatively low ($n+v=36$), the graph based methods result in poor approximation of the implicit manifold. The GL-GP's performance is similar to the $\mathbb{R}^{30}$GP. The mean and standard deviation of the root mean square errors of the twenty testing sets are also calculated in Table \ref{tb:wifi} for all five models. GPUM  significantly outperforms the other four models and achieves the minimum mean RMSE.

\begin{table}[t]
\caption{\label{tb:wifi} Comparison of the root mean squared errors of five methods on WiFi signal data. Values in parentheses show the standard deviation. }
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
    & $\mathbb{R}^{30}$GP & $\mathbb{R}^2$GP & GPUM  &  GL-GP &  GM-GP \\
\hline
mean RMSE & 5.57(1.43)  & 4.83(1.83) & 4.11(0.88) & 5.6(1.15) & 6.04 (0.57) \\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\begin{figure}[t]
    \centering
    \subfigure[Ground truth]{\label{wifi:true}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_true.png}}
     \subfigure[ \small{ $\mathbb{R}^{30}$ GP prediction }]{\label{wifi:hd}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_30d.png}}
    \subfigure[GPUM   prediction]{\label{wifi:in-GP}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_mani.png}}

    \caption{\label{fig:wifiGP}
   {\footnotesize
    Comparison of GPUM  and Euclidean GPs in WiFi signal example. (a) The ground truth is plotted with color at the observed points in the latent space. (b) $\mathbb{R}^{30}$ GP prediction. $\mathbb{R}^{30}$ GP is constructed using the Euclidean distance of WiFi signals in $\mathbb{R}^{30}$. (c) GPUM  prediction.
    }
 }
\end{figure}

\section{Camera angle estimation from images} \label{sec:real}
In this section, we consider the problem of estimating the camera angle from images. We consider the setting in which an object is placed on a turntable and a set of images are taken at different angles with respect to the camera. We aim at recovering the camera angles associated with individual images by knowing the true camera angles for some of the images. We use the images from the COIL data set \citep{nene1996}. Here we consider 66 images of object-14 (a toy cat). The camera angles range from 15 degrees to 340 degrees.
The raw images are converted to grayscale and downscaled to $32 \times 32$.
The raw pixels of each image are flattened into a 1024 dimensional vector.
The dimension of the original image space is $p=1024$. Six image examples are given in Fig. \ref{coil:image}. We estimate a $q= 2$ dimensional latent space using Bayesian GPLVM. The two dimensional latent space is plotted with the variance of the mapping as the gray background in Fig. \ref{coil:var}. A dark background represents a high uncertainty in the corresponding region. The unlabeled data points are marked as blue triangles and the labeled data points as red stars. The overall shape of the latent points in Fig.\ref{coil:var} looks like a ring with a gap in the lower right. The magnification factor is plotted as the gray background in Fig. \ref{coil:mag}. A dark background represents a high magnification factor. The boundary of the implicit manifold $\partial \mathbb{M}$ is shown in Fig.\ref{coil:bound} in Appendix \ref{ax:coil}.
\begin{figure}[t]
    \centering
\subfigure[ COIL images ]{\label{coil:image} \includegraphics[width=0.31\textwidth,height=0.28\textwidth]{coilimage.png} }
\subfigure[ Latent space and variance of mapping]{ \label{coil:var} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{ncoilLatentVar.png} }
\subfigure[Latent space and magnification factor]{ \label{coil:mag} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{ncoilLatentMag.png} }
    \caption{ \label{fig:coil}
    \footnotesize{ (a) Examples of six COIL object images. (b) The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and red stars are labeled points. The background color represents the variance of the mapping. The dark color represents high uncertainty. (c) The same latent space visualization with the background color representing the magnification factor. The dark color represents a high magnification factor.
      }
  }
\end{figure}

The scaled camera angle $y$ is modelled as a function of the image.
 \begin{align*}
y_i = f(s_i) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^{1024} .
 \end{align*}
where $f$ is the unknown regression function, and $s_i$ is represented by a 1024 dimensional vector. n=13 images are randomly selected and used as the training set (labeled data). They are plotted as the red stars in the latent space in Fig.\ref{coil:var}. The remaining 53 images are used as the testing set (unlabeled data)  which are plotted as the blue triangles in the latent space in Fig. \ref{coil:var}. The true angle values are marked by different colors on the observed points in the latent space in Fig. \ref{coil:true}. Different methods have been applied to estimate the scaled camera angles. The $\mathbb{R}^{1024}$GP using the squared exponential kernel with $\mathbb{R}^{1024}$ Euclidean distance is applied to the image data first. The predictive means of the $\mathbb{R}^{1024}$GP are plotted in Fig.\ref{coil:hd}. The color coding of the prediction at the lower right of the plot is bright green which is different from the truth. Ignoring the boundary and the magnification factor in the latent space, the $\mathbb{R}^2$GP using the squared exponential kernel with $\mathbb{R}^2$ Euclidean distance in the latent space is applied to the image data. The predictive means are shown in  Fig.\ref{coil:2d} in Appendix \ref{ax:coil}. The predictive means of GPUM  are plotted in Fig.\ref{coil:in-GP}. The overall pattern of the GPUM  prediction is very similar to the true function. Since the boundary of the implicit manifold is defined in the lower right region of the latent space, the GPUM  does not smooth across the boundary and gives a better prediction. Ten training and testing
sets are generated from the random selection.
The mean and standard deviation of the root mean square errors are calculated for the ten testing sets for all methods in Table \ref{tb:coil}. The GPUM achieve the smallest mean RMSE. It is significantly better than  the GL-GP and GM-GP. The difference between the GPUM  and the $\mathbb{R}^{1024}$GP is not significant. The boxplot of the root mean squared errors are shown in Fig.\ref{coil:box} in Appendix \ref{ax:coil}.

 \begin{figure}[t]
    \centering
    \subfigure[True function]{\label{coil:true}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilTrue.png}}
     \subfigure[ \small{ $\mathbb{R}^{1024}$ GP prediction }]{\label{coil:hd}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilHd.png}}
    \subfigure[GPUM  on manifold prediction]{\label{coil:in-GP}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilIn-GP.png}}

    \caption{\label{fig:coilGP}
   { \footnotesize
     Comparison of GPUM  and Euclidean GPs in COIL images example. (a) The true function is plotted with color at the observed points in the latent space. (b) $\mathbb{R}^{1024}$ GP prediction. $\mathbb{R}^{1024}$ GP is constructed using the Euclidean distance of data points in $\mathbb{R}^{1024}$. (c) GPUM  prediction.
    }
 }
\end{figure}

\begin{table}[t]
\caption{\label{tb:coil}Comparison of the root mean squared errors of five methods on COIL images. Values in parentheses show the standard deviation.}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
&$\mathbb{R}^{1024}$GP &$\mathbb{R}^2$GP  & GPUM  &  GL-GP& GM-GP  \\
\hline
mean RMSE & 0.097(0.040) &0.094(0.042) & 0.060(0.038)  & 0.116(0.027) & 0.143(0.026) \\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\section{Discussion} \label{discuss}

Our work provides a novel framework for regression on implicit manifolds embedded in high dimensional point clouds. The geometry of the implicit manifold is learned by probabilistic latent variable models. This gives a distribution over a smoothly changing local metric at each point in the latent space.  The boundary of the implicit manifold is defined according to the uncertainty of the mapping. The expression for the expected local metric is derived and used to simulate BM sample paths on manifolds. We have proved the BM transition density estimation is coordinate independent in section \ref{sec:density}. This allows us to compare the transition density estimates from different parameterisations of the same manifold. The BM simulation using the B-GPLVM metric gives similar transition density estimation results as the BM simulation using the analytical metric. GPUM  is constructed by employing the equivalence relationship between heat kernels and the transition density of BM on manifolds. This allows the GPUM  to incorporate the intrinsic geometry of the implicit manifold for inference while respecting the interior constraints and boundary. The experiment results in section \ref{sec:simdata}, \ref{sec:wifi} and \ref{sec:real} indicate that GPUM  achieves significant improvements over Euclidean GPs and Graph Laplacian based GPs. Although the simulation of BM paths can be easily parallelised, Algorithm \ref{alg:bm} can still be computational expensive if the sample size is large. The number of the sample paths can be reduced by leveraging the idea of pseudo data, which has been widely studied in the literature of sparse GPs \citep{QuioneroCandelaRasmussen2005,niu2019}.

\section{Introduction}

Gaussian Processes (GPs) have been widely used as data efficient modelling approaches that produce good uncertainty estimates. They also power many decision making approaches such as Bayesian optimisation, multi-armed bandits and experiment design. Characteristics of the function prior such as differentiability, periodicity and symmetry can be easily controlled by constructing a specific covariance function. Most widely used GP covariance functions are defined over Euclidean space. However, in real world applications, data often lies on a manifold within the original space. Predictions are only valid on the manifold and should only be extrapolated from observations along the manifold. For example, traffic flows can only be measured over networks of roads, surface tension is only measured on the surface of a specific object and the joints of a robot arm can only be moved safely within a manifold of the joint space of all the joints. The traffic connectivity can also easily differ from road network due to road maintenance or traffic congestion. The inference of the manifold based on data brings the accuracy and flexibility.
Directly applying a GP with a covariance function defined over Euclidean space would not be ideal as the extrapolation would not respect the geometry of the manifold.

Previously, several GP on manifold methods \citep{niu2019,lin2019,borovitskiy2020} have been proposed under the assumption that the geometry of the manifold is known. \cite{niu2019} defined a heat kernel on a manifold and constructed an intrinsic Gaussian Process on the manifold. The intrinsic GP refers to a GP that employs the intrinsic Riemannian geometry of the manifold, including the boundary features and interior conditions. \cite{lin2019} proposed extrinsic framework for GP modeling on manifolds, which relies on embedding of the manifold into a Euclidean space and then constructing extrinsic kernels for GPs on their images. \cite{dunson2020diffusion} , \cite{borovitskiy2021} and \cite{bolin2022} focused on developing GPs on graphs and metric graphs formed by data observed on the manifold.

In this work, we study GP modelling on manifolds of which the geometry is unknown. This is a more realistic setting for real world applications because measuring the exact geometry of data manifolds is often impossible or overly expensive. In particular, we focus on the scenario where only a sparse set of data points can be collected. As shown in our experiments, Graph Laplacian based methods perform poorly in this scenario due to the poor graph approximation. In contrast, we estimate the probabilistic parameterisation of the unknown manifolds using probabilistic latent variable models and propose a practical and general intrinsic GP on unknown manifolds (GPUM) methodology. This is the major novel contribution of the paper. Specifically, we investigate the geometrical structure of the unknown manifolds using Riemannian geometry. The distributions of the metric tensor and the boundaries of resulting manifolds are estimated using the Bayesian Gaussian process latent variable models (B-GPLVM). Brownian Motion(BM) sample paths on the unknown manifold are simulated using the probabilistic metric and respecting the boundaries.
The covariance kernel on the unknown manifold is estimated by employing the equivalence relationship between the heat kernel and the transition density of BM on the unknown manifold. We prove this estimator is coordinate independent.
 Our method can incorporate the intrinsic geometry and the uncertainty of the unknown manifold for inference and respect the interior constraints. With numerical experiments on synthetic and real world datasets, we compared our method against Graph Laplacian based methods and GP regressions without manifold estimations. We showed that our method is significantly better than these other methods on all the datasets.

In the following sections, concepts of Riemannian geometry are introduced in Section \ref{sec:concepts}. The metric learning algorithms are explained in Section \ref{sec:metric}. We prove the BM paths on manifolds simulated using different metrics have the same transition density in Section \ref{sec:heatBM}. The heat kernel estimates derived from the analytical metric and different metric learning methods such as Gaussian Processes Latent Variable Model (GPLVM) and B-GPLVM are compared in Section \ref{sec:density}. Applications of GPUM  on a synthetic dataset on Swiss roll, high dimensional real datasets of WiFi signals \citep{ferris2007} and COIL images \citep{nene1996} are illustrated and compared to Graph Laplacian GPs \citep{dunson2020diffusion}, Graph Mat\'{e}rn GPs \citep{borovitskiy2021} and Euclidean GPs in Section \ref{sec:simdata}, \ref{sec:wifi} and \ref{sec:real}.

\section{Concepts of Riemannian Geometry and Theoretical Background}
\label{sec:concepts}

One way of representing a high dimensional dataset is to relate it to a lower dimensional set of latent variables through a set of mapping functions (potentially nonlinear). \cite{tosi2014metric,arvanitidis2019} investigated the geometrical structure of probabilistic generative dimensionality reduction models (or latent variable models) using Riemannian metrics and computed the geodesic distances on the manifolds learned from data. A manifold embedded in a high dimensional Euclidean space can be well described by a probabilistic mapping function and the corresponding latent space. If the dimension of the latent space is the same as the intrinsic dimension of the manifold, the latent space can be interpreted as the chart of the learned manifold. Intuitively, the chart provides a distorted view of the manifold.
An illustration is shown in Fig.\ref{fig:mapping}. The $x^1$ and $x^2$ coordinates in the chart represent the radius and width of the Swiss roll in $\mathbb{R}^3$. The blue triangles in the chart can be mapped to the black points in the embedded space (Swiss roll in $\mathbb{R}^3$)  through $\phi$. In other words, the Swiss roll can be parameterised by the chart. Measurements on the manifold can be computed in the chart locally, and integrated to provide global measures. This gives rise to the definition of a local inner product, known as a Riemannian metric tensor.

Let $\mathbb{M}$ be a $q$-dimensional complete and compact  Riemannian manifold with the Riemannian metric $\mathbf{g}$, and $\partial \mathbb{M}$ its boundary.  The Riemannian metric $\mathbf{g}$ can be represented as a symmetric, positive definite matrix-valued function, which defines a smoothly varying inner product in the tangent space of $\mathbb{M}$. Let $\mathcal J$ denote the Jacobian of $\phi$. We have
\begin{align} \label{eqn:deter_metric}
\mathbf{g} = \mathcal{J}^T \mathcal{J}, \ \ \  \mathcal{J}_{i,j} = \frac{\partial \phi^i }{ \partial  x^j}
\end{align}
The superscript indicates the $j_{th}$ dimension of the chart and the $i_{th}$ dimension of the observation space.

\begin{figure}[t]

    \subfigure[Chart of $\mathbb{M}$(Swiss roll)]{\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{map.pdf} \label{fig:mapping} }
    \subfigure[Regression on $\mathbb{M}$]{\includegraphics[width=0.2\textwidth,height=0.2\textwidth]{3dregSw.png} \label{fig:regexample} }
    \caption{
    \footnotesize{
    { An illustration of a parameterisation of the manifold $\mathbb{M}$(Swiss roll) estimated from a data cloud in (a). $\phi$ maps the chart into $\mathbb{M} \subset \mathbb{R}^3$. The blue triangles in the chart are mapped to the black dots on $\mathbb{M} \subset \mathbb{R}^3$. The Riemannian geometry can help to learn the regression function in (b) (the color indicates function values)
    .}}
    }
\end{figure}

 Moreover, based on its metric tensor, $\mathbb{M}$ has an associated Laplace-Beltrami operator, which is an intrinsically defined differential operator denoted $\Delta_s$. In local coordinates, the Laplacian-Beltrami operator i
\begin{align}
\label{eqn:delta_s}
\Delta_s f = \frac{1}{\sqrt{G}} \frac{\partial}{\partial x^j} \left (    \sqrt{G} \mathbf{g}^{ij} \frac{\partial f}{\partial x^i}  \right ),
\end{align}
where $G$ is the determinant of the metric $\mathbf{g}$, $\mathbf{g}^{ij}$ is the $(i,j)$ element of its inverse and $f$ is a smooth function on $\mathbb{M}$. Take the special case where $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$, $\mathbf{g}$ becomes an identity matrix. The Laplace-Beltrami operator $\Delta_s$ becomes the Laplace operator $\Delta$ (the sum of second partial derivatives).

Consider the heat equation on $\mathbb{M}$, given by

\begin{align*}
\frac{\partial}{\partial t}K_{heat}(s_0,s,t)=\frac{1}{2}\Delta_s K_{heat}(s_0,s,t),\qquad  s_0, s\in \mathbb{M},
\end{align*}
where $s \in \mathbb{M}$, $\Delta_s$ is the Laplacian-Beltrami operator on $\mathbb{M}$, and $t \in \mathbb{R}^{+}$ is the diffusion time. A heat kernel of $\mathbb{M}$ is a smooth function $K(s_0,s,t)$ on $ \mathbb{M} \times \mathbb{M} \times \mathbb{R}^{+}$ that satisfies the heat equation. It can be interpreted as the amount of heat that is transferred from $s_0$ to $s$ in time $t$ via diffusion.

The heat kernel satisfies the initial condition $\lim_{t\rightarrow 0}K_{heat}(s_0,s,t)=\delta(s_0,s)$ with $\delta$ the Dirac delta function. The heat kernel becomes unique when we impose a suitable condition along the boundary $\partial \mathbb{M}$, such as the Neumann boundary condition: $ \partial K / \partial \mathbf{n} = 0$ along $\partial \mathbb{M}$ where $\mathbf{n}$ denotes a normal vector of $\partial \mathbb{M}$. The Neumann boundary condition can be understood as allowing no heat transfer across the boundary. If $\mathbb{M}$ is a Euclidean space $\mathbb{R}^q$,  The heat kernel has a closed form corresponding to a time-varying Gaussian function:
\begin{align*}
K_{heat}(\mathbf{s}_0,\mathbf{s},t)
=\frac{1}{(2\pi t)^{q/2}}\,
  \exp\left\{-\frac{||\mathbf{s}_0-\mathbf{s}||^2}{ 2t }\right\}, \; \mathbf{s}\in \mathbb R^q.
\end{align*}
The diffussion time $t$ controls the rate of decay of the covariance. In the following, we will also write $K_{heat}(\mathbf{s}_0,\mathbf{s},t)$ as $K_{heat}^t(\mathbf{s}_0,\mathbf{s})$.

For arbitrary Riemannian manifold, the construction of the heat kernel associated with the Laplace-Beltrami operator is not a trivial task and belongs to the fields of partial differential equations and differential geometry \citep{chavel1984}. To circumvent solving the heat equation on manifolds, \cite{niu2019} estimates the heat kernel as the BM transition density by simulating BM paths on a known manifold. In the case of point clouds, both the metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$ are unknown. In this work, we will use probabilistic latent variable models to learn the distribution of the metric $\mathbf{g}$ and define the boundary $\partial \mathbb{M}$ based on the uncertainty quantification of the mapping. The Laplace-Beltrami operator in \eqref{eqn:delta_s} is the infinitesimal generator of BM on manifolds \citep{hsu1988}. The BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations in the Ito form \citep{hsu1988,hsu2008}:
\begin{align}
\label{eqnBMhsu}
dx_i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x_j} \left(  \mathbf{g}^{ij}G^{1/2} \right) dt + \left( \mathbf{g}^{-1/2} dB(t)\right)_i,
\end{align}
where $\mathbf{g}$ is the metric tensor of $\mathbb{M}$, $G$ is the determinant of $\mathbf{g}$ and $B(t)$ represents an independent BM in the Euclidean space.

\section{Intrinsic Gaussian Processes on unknown manifolds}
\label{sec:GPmani}

In this work, we focus on the following model,
\begin{align} \label{eqn:reg}
y_i &= f(s_i) +\epsilon_i, \ \ \  \epsilon_i \sim \mathcal{N}(0, \sigma_{\epsilon}^2)
\end{align}
where  $f:\mathbb{M} \rightarrow \mathbb{R}$ is an unknown regression function. $y_i \in \mathbb{R}$ is a response variable. $s_i= (s_{i}^1,\ldots,s_{i}^p) \in \mathbb{M} \subset \mathbb{R}^p$ is an observed predictor on a complete and compact Riemannian manifold $\mathbb{M}$ embedded in $\mathbb{R}^p$.
$\mathbb{M}$ can be parameterised by a $q$ dimensional local coordinate system (chart), and $q<p$. For example, let $\mathbb{M}$ be a unit sphere and we have a point $i$ on $\mathbb{M}$. $s_i \in \mathbb{R}^3$ is a three dimensional vector representing the Cartesian coordinates of $i$. The two dimensional chart of $\mathbb{M}$ is the longitude and latitude coordinate system. The point $i$ can also be represented as a two dimensional vector $x_i$ in the chart. There exists a mapping function $\phi$ from the Cartesian coordinates to the spherical coordinates. However in the case of point clouds where $s_i$ is observed in $\mathbb{R}^p$, the parameterisation $\phi$ of $\mathbb{M}$ is unknown and $x_i$ becomes a latent variable.
We propose to infer how the output $y$ varies with the input $s$, including predicting $y$ values at new locations not represented in the training set.

A GP prior can be assigned to $f$ with a covariance function.
The choice of
covariance kernel has a fundamental impact on the results. Most common choices of covariance kernels such as the squared exponential kernel and Mat\'{e}rn kernel depend critically on the Euclidean distance between $s_i$ and $s_j$ and ignore the intrinsic geometry of $\mathbb M$. By contrast, the heat kernel depends only on the Riemmanian metric and the intrinsic geometry of $\mathbb{M}$. It provides a natural generalisation of the RBF kernel on manifolds. The heat kernel represents the diffusion of the heat on a Riemannian manifold. However it is analytically intractable to directly evaluate \citep{hsu1988}. %
\cite{niu2019} proposed a computational framework to estimate the heat kernel on the manifold in which the analytical parameterisation is known.
However, data represented as point clouds is often high dimensional and concentrated around some unknown lower dimensional structures, \cite{niu2019} is not applicable due to the lack of analytical manifold parameterisation.
In this paper, we address this problem by using the probabilistic generative dimension reduction models to learn the geometry of the implicit manifold and define the boundary. We propose to construct GPUM  by using the heat kernel of the implicit manifold as the covariance kernel. With the help of Riemannian geometry, the Brownian motion sample paths can be simulated on $\mathbb{M}$ and the heat kernel can be estimated as the transition density of the BM. With this construction we can learn the regression function on $\mathbb{M}$ as in Fig.\ref{fig:regexample}.

Let $\mathcal{D} =  \{ (s_i, y_i), i = 1, . . . , n \} $ be the data, with $n \geq 1$ the number of labeled observations, $s_i \in \mathbb{M} \subset \mathbb{R}^p$ is the predictor and $y_i \in \mathbb{R}$ is the corresponding response. Suppose we are also given an unlabeled dataset $\mathcal{V} = \{ s_i, i=n+1 . . . ,n+v \} $ where $v \geq 1$. Consider the regression model in \eqref{eqn:reg}, we would like to make inferences about $f$ with the labeled dataset $\mathcal{D}$ and predict $y$ values for the unlabeled dataset $\mathcal{V} $. Under an GPUM  prior for the unknown regression function as $f \sim GP(0,K_{heat}^t(.,.) )$, we have
\begin{align*}
p( \mathbf{f} | s_1, s_2, ..., s_n) \sim \mathcal{N} (0, \Sigma_{\mathbf{f} \mathbf{f}}),
\end{align*}
where $\mathbf{f}$ is the discretization of $f$ over $s_1,s_2,\cdot\cdot\cdot,s_n$ so that $f_i = f(s_i)$. $\Sigma_{\mathbf{f} \mathbf{f}} \in \mathbb{R}^{n\times n}$ is the covariance matrix induced from the heat kernel. The $(i,j)$ entry of $\Sigma_{\mathbf{f}\f}$ corresponds to $\Sigma_{\mathbf{f}\f _{i,j}} = \sigma_h^2 K_{heat}^t(s_i,s_j)$. We introduce the rescaling hyperparameter $\sigma_h^2$ to add extra flexibility to the heat kernel.

\section{Related work}
One of the paramount challenges in developing GP models on manifolds is the difficulty in specifying the covariance structure via constructing valid covariance kernels on manifolds. One might hope to achieve this by replacing Euclidean norms in the squared exponential kernel with geodesic distances. However \cite{feragen2015} proved this is not generally a well-defined kernel. \cite{extrinsicGP} proposed extrinsic Gaussian Processes on manifolds by embedding the manifolds onto a higher dimensional Euclidean space. The squared exponential kernel is applied on the images of the manifold after embedding. However, such embeddings are not easy to obtain and only available for certain manifolds when the geometry is known.
\cite{yang2016} proposed a model bypassing the need to estimate the manifold, and can be implemented using standard algorithms for posterior computation in GPs.
They show that by imposing a GP prior on the regression function with a covariance kernel defined directly on the ambient space (the embedding of the manifold in a high dimensional Euclidean space), the posterior distribution yields a posterior contraction rate depending on the intrinsic dimension of the manifold. They assume that the unknown lower dimensional space where the predictors center around are a class of submanifolds of Euclidean space.  They focus on compact manifolds without boundary.

\cite{niu2019} proposed to use heat kernels to construct the intrinsic Gaussian Processes on complex constrained domains. Heat kernel can be interpreted as the transition density of Brownian Motion on the manifold \citep{hsu1988}. It is estimated by simulating BM paths on manifolds in which the analytical parameterisation is known in \cite{niu2019}.

Alternatively if the eigen-paires of the Laplacian-Beltrami operator of the manifold are available, \cite{borovitskiy2020} approximated the heat kernel with the sum of finite-many eigen-paires of the Laplacian-Beltrami. Both \cite{borovitskiy2020} and \cite{niu2019} are only applicable when the geometry of the manifold is known and the dimensions of the observation(or embedding) space are low such as $\mathbb{M} \subset \mathbb{R}^2$ or $\mathbb{M} \subset \mathbb{R}^3$ ($p= 2$ or $3$).

Most recent research in \cite{dunson2020diffusion} tackled this problem by approximating the heat kernel kernel of a compact Riemannian manifold with finitely-many eigenpairs of the Graph Laplacian using the labeled and unlabeled predictor values. We refer the Gaussian processes constructed by this approximation as the Graph Laplacian Gaussian processes (GL-GPs). Let $\Delta$ be the Laplace-Beltrami operator of a manifold and $ \{ \lambda_i \}_{i=0}^{\infty}$ be the spectrum or eigenvalues of $-\Delta$. Denote $\varphi_i$ the corresponding eigenfunction, for each $i \in N$, we have $\Delta \varphi_i = - \lambda_i \varphi_i$. If the eigen-decomposition of $\Delta$ is known, the corresponding heat kernel of the manifold has the following expression: $K_{heat}(s,s',t) = \sum_{i=0}^{\infty} e^{-\lambda_i t} \varphi_i(s) \varphi_i(s')  $. If the geometry of the manifold and the corresponding $\Delta$ are unknown, $\Delta$ is approximated  by the Graph Laplacian matrix $\bm{L}$ in \cite{dunson2020diffusion}.  $\bm{L}$ is constructed from a point cloud whose adjacency matrix is computed by using a Gaussian function with pair-wised Euclidean distance. The heat kernel is approximated as the summation of finite eigenpairs of $\bm{L}$,  $\sum_{i=0}^{n_G} e^{-\mu_i t} v_i v_i^T $, where $\mu_i$ and $v_i$ are the $i_{th}$ normalised eigenvalue and eigenvector of $\bm{L}$. The implementation of the GL-GP is provided in  Appendix \ref{ax:GL_kern}.

\cite{borovitskiy2021} leveraged the stochastic partial differential equation characterization of Mat\'{e}rn kernel to study their analog for undirected graphs and developed the Graph Mat\'{e}rn Gaussian processes (GM-GPs). The Graph Mat\'{e}rn kernel is constructed with the sum of finitely-many eigenpairs of the normalised Graph Laplacian $\bm{L}$.
$\bm{L}$ is computed from the adjacency matrix of a predefined graph. We implemented the GM-GP by following the instructions in the Github repository in \cite{borovitskiy2021}. The key component of GL-GPs and GM-GPs is the Graph Laplacian. If the true graph connections are not known, the graph constructed based on local distances such as Delaunay triangulation can be error-prone when observations are sparse \citep{hjelle2006}. The graph based methods often result in poor approximation of the manifold when the number of the observations are low. GL-GPs and GM-GPs are applied to the simulation studies and real datasets and compared to GPUM in the later sections.

There are still some critical gaps in current practices. In particular, the lack of robust methods for carrying out intrinsic statistical inference and  effective models for regressions with manifold-valued data embedded in a point cloud.
In this work we focus on learning the manifold structure using probabilistic dimension reduction methods such as Bayesian GPLVM and constructing the GPUM by using the heat kernel of the learned manifold. Other related methods such as Auto-encoders \citep{kramer1991} and Variational Auto-encoders (VAE)\citep{kingma2019}, can also be considered in this two stage approach. Auto-encoders provide a neural network based framework for learning deep latent variable models.  However, the resulting mapping function is deterministic and the model does not have a built-in quantification of its uncertainty.
There is lack of uncertainty quantification in the learned manifold. Alternatively, VAE
address this concern directly through an explicit likelihood model and a variational approximation of the representation posterior. It learns a generative model by specifying a likelihood of observations conditioned on latent variables and a prior over the latent variables.  Both the likelihood and the variational distributions have parameters predicted by
neural networks that act similarly to the encoder–decoder pair of the classic autoencoder.

\section{Learning metrics and boundaries}\label{sec:metric}
Let $\mathcal S = \{ s_i |  i=1,\cdot \cdot \cdot, n+v \}$,  $s_i \in \mathbb{R}^p$, include all predictors in the labeled dataset $\mathcal{D}$ and unlabeled dataset $\mathcal{V}$. We have $\mathcal{S} \in \mathbb{R}^{ (n+v) \times p }$. Suppose we perform probabilistic nonlinear dimensionality reduction by defining a latent variable model that introduces a set of latent (unobserved) variables $\mathcal{X}= \{x_i | i= 1,\cdot \cdot \cdot, n+v \}$, $x_i \in \mathbb{R}^q$, $\mathcal{X} \in \mathbb{R}^{(n+v)\times q}$ and $q < p$. $\mathcal X$ is related to $\mathcal S$ which is observed in a higher dimensional space. A prior distribution is placed on the latent space which induces a distribution over $\mathcal S$ under the assumption of the probabilistic mapping
\begin{align} \label{eqn:mapmodel}
s_{i}^j = \phi^j(x_i) + e_{i}^j,
\end{align}
where $x_i \in \mathbb{R}^q$ is the latent point associated with the $i_{th}$ observation $s_i \in \mathbb{R}^p$ in the original observation space and $i\leq n$. $j$ is the index of the features (dimensions) of $s$ in the observation space and $j\leq p$. $e_{i}^j$ is a Gaussian distributed noise term,   $e_i^{j} \sim \mathcal N(0,\beta^2)$.
An illustration is given in Fig. \ref{fig:mapping} where $q=2$ and $p=3$. If the mapping $\phi$ is linear and the prior $p( \mathcal X)$ is Gaussian, this model is known as probabilistic principal component analysis \citep{tipping1998}. In this work, we do not restrict to this linear assumption and consider some nonlinear dimensionality reduction methods such as GPLVM and B-GPLVM.

When the $\phi$ function in Fig.\ref{fig:mapping} is differentiable, it can be interpreted as the mapping between the latent space and the manifold $\mathbb{M}$. If the dimensions of $\mathbb{M}$ are known, by letting $q$ equal the dimensions of $\mathbb{M}$, the latent space can be interpreted as the chart of $\mathbb{M}$. If the dimensions of $\mathbb{M}$ are unknown, $q$ is estimated by using the so called Automatic Relevance Determination (ARD) in GPLVM and B-GPLVM \citep{zwiessele2017}. ARD allows assigning scaling parameters for each dimension. These scaling parameters can be incorporated to kernels such as the RBF kernel as the inverse of the squared lengthscales.

In general
a manifold may need more than one chart to be parameterised. In this work, we focus on examples of single chart. But our method can be extended to learn multiple charts from multiple datasets.
The Riemannian metric of the given model can be computed as in \eqref{eqn:deter_metric}. In the case of probabilistic LVMs, we place a Gaussian prior over the mapping function $\phi(x)|x$. The conditional probability over the Jacobian also follows a Gaussian distribution, this naturally induces a distribution over the metric tensor $\mathbf{g}$. We denote the distribution of the Jacobian as $\mathbf{J}$ given the set $\mathcal{S}$ and the mapping $\phi$ in \eqref{eqn:mapmodel}. Assuming independent rows of $\mathbf J$ \citep{lawrence2005,titsias2010},
\begin{align}
p(\mathbf J | \mathcal{S}, \Phi) = \prod_{j=1}^p \mathcal{N}( \mu_{\mathbf{J}}^j, \Sigma_{\mathbf J}  ).
\end{align}
This independent row assumption is for the dimensions in the original
observational space. The dimensions in the learned latent space are
not independent. This assumption can be relaxed by using a multi-output GP which is more compuationaly expensive \citep{alvarez2011}. The expressions of the mean $\mu_{\mathbf J}$ and variance $\Sigma_{\mathbf J} $ of the Jacobian are model specific and given in section \ref{gplvm} and \ref{bgplvm}. The resulting metric $\mathbf{g}$ follows a non-central Wishart distribution \citep{anderson1946}
\begin{align} \label{eqn:metric}
\mathbf{g} \sim \mathcal{W}_q \left(  p, \Sigma_{\mathbf{J}}, \E(\mathbf{J}^T)  \E(\mathbf{J})  \right).
\end{align}
From this distribution, the expected metric tensor can be computed as
\begin{align} \label{eqn:exmetric}
\mathcal{G} = \E( \mathbf{g} )  = \E( \mathbf{J}^T) \E( \mathbf{J} ) + p \Sigma_{\mathbf{J}}.
\end{align}
We denote the expectation by $\mathcal{G}$. Note that the variance term $\Sigma_{\mathbf{J}}$ is included in $\mathcal{G}$. It implies that the metric tensor expands as the uncertainty over the mapping increases.
Hence the BM simulation steps in the SDE in (\ref{eqnBM}) will travel `slower' in the region of the latent space where the uncertainty is high.

We also need the gradient of the expected metric to simulate BM as in section \ref{trandens}.
\begin{align}
\frac{\partial \mathcal{G}}{\partial x^l} = \frac{ \partial \E[ \mathbf{g}]}{\partial x^l} = \frac{\partial \E[\mathbf{J}^T] } { \partial x^l }  \E[\mathbf{J}] + \E[\mathbf{J}^T] \frac{ \partial  \E[\mathbf{J}]  }{\partial x^l}  + p \frac{ \partial \Sigma_\mathbf{J} }{ \partial x^l }
\end{align}
The estimates of the mapping  $\phi(x)$ become highly unreliable when $x$ is far from the data points. The corresponding metric tensor estimates and BM simulations can also be ill defined and violate the manifold geometry. To avoid these poorly estimated region due to lack of data, the boundary of the learned manifold can be defined by $\mathbb{V} ar( \phi(x) | x)$, the variance of the mapping at $x$. $\mathbb{V} ar( \phi(x) | x)$ is also a smooth function in B-GPLVM. In most cases, the $\partial \mathbb{M}$ defined here is also a $q-1$ smooth manifold. Any $x$ outside of the boundary has $\mathbb{V} ar( \phi(x) | x ) > \alpha$.
\begin{align}\label{eqn:bound}
\partial \mathbb{M} = \{  x \in \mathbb{R}^q \  \vert \ \  \mathbb{V} ar( \phi(x) | x) = \alpha   \}.
\end{align}
where $\alpha$ is defined as the maximum variance of the mapping at the shifted latent points, $\max \left( \mathbb{V} ar\left( \phi(h) | h \right) \right)$. The shifted latent points are randomly sampled from the set $\mathbb{H}=\{ h |  \| h -x_i\| = \delta_{\mathcal{X}}, x_i\in \mathcal{X}  \}$.  $\delta_{\mathcal{X}}$ can be chosen according to the maximum distance of two neighbouring latent points. In practice, we create $h$ by moving all data points with $\delta_{\mathcal{X}}$ in random directions and let $\alpha$ equal to the maximum variance of the mapping at these relocated points. The samples are chosen to make sure that $\alpha > \max \left( \mathbb{V} ar\left( \phi(x_i) | x_i \right) \right)$, for all $x_i\in \mathcal{X}$.

\subsection{ GPLVM metric}\label{gplvm}
Gaussian Process Latent Variables Model \citep[GPLVM; ][]{lawrence2005} is a nonlinear probabilistic generative model. In this section, we will derive the distribution of the metric from GPLVM.
A sample from GPLVM defines a generative mapping from $x\in \mathbb{R}^q$ in the latent space  to $s \in \mathbb{M} \subset \mathbb{R}^p$ in the observation space.
GPs define a prior over the mapping $\phi$ in \eqref{eqn:mapmodel}. A zero mean prior is used as a default choice. If the domain
knowledge of where the prior should be centred is available, such as the embedding of $\mathbb{R}^q$ into $\mathbb{R}^p$, it can also be encoded into
the mean function. The individual dimensions of the $p$ dimensional observation space are modeled independently in the GP prior sharing the same hyperparameters. Given the construction outlined above, the probability of the observed data $\mathcal{S} =\{ s_i^j| i\in \{1,\cdots,n+v \},j\in \{1,\cdots,p\}  \}$ conditioned on all latent variables $\mathcal{X} = \{ x_i^j| i\in \{1,\cdots,n+v \},j\in \{1,\cdots,q\}   \}$ is written as follows:
\begin{align}\label{eqn:gplvmJoint}
p(\mathcal{S}, \Phi |\mathcal{X} , \beta) =   p(\mathcal{S} | \Phi ,\beta ) p(\Phi | \mathcal{X})   = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j, \beta ) p( \bm{\phi}_{:}^j | \mathcal{X}),
\end{align}
 where $\Phi = \{ \bm{\phi}_i^j | i\in \{ 1,\cdot\cdot\cdot, n+v \} , j\in \{ 1,\cdot\cdot\cdot, p \} \} $, $\Phi \in \mathbb{R}^{(n+v)\times p}$ and $ \bm{\phi}_i^j = \phi(x_i)^j$.  $\mathbf{s}_{:}^j$ represents the $j_{th}$ dimension of all points in $\mathcal S$ and $\mathcal S \in \mathbb{R}^{ (n+v)\times p}$. The likelihood $p( \mathcal S | \mathcal{X} )$ is computed by marginalising out $\Phi$ and optimising the latent variables $\mathcal{X}$
\begin{align}
p(\mathcal{S} | \mathcal{X} ) = \prod_{j=1}^p \mathcal N( \mathbf{s}_{:}^j, K_{\mathcal{X} \mathcal{X}} + \beta^2 I ).
\end{align}
$K_{\mathcal{X} \mathcal{X}} $ is the $(n+v)\times (n+v)$ covariance matrix defined by the squared exponential kernel.
\cite{lawrence2005} estimated all the latent variables $\mathcal X$ and the kernel hyper-parameters of GPLVM with the maximum likelihood estimate. If the covariance kernel is differentiable, the derivative of a GP is again a GP \citep{Rasmussen2004}. This property allows us to compute the derivative of GP. The Jacobian $\mathbf{J}$ of the GPLVM mapping can be computed as the partial derivative $\frac{\partial \phi_*}{ \partial x^l}$ with respect to the $l_{th}$ dimension for any point $x_*$  in the latent space,
\begin{align} \label{eqn:Jacobian}
\mathbf{J}^{T} = \frac{\partial \phi_*}{ \partial x} = \begin{bmatrix}
 \frac{ \partial \phi(x_*)^1 }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^1}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^1}   \\
  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot  &  \cdot \cdot \cdot \\
  \frac{ \partial \phi(x_*)^1 }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^j }{\partial x^q}   &   \cdot \cdot \cdot  &    \frac{ \partial \phi(x_*)^p }{\partial x^q}
\end{bmatrix},
\end{align}
where $\frac{\partial \phi_*}{ \partial x}$ is a $q \times p$ matrix. Considering the independence across the dimensions of the observation space, the joint distribution of the $j_{th}$ dimension of the mapping $\phi$ and the $j_{th}$ column of the Jacobian can be written as
\begin{align}
\begin{bmatrix}
 \phi( \mathcal{X})^j\\
 \frac{ \partial \phi(x_*)^j }{\partial x}
\end{bmatrix}, \
\sim \mathcal N \left( 0, \begin{bmatrix}
K_{\mathcal{X},\mathcal{X}} & \partial K_{\mathcal{X},*} \\
 \partial K_{\mathcal{X},*}^T & \partial^2 K_{*,*}
\end{bmatrix}  \right ).
\end{align}
The expressions of $K_{\mathcal{X},\mathcal{X}}$ ,  $\partial K_{\mathcal{X},*}$ , and  $\partial^2 K_{*,*} $ are given in Appendix \ref{ax:GPLVM}. GPLVM provides an explicit mapping from the latent space to the observation space. This mapping defines the support of the observed data $\mathcal{S}$ as a $q$-dimensional manifold embedded in $\mathbb{R}^p$. The distribution of the Jacobian of GPLVM is the product of $p$ independent Gaussian distributions (one for each dimension of the observation space). For a point $x_*$ in the latent space, the  distribution of the Jacobian takes the form
\begin{align}
p(  \mathbf{J} | \mathcal{X}, \mathcal{S}) &= \prod^p_{j=1} \mathcal N( \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \partial K^T_{\mathcal{X},*} K^{-1}_{\mathcal{X},\mathcal{X}} \mathbf{s}_{:}^j  \ , \partial^2K_{*,*} - \partial K^T_{\mathcal{X},*} K_{\mathcal{X},\mathcal{X}}^{-1} \partial K_{\mathcal{X},*}  ). \nonumber
\end{align}
From this distribution, the expected metric tensor can be computed as in \eqref{eqn:exmetric}. The boundary $\partial \mathbb{M}$ is defined by $\mathbb{V} ar( \phi(x_*) | x_*)$, the variance of mapping in \eqref{eqn:bound}. For any $x_*$ in the latent space
\begin{align}\label{eqn:varb}
\mathbb{V} ar(\phi(x_*) | x_*) = K_{*,*} - K_{*,\mathcal{X}}K_{\mathcal{X},\mathcal{X}}^{-1} K_{\mathcal{X},*}.
\end{align}

\subsection{Bayesian GPLVM metric }\label{bgplvm}
The maximum likelihood estimation of the latent inputs  $\mathcal{X}$ in GPLVM often leads to overfitting due to its high dimensionality \citep{damianou2016}.
This overfitting can be avoided by applying a Bayesian treatment to the latent inputs. By introducing a prior distribution to the latent inputs, the marginal likelihood takes the form
\begin{align*}
p(\mathcal S)=\int p(\mathcal{S} \mid \mathcal{X}) p(\mathcal{X}) d \mathcal{X}.
\end{align*}
The  integral is intractable as the inputs $\mathcal{X}$ to the latent variable model go through a non-linear calculation in the inverse of the covariance matrix. Bayesian Gaussian Process Latent Variable Model \citep[B-GPLVM;][]{titsias2010} introduces a variational inference framework for training the latent variable model. It variationally integrates out the input variables and computes a lower bound on the exact marginal likelihood of the nonlinear latent variable model. The maximization of the variational lower bound provides a Bayesian training procedure that is robust to overfitting.

The key to the tractable variational Bayes approach is the application of variational inference to an augmented GP formulation, known as sparse GP, where the GP prior on $\phi$ is augmented to include auxiliary variables.
More specifically, we expand the conditional probabilistic model in \eqref{eqn:gplvmJoint} by including $m$ extra samples (inducing points) to the GP latent mapping
e.g. $u_i = \phi(x_{ui}) \in \mathbb{R}^p$ is such a sample \citep{lawrence2007}.
These inducing points are denoted by $\mathcal{U} = \{u_i^j |  i \in \{ 1,\cdot\cdot\cdot,m \},  j \in \{  1,\cdot\cdot\cdot,p \}    \}$, $\mathcal{U} \in \mathbb{R}^{m\times p}$ and constitute latent function evaluations at a set of pseudo-inputs $\mathcal{X}_u = \{ x_{ui}^j | i \in \{1,\ldots,m \} , j \in \{ 1,\ldots,q\} \} $, $\mathcal{X}_u \in \mathbb{R}^{m \times q}$.
The inducing inputs $\mathcal{X}_u$ are variational parameters.
$\Phi$ is defined as in \eqref{eqn:gplvmJoint}. The augmented joint probability and the marginal likelihood can be written as \eqref{eqn:full-marginal}
\begin{align}\label{eqn:full-marginal}
p(\mathcal{S} , \Phi , \mathcal{U} ,  \mathcal{X} | \beta) &=   p(\mathcal{S} | \Phi ) p(\Phi | \mathcal{U} ,\mathcal{X}) p(\mathcal {U}) p(\mathcal{X})  = \prod_{j=1}^p p(\mathbf{s}_{:}^j | \bm{\phi}_{:}^j) p( \bm{\phi}_{:}^j | \mathbf{u}_{:}^j, \mathcal{X}) p(\mathbf{u}_{:}^j) p(\mathcal{X}), \\
p(\mathcal{S}) &= \int  \int  \int  p(\mathcal{S}, \Phi, \mathcal{U}, \mathcal{X}) \ d \mathcal{U}     \  d\Phi \ d \mathcal{X}.
\end{align}
With variational inference, $\log \left(p(\mathcal{S} )\right)$ can be lower bounded by applying Jensen’s inequality \citep{damianou2016}.
The resulting lower bound can be computed analytically. The full technical details of the lower bound are given in \eqref{apxeqn:vlowbou} in Appendix \ref{ax:BGPLVM}. The Jacobian has the same shape as in \eqref{eqn:Jacobian}. For a point $x_*$ in the latent space, the distribution of the Jacobian for B-GPLVM takes the form:
\begin{align}
p(  \mathbf{J} | \mathcal{X},\mathcal{U},\mathcal{S}) &= \prod^p_{j=1} \mathcal N( \  \mu_{\mathbf{J}}^j , \Sigma_{\mathbf{J}} ) \\
   &=  \prod^p_{j=1} \mathcal N( \  \partial K^T_{\mathcal{X}_u,*} K^{-1}_{\mathcal{X}_u,\mathcal{X}_u} \mu_{qu}^j , \  \partial^2K_{*,*} - \partial K^T_{\mathcal{X}_u,*} \Lambda \partial K_{\mathcal{X}_u,*}  )  \nonumber   \\
 \Lambda &= K_{\mathcal{X}_u\mathcal{X}_u}^{-1} - K_{\mathcal{X}_u\mathcal{X}_u}^{-1} \Sigma_{qu} K_{\mathcal{X}_u\mathcal{X}_u}^{-1} \nonumber  \\
\mathbb{V} ar(\phi(x_*) | x_*) &= K_{*,*} - K_{*,\mathcal{X}_u} \Lambda K_{\mathcal{X}_u,*}\label{eqn:varb}
\end{align}
where $\mu_{qu}^j$ and $\Sigma_{qu}$ are the mean and variance of the variational distribution of the inducing points $\mathcal{U}$. The expressions of $\mu_{qu}$ and $\Sigma_{qu} $ are given in \eqref{apxeqn:vquqsig} in Appendix \ref{ax:BGPLVM}. The expectation of the metric tensor can be computed as in \eqref{eqn:exmetric}. The boundary $\partial \mathbb{M}$ of the implicit manifold can be defined by computing the variance of the mapping as in \eqref{eqn:varb}.

\section{Estimate heat kernel as BM transition density} \label{sec:heatBM}

 In this section, we will estimate the BM transition density by simulating BM sample paths on the implicit manifold $\mathbb{M}$.

\subsection{ Simulating Brownian motion on implicit manifolds }
\label{trandens}
From section \ref{sec:metric}, we learned the probabilistic parameterisation $\phi$ of $\mathbb{M}$, the distribution of the associated metric tensor $\mathbf{g}$ and the boundary $\partial \mathbb{M}$. In order to estimate the BM transition density on $\mathbb{M}$, we first need to simulate BM trajectories. Simulating the sample paths of BM on $\mathbb{M} \subset \mathbb{R}^p$ is equivalent to simulating the stochastic processes in the latent space (or chart) in $\mathbb{R}^q$, $q<p$. An illustrated example is shown in Fig. \ref{fig:BM}.

\begin{figure}[t]
    \centering
\includegraphics[width=0.7\textwidth,height=0.3\textwidth]{bmlatent1.pdf}
    \caption{ \label{fig:BM}
    \footnotesize{
    {A BM sample path (blue line, right panel) on $\mathbb{M}$ (Swiss roll in $\mathbb{R}^3$) and its equivalent stochastic process (purple line, left panel) in the chart in $\mathbb{R}^2$. $\phi: \mathbb{R}^2 \rightarrow \mathbb{M} \subset \mathbb{R}^3$ is a parametrization of $\mathbb{M}$. The red dot is the starting location of the BM trajectory. The horizontal axis of the latent space represents the radius of Swiss roll. The vertical axis is for the width. The gray color in the latent space indicates the magnification factor. When the latent space is mapped to $\mathbb{M} \subset \mathbb{R}^3$, the darker region will be stretched more.}}
    }
\end{figure}
BM on a Riemannian manifold in a local coordinate system is given as a system of stochastic differential equations (SDE) in Ito form \citep{hsu1988}. We use the expected metric $\mathcal{G} = \E[\mathbf{g}]$ to construct the SDEs
\begin{align}
\label{eqnBM}
dx^i(t) = \frac{1}{2}G^{-1/2} \sum^{q}_{j=1}\frac{\partial}{\partial x^j} \left(  \mathcal{G}^{ij}G^{1/2} \right) dt + \left( \mathcal{G}^{-1/2} dB(t)\right)_i
\end{align}
where $x^i$ represents the $i_{th}$ dimension of the latent space (chart). $\mathcal{G}$ is defined in eqn (\ref{eqn:exmetric}),
$G$ is the determinant of $\mathcal{G}$ and $B(t)$ represents an independent BM in Euclidean space. The discrete form of \eqref{eqnBM} is derived in \eqref{disBM}.
\begin{align}
\label{disBM}
x^i(t) &= x^i(t-1) + \frac{1}{2} \sum^{q}_{j=1} \left(  -\mathcal{G}^{-1} \frac{\partial \mathcal{G}}{\partial x^j} \mathcal{G}^{-1} \right)_{ij} \Delta t + \frac{1}{4} \sum_{j=1}^q(\mathcal{G}^{-1})_{ij} tr\left(\mathcal{G}^{-1}\frac{\partial \mathcal{G}}{\partial x^j}\right) \Delta t + \left( \mathcal{G}^{-1/2} dB(t)\right)_i \nonumber \\
         &= \mu( x^i(t-1),\Delta t)_i + \left( \sqrt{\Delta t} \mathcal{G}^{-1/2} z^q\right)_i
\end{align}
where $\Delta t$ is the diffusion time of the BM simulation step and $z^q$ is a $q$-dimensional standard normal random variable. The discrete form of the SDE also defines the proposal mechanism of the BM
\begin{align}
\label{eqn:disProp}
q\left(x(t)|x(t-1)\right) = \mathbb N \left( x(t) |  \mu \left( x\left(t-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right ).
\end{align}

\begin{theorem} \label{theorem:sde}
The stochastic process defined in \eqref{eqnBM} is coordinate independent. With a given $\delta t$, simulations in any choice of local coordinates (or metric) as above are equivalent to the same step in $\mathbb{M}$.
\end{theorem}

The proof of Theorem \ref{theorem:sde} is given in Appendix \ref{ax:independent}. This implies that the BM sample paths simulated from different choices of metric $\mathcal{G}$ should have the same properties. The BM steps are sampled from the proposal distribution in \eqref{eqn:disProp}, which is defined by the metric tensor. The boundary $\partial \mathbb{M}$ of the manifold is also quantified by the uncertainty of the mapping. We apply the Neumann boundary condition as in Section \ref{sec:GPmani}. As a result the simulated sample paths stay within the boundary. An example of BM trajectory on Swiss roll is given in Fig. \ref{fig:BM}. The latent space and the associated metrics are learned from B-GPLVM. The gray color in the latent space indicates the square root of the determinant of the metric. It is also known as the magnification factor \citep{bishop1997,zwiessele2017,tosi2014}.
\begin{align*}
\mathcal{MF} = \sqrt{det(\mathcal{G}) }
\end{align*}
The geometric interpretation of the magnification factor is how much a small piece of the latent space in $\mathbb{R}^q$ will be stretched or compressed when it is mapped to $\mathbb{M} \subset \mathbb{R}^p$. For example in Fig.\ref{fig:BM} left panel, the horizontal axis of the latent space can be interpreted as the scaled radius of the Swiss roll, the vertical axis as the width of the Swiss roll. When the radius is bigger, the corresponding area in the latent space is darker and the magnification factor is larger. When the latent space is mapped back to the manifold, the darker region will be stretched more. The purple trajectory of the stochastic process in the latent space (Fig. \ref{fig:BM} left panel) is denser in the darker area and more spread out in the bright area. The stochastic process travels with smaller steps when the magnification factor is large and vice versa.
Note that, due to the BM being coordinate independent as in Theorem \ref{theorem:sde}, the BM trajectory is evenly spread out in the manifold of the Swiss roll in $\mathbb{R}^3$ (see Fig. \ref{fig:BM} right panel).

\subsection{Estimate the transition density of BM} \label{sec:density}
Considering the BM $\{ \mathbf{S}(t) | t>0 \}$ on $\mathbb{M} \subset \mathbb{R}^p$. The BM starts from $\mathbf{S}(0)= s_0$ at time $0$. We simulate $N_{BM}$ sample paths. Given a point $s \in \mathbb{M}$, we define a small neighbourhood of $s$ as $\mathbf{A}_s \in \mathbb{M}$. For any $t > 0$, the probability of $\mathbf{S}(t)$ reaching $\mathbf{A}_s$ at time t, $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$, can be approximated by
\begin{align}
\label{eqn:BMdensity}
p( \mathbf{S}(t)\in \mathbf{A}_s\,| \ \mathbf{S}(0)=s_0)  \approx  \frac{N_{A_s}}{N_{BM}}
\end{align}
where $N_{A_s}$ is the number of sample paths reaching $\mathbf{A}_s$ at time $t$. An illustrative diagram is shown in Fig. \ref{fig:BM-density}. This BM transition probability is defined as the integral of the BM transition density over $\mathbf{A}_s$ which is also the heat kernel $K_{heat}^t$.
Since we do not have the analytical expression of the transition probability,  the transition density cannot be derived by taking the derivative of $p( \mathbf{S}(t)\in \mathbf{A}_s | \mathbf{S}(0)=s_0 )$. Instead, $K_{heat}^t$ can be numerically approximated as $ \hat{K}_{heat}^t $
\begin{align}
\label{k:heat}
K_{heat}^t(s_0,s) \approx  \frac{p( \mathbf{S}(t) \in \mathbf{A}_s \ | \ \mathbf{S}(0)=s_0 )}{V(\mathbf{A}_s)} \approx \frac{1}{V(\mathbf{A}_s)}\cdot \frac{N_{A_s}}{N_{BM}} = \hat{K}_{heat}^t
\end{align}
where  $V(\mathbf{A}_s)$ is the Riemannian volume of $\mathbf{A}_s$. $V(\mathbf{A}_s)$ is parameterised by its radius $\omega$. \cite{niu2019} has provided a indication of the optimal order of magnitude of $\omega$ by minimizing the error of the kernel estimator. When $V(A_s)$ is large, the error of estimating the transition probability becomes smaller. But the error of approximating the transition density become bigger. The former is called Monte Carlo error and the later is called numerical error in Niu et. al (2019). The optimal order of magnitude of $\omega$ can be derived by balancing these two errors.

\begin{figure}[t]
 \centering
\includegraphics[width=0.45\textwidth,height=0.4\textwidth]{BM-density-example.png}
    \caption{ \label{fig:BM-density}
       {\footnotesize
    Three BM paths (red,blue,green) are simulated on the Swiss Roll, with the starting point at $s_0$ (black ball). $s$ (red ball) is the target point. $\bm{A}$ is the neighbourhood of $S$. Only the red path reach $\bm{A}$ at time $t$. The transition probability $p\{S(t)\in \bm{A} | S(0)=s_0 \}$ is 1/3.
     } }
\end{figure}

The Neumann boundary condition corresponds to BM reflecting at the boundary. This can be approximated by pausing time and resampling the next step until it falls into the interior of $\mathbb{M}$. This estimator is asymptotically unbiased and consistent \citep{niu2019}.  Note that $t$ is the BM diffusion time.  If $t$ is large, the BM paths have a higher probability to reach the neighbourhood of the target point leading to higher covariance and vice versa. The transition density can be estimated using Algorithm \ref{alg:bm}.

\begin{corollary} \label{corollary:metric}
The transition density estimate in \eqref{k:heat} is coordinate independent.
\end{corollary}
By Theorem \ref{theorem:sde}, it is straightforward to have Corollary \ref{corollary:metric}. The transition density estimate in \eqref{k:heat} of the stochastic process defined in \eqref{eqnBM} is also coordinate independent.
Based on Corollary \ref{corollary:metric},  we can evaluate the heat kernel of the implicit manifold estimated by different LVMs, independent of the specific parameterisations of their latent spaces.
If the estimated manifold is close to the true manifold, we expect the resulting BM transition density estimates to be similar to the ones estimated using the analytical parameterisation.

Here we take the Swiss roll as an example. Assuming the geometry of the Swiss roll is known, we can follow \cite{niu2019}'s approach and evaluate the heat kernel $K_{heat}^t(s_0,s)$ by simulating BM paths with the analytical metric tensor. The derivations of the analytical metric and parameterisation of the Swiss roll are shown in Appendix \ref{ax:swisspar}. Let $s_0$ with $radius=6$ and $width=3$ be the starting point of the BM. $N_{BM}=20000$ BM paths are simulated. The BM transition density is evaluated using \eqref{k:heat} at twenty nine target points $\{ s_j \in \mathbb{M} \subset \mathbb{R}^3| j \in \{1, \cdots,29 \} \}$ in the observation space. These target points are centred on $s_0$ and equally spaced. The diffusion time is fixed at 50. The results are plotted as the red solid line in Fig. \ref{fig:swiss-gl-kern}. The horizontal axis is the radius of the Swiss roll and the vertical axis is the transition density. The red density plot is asymmetric. When the radius is large, the transition density estimate decreases more quickly.

If the geometry of the Swiss roll is unknown,  Bayesian GPLVM and GPLVM can be applied to learn the metrics from 250 grid points on the Swiss roll. The derivations of B-GPLVM metrics and GPLVM metrics are given in section \ref{bgplvm} and \ref{gplvm}. The BM trajectories are simulated using the estimated metric tensors. The heat kernel estimates using B-GPLVM metrics are plotted as the green dashed line in Fig.\ref{fig:swiss-gl-kern}. It is clear that the B-GPLVM results are very close to the red solid line.
The estimates using GPLVM metrics are plotted in brown dashed line. It is also close to the solid red line, but not as good as the B-GPLVM estimates. GPLVM performs point estimate of $x$ in the latent space while B-GPLVM estimates a Bayesian posterior of $x$. As a result, B-GPLVM is more robust in terms of estimating the latent variables \citep{damianou2016}.

The heat kernel estimates from the Graph Laplacian (GL) approach (\cite{dunson2020diffusion}) in different data regime are also provided in Fig. \ref{fig:swiss-gl-kern}. When the number of points on the Swiss roll is 250, the GL kernel estimates are plotted as the blue dashed line. It is far from the solid red line and do not match the overall pattern of the analytical kernel estimates. When the number of grid points is increased to 1000 and 10000, the GL kernel estimates are plotted as the purple dashed line and black dash dotted line. The estimates are closer to the solid red line as the number of grid points increases. Comparing to the GL approach, the kernel estimates using B-GPLVM metrics achieve the best performance with much fewer points on the Swiss roll. B-GPLVM is used in the simulation study and real data applications in later sections.

\begin{figure}[h]
 \centering
\includegraphics[width=0.6\textwidth,height=0.5\textwidth]{GLkernel-c.pdf}
    \caption{ \label{fig:swiss-gl-kern}
       {\footnotesize
Comparison of heat kernel estimates using the analytical metric, B-GPLVM metric, GPLVM metric and Graph Laplacian.  The red solid line represents heat kernel estimates using the analytical metric. The green dashed line represents estimates using B-GPLVM metric. The brown dash dotted line represents estimates using GPLVM metric. The heat kernel estimates from the Graph Laplacian approach are plotted as the blue dashed line. We increase the number of the grid points on the Swiss roll from 250 to 1000 and 10000. The GL estimates are plotted as the purple dotted line and black dashed line.
     } }
\end{figure}

\begin{algorithm}[tb]
   \caption{ Simulating BM sample paths on $\mathbb{M}$ for estimating $K_{heat}^t$ }
   \label{alg:bm}
\begin{algorithmic}
  \STATE Learn the metric $\mathcal{G}$ from the point cloud $\mathcal{S} = \{  s_i | i = 1,\ldots,n+v\ \} $. \COMMENT {\footnotesize use eqn \eqref{eqn:exmetric}  }
   \STATE  1.1 Generate BM trajectories on implicit manifold.
   \FOR{$i = 1,\ldots,n$   \COMMENT{ \footnotesize $n$ is the size of data points } }
     \FOR{$j = 1,\ldots,N_{BM}$   \COMMENT{\footnotesize $N_{BM}$ is No. of trajectories }  }
    \FOR{$ l = 1,\ldots,N_t$   \COMMENT{\footnotesize $N_t$ steps Brownian motion, $N_t\times \Delta t$  $\rightarrow$ max diffusion time} }

 \STATE {\bf do}   \COMMENT{\footnotesize keep proposing x until it falls inside of the boundary }
   \STATE $q\left(x_{i,j}(l)|x_{i,j}(l-1)\right) = \mathbb N \left( x_{i,j}(l) |  \mu \left( x_{i,j}\left(l-1\right),\Delta t \right)  , \Delta t \mathcal{G}^{-1} \right )$  \COMMENT{ \footnotesize  use eqn \eqref{eqn:disProp} }
\STATE {\bf While} $\mathbb{V} ar( \phi(x_{i,j})  ) > \alpha$, $x_{i,j}$ is outside of the boundary $\partial \mathbb{M}$.   \COMMENT{ \footnotesize use eqn \eqref{eqn:varb}  }

    \ENDFOR
     \ENDFOR
    \ENDFOR
    \RETURN $\bm x$
\STATE 1.2 Given a discrete choice of the diffusion time $t\in \{ \Delta t, 2\Delta t,\cdot \cdot \cdot, N_t\Delta t \}$, the covariance matrix $\Sigma^{t}$ is estimated based on the BM simulation from Algorithm 1.1.
\FOR{$i = 1,\ldots,n$}
\FOR{$j = 1,\ldots,n$}
\STATE $N_{\mathbf{A}_j}$ = which( x(t) $\in$ $\mathbf{A}_j$ )  \COMMENT{ {\small counting how many BM paths reach $\mathbf{A}_j$}  }
\STATE $K_{heat}^t(s_i,s_j)= \frac{ N_{\mathbf{A}_j}}{N_{BM}*V( \mathbf{A}_j)}$  \COMMENT{ {\small use eqn \eqref{k:heat} } }
\STATE $\Sigma_{ij}^t = \sigma_h^2 K_{heat}^t(s_i,s_j)$
\ENDFOR
\ENDFOR
\RETURN $\Sigma^t$
\end{algorithmic}
\end{algorithm}

From  \eqref{k:heat} we can see the construction of GPUM and the heat kernel requires simulating BM sample paths at each data point. Although the BM
simulations are trivially parallelizable, the computational cost can be high when the
number of data points is large. \cite{niu2019} proposed the sparse intrinsic GP on known manifold by introducing some inducing points. The number of inducing
points is much smaller than the number of data points. BM paths only need to be
simulated starting at the inducing points instead of every data point. The inducing point approximation summarizes the training data into a small set of inducing points, so that
inference could be done more efficiently \citep{QuioneroCandelaRasmussen2005}. Similar approach can be applied in the
GPUM when the manifold is unknown. Small number of inducing
points can be introduced in the learned latent space. The focus of this paper has been on developing the GPUM. The development of the sparse GPUM is for the future research.

\subsection{Optimising the kernel hyperparameters} \label{sec:optimpar}
Given a diffusion time $t$, we can generate the covariance matrix $\Sigma_{ {\text {\bf ff}} }^{t}$ for training data $\mathcal{D}$ using Algorithm \ref{alg:bm}. $\Sigma_{ \text{\bf ff} }^t$ can be obtained as follows: with the $i_{th}$ data point as the starting point, $N_{BM}$ trajectories are simulated to generate the $i_{th}$ row of $\Sigma_{ \text{\bf ff} }^t$. For each element of $\Sigma_{ \text{\bf ff}}^t$, $\hat{K}^t(s_i,s_j)$ is then estimated using  \eqref{k:heat}. The hyperparameters can be obtained by maximizing the log of the marginal likelihood  (over $f$) in \eqref{loglike}.
The maximum BM diffusion time is set as $N_t \times \Delta t$. $\Delta t$ is the BM simulation time step as defined in \eqref{eqn:disProp}. $N_t$ is the number of simulation steps. $N_t$ covariance matrices $\Sigma_{ {\text {\bf ff}} }^{1 \ldots N_t}$ can be generated based on the BM simulations. Optimisation of the diffusion time $t$ can be done by selecting the corresponding $\Sigma_{ {\text {\bf ff}} }^{t}$ which maximizes the log marginal likelihood. Estimation of $\sigma_h$ follows standard optimisation routines, such as quasi-Newton.
\begin{align}
\label{loglike}
\log p( \bm{y}| s) &= \log \int p( \bm{y}| {\text {\bf f} }) p({\text {\bf f} }|s) d{\text {\bf f} }   \nonumber \\
&= -\frac{1}{2} \bm y^T (\Sigma_{ {\text {\bf ff}} }^{t} + \sigma_{noise}^2 I)^{-1} \bm y - \frac{1}{2} \log|\Sigma_{ {\text {\bf ff}} }^t +\sigma_{noise}^2I | - \frac{n}{2}\log2\pi.
\end{align}

Let ${\bf f}_*$ be a vector of values of $f(\cdot)$ at the unlabeled points in $\mathcal{V}$. Under the regression model in \eqref{eqn:reg} and the GPUM  prior, we have the joint distribution of ${\bf y} $ and ${\bf f}_*$:
\begin{align}
p( {\text {\bf y}}, {\text {\bf f}}_*) = \mathcal{N} \left ( 0,
 \left [ \begin{array}{cc}
 \Sigma_{{\text {\bf f}} {\text {\bf f}}} + \sigma_{noise}^2 I_n &\Sigma_{ {\text {\bf f}}  {\text {\bf f}}_*} \\[0.3em]
 \Sigma_{ {\text {\bf f}}_* {\text {\bf f}}} &\Sigma_{ {\text {\bf f}}_* {\text {\bf f}}_*} \\[0.3em]
                \end{array} \right ]
                \right ),
\end{align}
where $\Sigma_{{\text {\bf f}}_*{\text {\bf f}}}$ is the covariance matrix for training(labeled) and unlabeled data points. The predictive distribution is derived by marginalising out {\bf f}:
\begin{align}
p( {\text {\bf f}}_*| \bm y) = \int p( {\text {\bf f}}_*{\text {\bf f}}| \bm y) d{\text {\bf f}} = \mathcal N \left ( \Sigma_{ {\text {\bf f}}_*{\text {\bf f}}}  \left(\Sigma_{{\text {\bf f}}{\text {\bf f}}}  + \sigma_{noise}^2 I \right)^{-1}\bm y , \ \  \Sigma_{{\text {\bf f}}_*{\text {\bf f}}_*}-\left(\Sigma_{{\text {\bf f}}{\text {\bf f}}} +\sigma_{noise}^2 I \right)^{-1} \Sigma_{{\text {\bf f}}{\text {\bf f}}_*} \right)
\end{align}

\section{Simulation study on Swiss roll} \label{sec:simdata}

 In this section, we carry out a simulation study for a regression model with synthetic data on the Swiss roll which is a two dimensional manifold depicted by a point cloud in $\mathbb{R}^3$. The point cloud is plotted in Fig. \ref{swiss:point3d}. It is comprised of the set of labeled points $n=24$ and the set of unlabeled points $v=450$. Both labeled and unlabeled observed points are used in B-GPLVM to learn the latent space. The point cloud is unfolded into a flat surface as in Fig. \ref{swiss:var}. The unlabeled points are plotted as blue triangles and the labeled points are in red in the latent space. The variance of the mapping is plotted as the gray background of the latent space. It is clear from Fig. \ref{swiss:var} that the background color gets darker in the region further away from the blue triangles. This indicates the variance of the mapping
gets bigger in the region which is far from the observations.

 Based on the definition in \eqref{eqn:bound}, the boundary $\partial \mathbb{M}$ is plotted in Fig. \ref{swiss:boud}. The black regions on the left and right sides of Fig. \ref{swiss:boud} are outside of $\partial \mathbb{M}$. The white area in the middle is within $\partial \mathbb{M}$. The magnification factor is plotted as the background color in Fig. \ref{swiss:mag}. Similar to Fig. \ref{fig:BM}, the horizontal axis can be interpreted as the scaled radius of the Swiss roll.
 When the radius is bigger, the corresponding magnification factor is larger (the color is darker).  However, as there are fewer data points available at the tail of the Swiss roll, the estimation of the implicit manifold from B-GPLVM becomes less accurate. This can be observed on the right end of Fig. \ref{swiss:mag} where the last column of the observed latent points is further away from the rest. Once the metric $\mathcal{G}$ and the boundary $\partial \mathbb{M}$ are learned from B-GPLVM, we can estimate the heat kernel by simulating BM paths on the implicit manifold. An example of a BM sample path on Swiss roll is shown in Fig. \ref{fig:BM}.

For the labeled points, the response variables are
 \begin{align*}
y_i = f(s_i^1,s_i^2,s_i^3) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^3
 \end{align*}

where $f$ is the unknown regression function, $s_i$ is the coordinate of the observed point in $\mathbb{M} \subset \mathbb{R}^3$. For better visualisation, the true function is plotted in the unfolded Swiss roll in Fig.\ref{swiss:true} using a two dimensional analytical parameterisation. The coordinates are radius and width. 24 labeled observations are marked as black crosses. The true function values are indicated by the color codes and contours. The regression function varies slowly when the radius is small and rapidly when the radius is big.

We first apply the Euclidean $\mathbb{R}^3$ GP (the standard GP as in \cite{Rasmussen2006} Chapter 4) constructed by the squared exponential kernel with $\mathbb{R}^3$ Euclidean distance. With this kernel setting, the $\mathbb{R}^3$ GP completely ignores the interior structure of the manifold and lets the inner layer and outer layer of the Swiss roll interact. The predictive means on the unlabeled points are shown in Fig. \ref{swiss:3d}. Compared to the true function in Fig.\ref{swiss:true}, the overall shape of $\mathbb{R}^3$ GP prediction contours is more wiggly. The color coding of Fig.\ref{swiss:3d} is also very different in the regions where the radius (horizontal axis) is 6, 8 and 10. In the second model, the $\mathbb{R}^2$ Euclidean distance in the latent space is used in the squared exponential kernel to construct the $\mathbb{R}^2$ GP. The geometric properties such as the metric and magnification factors are ignored. The color coding of the predictive means for $\mathbb{R}^2$ GP is shown in Fig. \ref{swiss:r2} in Appendix \ref{ax:swisspar}. Since the regression function is nonstationary in the latent space, the $\mathbb{R}^2$ GP is underfitting. The GPUM  predictive mean is shown in Fig. \ref{swiss:gp}. The overall pattern of the  GPUM  prediction is similar to the true function. The shape of the contours and the color coding of Fig.\ref{swiss:gp} are consistent with Fig.\ref{swiss:true}. The prediction of the GL-GP is shown Fig.\ref{swiss:gl} in Appendix \ref{ax:swisspar}. Since the GL kernel estimates in Fig. \ref{fig:swiss-gl-kern} are far from the analytical kernel, the GL-GP prediction is also poor. A one dimensional comparison is generated by plotting the predictions of all methods at $z=4$ and varying the radius from 2 to 12.5 in Fig.\ref{swiss:compare1d}. It is clear the $\mathbb{R}^2$ prediction(brown dotted line) is under fitting. The $\mathbb{R}^3$ prediction(green dashed line) is oscillating in the opposite direction of the ground truth (blue solid line). The GL-GP prediction (black dashed line) is also oscillating around the blue solid line. The GPUM  prediction (red dashed line) achieves the best performance and follows the overall pattern of the ground truth.

\begin{figure}[H]
    \centering
\subfigure[Swiss roll point cloud in $\mathbb{R}^3$ ]{\label{swiss:point3d}  \includegraphics[width=0.44\textwidth,height=0.3\textwidth]{swiss1.png} }
\subfigure[ Latent space and variance of mapping]{ \label{swiss:var} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{swissVar1.png} }
\subfigure[Latent space and magnification factor]{ \label{swiss:mag} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{swissMag1.png} }
\subfigure[Boundary]{ \label{swiss:boud} \includegraphics[width=0.49\textwidth,height=0.35\textwidth]{swissBound1.png} }
    \caption{ \label{fig:swissDR}
    \footnotesize{ (a) The Swiss roll represented as a point cloud in $\mathbb{R}^3$. (b) The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and red crosses are labeled points. The background color represents the variance of the mapping, a dark background represents high uncertainty. (c) The same latent space visualization with the background color representing the magnification factor, a dark background represents high magnification factor. (d) The same latent space visualization highlighting the boundary. The dark region is outside of the boundary and the white region is inside of the boundary.}
  }
\end{figure}

In order to evaluate the performance of different methods in different data regimes, we consider three scenarios with different numbers of unlabeled grid points on the Swiss roll from $v=250$ to $v=450$ and $v=800$. These unlabeled grid points are used to estimate the manifold for GPUM , GM-GP, GL-GP and are also used as the test set of the regression accuracy. We constructed twenty sets of training points by randomly selecting $n=23$ labeled points from the labeled set. For each training set, GPUM , GM-GP, GL-GP, $\mathbb{R}^2$GP and $\mathbb{R}^3$GP have been applied to make predictions at the points in the test set. The root mean square errors (RMSE) are calculated between the true function values and the predictive means of all five models. The mean and standard deviation of the root mean square errors are reported in Table \ref{tb:swiss}. GPUM’s results are consistent
across all three scenarios and significantly better than all other methods. The GL-GP and GM-GP have similar performance. Both models perform better when the number of the grid points is large. They are significantly better than the Euclidean GPs. Comparing to the graph based approaches, GPUM achieves the minimum mean RMSE with fewer points on the manifold.

 \begin{figure}[H]
    \centering
    \subfigure[True function and data points]{\label{swiss:true}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPtrue.pdf}}
     \subfigure[ \small{$\mathbb{R}^3$ GP  prediction }]{\label{swiss:3d}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissGPr3.pdf}}
    \subfigure[GPUM  prediction]{\label{swiss:gp}  \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{swissmanifoldGP.pdf}}
    \subfigure[ Comparison of predictions when fixing z ]{ \label{swiss:compare1d} \includegraphics[width=0.45\textwidth,height=0.4\textwidth]{fig/SwissReg-z4new.pdf}
    \caption{\label{fig:swiss2d}
   { \footnotesize
(a) The true function in the unfolded Swiss roll. The labeled points are marked with black crosses. (b) The prediction of $\mathbb{R}^3$ GP. $\mathbb{R}^3$ GP is constructed using the Euclidean distance in $\mathbb{R}^3$. (c) The prediction of GPUM . (d) Comparing the predictions of all methods on a spiral in the swiss roll by fixing the z coordinate at 4 and changing the radius from 2 to 12.5. The vertical axis represent the value of the prediction. The horizontal axis represents the radius.
    }
 }
\end{figure}

\begin{table}[h]
\caption{\label{tb:swiss}Comparison of the root mean squared errors of five methods on Swiss roll. Values in parentheses show the standard deviation.}
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
    &   $\mathbb{R}^3 GP$ & $\mathbb{R}^2 GP$ &  GPUM  & GL-GP & GM-GP\\
\hline
mean RMSE $v=250$ &0.284(0.006) & 0.293(0.005)&0.163(0.020)& 0.243(0.003)&0.231(0.001) \\
mean RMSE $v=450$ & 0.298(0.007)  & 0.290(0.005) & 0.162(0.003) &0.220(0.002) & 0.207(0.002)\\
mean RMSE $v=800$ &0.287(0.006)&0.282(0.005)& 0.164(0.002)&0.216(0.001) & 0.206(0.001)\\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\begin{figure}[H]
    \centering
\subfigure[ Latent space and uncertainty of mapping ]{\label{wifi:latent}  \includegraphics[width=0.44\textwidth,height=0.3\textwidth]{wifi_Latent.png} }
\subfigure[ Latent space and boundary]{ \label{wifi:bound} \includegraphics[width=0.51\textwidth,height=0.33\textwidth]{wifi_bound.png} }
\subfigure[Latent space and magnification factor]{ \label{wifi:mag} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{wifi_mag.png} }
\subfigure[A BM trajectory]{ \label{wifi:bm} \includegraphics[width=0.45\textwidth,height=0.3\textwidth]{wifi_BMmag.png} }
    \caption{ \label{fig:wifi}
    \footnotesize{ (a)
The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and the red crosses are labeled points. The background color represents the variance of the mapping, a dark background represents high uncertainty. (b) The same latent space visualization highlighting the boundary. The dark region is outside of the boundary and the white region is inside of the boundary. (c) The same latent space visualization with the background color representing the magnification factor, a dark background represents a high magnification factor. (d) A BM trajectory is shown in purple. The red star is the starting location.}
  }
\end{figure}

\section{Location estimation from WiFi signal} \label{sec:wifi}

Indoor location estimation has tremendous value but standard location estimation techniques such as Global Positioning System (GPS) do not work indoors. Instead, indoor wireless signals from devices such as WiFi access points can be exploited for location estimation. In this section, we consider the problem of indoor 2D location estimation from WiFi access point signal strengths. We use the WiFi data collected by \cite{ferris2007}, in which a series of WiFi signal strength traces are collected by a mobile device which travels in a one floor university building. The 2D location coordinates of the mobile device are also recorded by a click-to-map based annotation program. The total number of WiFi access points in this dataset is 30. As it is often expensive to collect labeled data, we mimic a low data indoor location estimation scenario by assuming that the true locations of only three points are known and we aim at predicting the indoor locations of the remaining points based on the WiFi signals.

This location estimation problem can be treated as a regression problem, in which the location coordinate $y$ of the mobile device can be modelled as a function of the high dimensional WiFi signal
\begin{align*}
y_i = f(s_i) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^{30}
\end{align*}
where $f$ is the unknown regression function and WiFi signal $s_i$ is represented by a 30 dimensional vector. Here we consider the WiFi signal measurements at $n+v=36$ locations. Only $n=3$ of the locations are labeled with one dimensional location coordinates of the mobile device. To avoid selecting the three points clustered together, we randomly pick the points from three different regions respectively, where the union of the regions cover the whole dataset. Different methods are applied to estimate the coordinates of the mobile device in the testing sets (the unlabeled $v=33$ locations). Twenty training and testing sets are generated from this random selection.

Unlike the simulation study of the Swiss roll, we cannot plot the high dimensional point cloud of the WiFi signals. The implicit manifold is also unknown. We first estimate a $q=2$ dimensional latent space (the chart of the underlying manifold) using B-GPLVM. The value of $q$ is determined by the ARD contributions which measure how much each dimension is contributing to the latent space. The input dimensions are sorted based on the relevance assigned by the scaling in B-GPLVM. The plot of ARD contributions is shown in Appendix \ref{ax:wifi}. The latent space is plotted with the variance of the mapping as the gray background in Fig. \ref{wifi:latent}. The 36  WiFi signal strength measurements are represented by the blue triangles. The training set (the labeled points) is marked by the red crosses. The dark color value is for high uncertainty. Since the mobile device moves in a loop closure, the latent points forms a closed loop in the latent space. The boundary of the implicit manifold $\partial M$ is shown in Fig. \ref{wifi:bound} which is defined by \eqref{eqn:bound}. The magnification factor is  plotted as the gray background in Fig.\ref{wifi:mag}. The dark color value is for high magnification factor. A sample path of BM in the implicit manifold is plotted in Fig. \ref{wifi:bm}. With the Neumann boundary condition the BM path can only exist within the boundary.

The coordinates of the mobile device are estimated by five different methods. The ground truth values are marked by different colors in the latent space in Fig. \ref{wifi:true}. If we ignore the interior structure of the manifold, a $\mathbb{R}^{30}$ GP using the squared exponential kernel with $\mathbb{R}^{30}$ Euclidean distance of the WiFi signals is applied. The predictive means of the $\mathbb{R}^{30}$ GP are plotted in Fig.\ref{wifi:hd}. It is clear the $\mathbb{R}^{30}$ GP prediction is very poor. The color coding of the prediction is different from the ground truth. The GPUM   predictive mean is shown in Fig. \ref{wifi:in-GP}. The overall pattern of the GPUM  prediction is similar to the ground truth. In the third case, the $\mathbb{R}^2$ Euclidean distance in the latent space is used with the squared exponential kernel to construct the $\mathbb{R}^2$ GP. The prediction results are shown in Fig.\ref{fig:wifi2d} in Appendix \ref{ax:wifi}. The Graph Laplacian based approaches such as GL-GP and GM-GP have also been applied. Since the number of the observations are relatively low ($n+v=36$), the graph based methods result in poor approximation of the implicit manifold. The GL-GP's performance is similar to the $\mathbb{R}^{30}$GP. The mean and standard deviation of the root mean square errors of the twenty testing sets are also calculated in Table \ref{tb:wifi} for all five models. GPUM  significantly outperforms the other four models and achieves the minimum mean RMSE.

\begin{table}[t]
\caption{\label{tb:wifi} Comparison of the root mean squared errors of five methods on WiFi signal data. Values in parentheses show the standard deviation. }
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
    & $\mathbb{R}^{30}$GP & $\mathbb{R}^2$GP & GPUM  &  GL-GP &  GM-GP \\
\hline
mean RMSE & 5.57(1.43)  & 4.83(1.83) & 4.11(0.88) & 5.6(1.15) & 6.04 (0.57) \\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\begin{figure}[t]
    \centering
    \subfigure[Ground truth]{\label{wifi:true}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_true.png}}
     \subfigure[ \small{ $\mathbb{R}^{30}$ GP prediction }]{\label{wifi:hd}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_30d.png}}
    \subfigure[GPUM   prediction]{\label{wifi:in-GP}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{wifi_mani.png}}

    \caption{\label{fig:wifiGP}
   {\footnotesize
    Comparison of GPUM  and Euclidean GPs in WiFi signal example. (a) The ground truth is plotted with color at the observed points in the latent space. (b) $\mathbb{R}^{30}$ GP prediction. $\mathbb{R}^{30}$ GP is constructed using the Euclidean distance of WiFi signals in $\mathbb{R}^{30}$. (c) GPUM  prediction.
    }
 }
\end{figure}

\section{Camera angle estimation from images} \label{sec:real}
In this section, we consider the problem of estimating the camera angle from images. We consider the setting in which an object is placed on a turntable and a set of images are taken at different angles with respect to the camera. We aim at recovering the camera angles associated with individual images by knowing the true camera angles for some of the images. We use the images from the COIL data set \citep{nene1996}. Here we consider 66 images of object-14 (a toy cat). The camera angles range from 15 degrees to 340 degrees.
The raw images are converted to grayscale and downscaled to $32 \times 32$.
The raw pixels of each image are flattened into a 1024 dimensional vector.
The dimension of the original image space is $p=1024$. Six image examples are given in Fig. \ref{coil:image}. We estimate a $q= 2$ dimensional latent space using Bayesian GPLVM. The two dimensional latent space is plotted with the variance of the mapping as the gray background in Fig. \ref{coil:var}. A dark background represents a high uncertainty in the corresponding region. The unlabeled data points are marked as blue triangles and the labeled data points as red stars. The overall shape of the latent points in Fig.\ref{coil:var} looks like a ring with a gap in the lower right. The magnification factor is plotted as the gray background in Fig. \ref{coil:mag}. A dark background represents a high magnification factor. The boundary of the implicit manifold $\partial \mathbb{M}$ is shown in Fig.\ref{coil:bound} in Appendix \ref{ax:coil}.
\begin{figure}[t]
    \centering
\subfigure[ COIL images ]{\label{coil:image} \includegraphics[width=0.31\textwidth,height=0.28\textwidth]{coilimage.png} }
\subfigure[ Latent space and variance of mapping]{ \label{coil:var} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{ncoilLatentVar.png} }
\subfigure[Latent space and magnification factor]{ \label{coil:mag} \includegraphics[width=0.31\textwidth,height=0.3\textwidth]{ncoilLatentMag.png} }
    \caption{ \label{fig:coil}
    \footnotesize{ (a) Examples of six COIL object images. (b) The latent space constructed from Bayesian GPLVM. The blue triangles are unlabeled points and red stars are labeled points. The background color represents the variance of the mapping. The dark color represents high uncertainty. (c) The same latent space visualization with the background color representing the magnification factor. The dark color represents a high magnification factor.
      }
  }
\end{figure}

The scaled camera angle $y$ is modelled as a function of the image.
 \begin{align*}
y_i = f(s_i) + \epsilon_i, \ \ i=1\ldots n, \ \ s_i\in \mathbb{R}^{1024} .
 \end{align*}
where $f$ is the unknown regression function, and $s_i$ is represented by a 1024 dimensional vector. n=13 images are randomly selected and used as the training set (labeled data). They are plotted as the red stars in the latent space in Fig.\ref{coil:var}. The remaining 53 images are used as the testing set (unlabeled data)  which are plotted as the blue triangles in the latent space in Fig. \ref{coil:var}. The true angle values are marked by different colors on the observed points in the latent space in Fig. \ref{coil:true}. Different methods have been applied to estimate the scaled camera angles. The $\mathbb{R}^{1024}$GP using the squared exponential kernel with $\mathbb{R}^{1024}$ Euclidean distance is applied to the image data first. The predictive means of the $\mathbb{R}^{1024}$GP are plotted in Fig.\ref{coil:hd}. The color coding of the prediction at the lower right of the plot is bright green which is different from the truth. Ignoring the boundary and the magnification factor in the latent space, the $\mathbb{R}^2$GP using the squared exponential kernel with $\mathbb{R}^2$ Euclidean distance in the latent space is applied to the image data. The predictive means are shown in  Fig.\ref{coil:2d} in Appendix \ref{ax:coil}. The predictive means of GPUM  are plotted in Fig.\ref{coil:in-GP}. The overall pattern of the GPUM  prediction is very similar to the true function. Since the boundary of the implicit manifold is defined in the lower right region of the latent space, the GPUM  does not smooth across the boundary and gives a better prediction. Ten training and testing
sets are generated from the random selection.
The mean and standard deviation of the root mean square errors are calculated for the ten testing sets for all methods in Table \ref{tb:coil}. The GPUM achieve the smallest mean RMSE. It is significantly better than  the GL-GP and GM-GP. The difference between the GPUM  and the $\mathbb{R}^{1024}$GP is not significant. The boxplot of the root mean squared errors are shown in Fig.\ref{coil:box} in Appendix \ref{ax:coil}.

 \begin{figure}[t]
    \centering
    \subfigure[True function]{\label{coil:true}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilTrue.png}}
     \subfigure[ \small{ $\mathbb{R}^{1024}$ GP prediction }]{\label{coil:hd}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilHd.png}}
    \subfigure[GPUM  on manifold prediction]{\label{coil:in-GP}  \includegraphics[width=0.315\textwidth,height=0.3\textwidth]{coilIn-GP.png}}

    \caption{\label{fig:coilGP}
   { \footnotesize
     Comparison of GPUM  and Euclidean GPs in COIL images example. (a) The true function is plotted with color at the observed points in the latent space. (b) $\mathbb{R}^{1024}$ GP prediction. $\mathbb{R}^{1024}$ GP is constructed using the Euclidean distance of data points in $\mathbb{R}^{1024}$. (c) GPUM  prediction.
    }
 }
\end{figure}

\begin{table}[t]
\caption{\label{tb:coil}Comparison of the root mean squared errors of five methods on COIL images. Values in parentheses show the standard deviation.}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccr}
\hline
&$\mathbb{R}^{1024}$GP &$\mathbb{R}^2$GP  & GPUM  &  GL-GP& GM-GP  \\
\hline
mean RMSE & 0.097(0.040) &0.094(0.042) & 0.060(0.038)  & 0.116(0.027) & 0.143(0.026) \\
\hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\section{Discussion} \label{discuss}

Our work provides a novel framework for regression on implicit manifolds embedded in high dimensional point clouds. The geometry of the implicit manifold is learned by probabilistic latent variable models. This gives a distribution over a smoothly changing local metric at each point in the latent space.  The boundary of the implicit manifold is defined according to the uncertainty of the mapping. The expression for the expected local metric is derived and used to simulate BM sample paths on manifolds. We have proved the BM transition density estimation is coordinate independent in section \ref{sec:density}. This allows us to compare the transition density estimates from different parameterisations of the same manifold. The BM simulation using the B-GPLVM metric gives similar transition density estimation results as the BM simulation using the analytical metric. GPUM  is constructed by employing the equivalence relationship between heat kernels and the transition density of BM on manifolds. This allows the GPUM  to incorporate the intrinsic geometry of the implicit manifold for inference while respecting the interior constraints and boundary. The experiment results in section \ref{sec:simdata}, \ref{sec:wifi} and \ref{sec:real} indicate that GPUM  achieves significant improvements over Euclidean GPs and Graph Laplacian based GPs. Although the simulation of BM paths can be easily parallelised, Algorithm \ref{alg:bm} can still be computational expensive if the sample size is large. The number of the sample paths can be reduced by leveraging the idea of pseudo data, which has been widely studied in the literature of sparse GPs \citep{QuioneroCandelaRasmussen2005,niu2019}.