Geometry-Aware Hierarchical Bayesian Learning on Manifolds

Our kernel derivation originates from the periodic potential diffusion process. It is a special case of the reaction diffusion process. A reaction diffusion process $T(v,t)$ is the solution to a reaction diffusion equation [18, 57]:

where $\Delta$ is the Laplace operator, constant $\alpha$ is 1. The initial condition is 0. $F$ is the reaction function that defines the property of the energy source [17]. Given a certain $F$, there exists a corresponding physical scenario [57, 42]. The reaction function $F$ can be expressed as the multiplication of a Dirac delta function at location $v$ and a temporal function $h(t)$: $F=h(t)\delta(v-v^{\prime})$. By defining the Green’s function of Laplace operator under the Dirichlet boundary condition as $G(v,v^{\prime},t)$ [3], $T$ is equal to:

Reminding that the Green’s function in an $\mathbb{R}^{d}$ diffusion problem has the standard form: $G=\frac{e^{-v^{2}/4t}}{(4\pi t)^{d/2}}$. When $h(t)$ is periodic: $h(t)=cos(\omega t)$, Eq. (26) is further derived as:

Eq. (3) is called the periodic potential diffusion process, which is the theoretical foundation of our kernel.

A GP regression (GPR) aims at learning a multivariate distribution that fits with the training data and predicts the observation $y_{n+1}$ when a testing sample $x_{n+1}$ arrives [48]. The Bayes’ theorem is used to transform the prior knowledge to posterior inference in the learning process. As known, every finite marginal distribution of a GP still follows a multivariate Gaussian distribution. Therefore, the predictive distribution $\mathcal{N}(m^{\prime},K^{\prime})$ can be uniquely determined by the standard rules for conditioning Gaussian distributions:

When new samples are continuously given, the GP is recursively updated. Later, we adopt the framework of GPR in a salient point selection algorithm. Each salient point is taken as a sample. We update the saliency map after adding the previous selection into the prior and then select the next one until a certain number of salient points are collected.

A DGP is a deep belief network that hierarchically concatenates multiple Gaussian process latent variable models together (GP-LVMs) [14]. It mimics the composition of restricted Boltzmann machines (RBMs) in NNs. The sparse variational inference is usually used in GPR to estimate the posterior and avoid the cubic complexity [52]. Suppose $M$ inducing points $Z=\{z_{1},...,z_{M}\}(M\ll N)$ are selected, the complexity is decreased to $\mathcal{O}(M^{2}N)$ in a single GPR. For a DGP, the doubly stochastic variational inference is often applied to estimate the posterior [49, 5]. Specifically, the sparse variational inference is used to simplify the correlations within layers and keep the correlations between layers unchanged. In a DGP with $L$ layers, the prior is recursively defined on a series of vector-valued stochastic functions $F=\{F^{1},...,F^{L}\}$. The $i^{th}$ row of $F^{l}$ is denoted as $f_{i}^{l}$. Function values at inducing points $Z$ are $U$. Each single function has an independent Gaussian prior and inducing points. A joint density of a DGP can be expressed as:

According to the theories of variational inference, a factorized form of the posterior joint density is defined as [49]:

where $q(U^{l})$ is a Gaussian with mean function $m^{l}$ and covariance function $S^{l}$ for layer $l$. Eq. (9) indicates that the prediction of the $l^{th}$ layer, $f^{l}$, depends on the previous prediction $f^{l-1}$ and the inducing points of the current layer. By marginalising the approximation $q(U^{l})$ from each layer, the $i^{th}$ factorized variational posterior of the final layer is the integral of all paths $(f_{i}^{1},...,f_{i}^{L})$ through the Gaussian distributions defined by parameters $m^{l}$, and $S^{l}$:

The objective function is the doubly stochastic evidence lower bound (ELBO):

where $\mathcal{KL}$ is the Kullback–Leibler divergence. The ELBO has the complexity $\mathcal{O}(M^{2}N(D^{1}+...+D^{L}))$ to compute, $D^{l}$ is the size of the $l^{th}$ layer. The variational expection likelihood $\mathbb{E}$ in Eq. 11 is computed using the Monte Carlo approximation. Please refer to [49] for more details.

In this section, we start by briefly explaining what motivates us to choose the periodic potential diffusion process as the theoretical basis. Then, we provide two implementations of the kernel for different applications. Furthermore, we introduce our theoretical analysis about the property of the kernel. In the end, a hierarchical Bayesian model is defined.

The idea of using periodic potential diffusion process comes from two observations:
The first observation is that the integral Laplace transform of function $f(t)=t^{d-1}e^{-\frac{1}{4}at}$ in $\mathbb{R}^{d}$ in Eq. (12) (which is a special upper incomplete gamma function, $a$ is a constant) has several similar terms with the Matérn kernel in Eq. (13):

Both equations have modified Bessel function of the second kind $\mathcal{K}_{d}$, and the rest parts are functions with the same dimension order $d$. This indicates the possibility of deriving a kernel from Eq. (12).

Summarizing these two observations and incorporating the background in Sec. 3.1, our goal is to derive a close-form expression from Eq. (3) and prove that this expression is PSD regarding its spatial variable $v$. Unfortunately, the integral in Eq. (3) has no explicit solution according to [50].

Alternatively, we apply the same approximation strategy in [21] to estimate the Eq. (3) as the combination of a cosine Fourier transform$\hat{f_{c}}(\omega)$ and a sine Fourier transform$\hat{f_{s}}(\omega)$:

By solving this approximated form, we have the closed-form periodic potential diffusion process:

For simplicity, we define a frequency term $\lambda=\sqrt{\frac{1}{2}\omega}$ and a phase term $\phi=\omega t$. $\lambda$ and $\omega$ are hyper-parameters that are determined by regular parameter tuning methods. Reminding that the diffusion process is dynamic, we can express it as the accumulation of values at $N$ time slots. Therefore, the final kernel definition is expressed as an additive kernel [15]:

For regression tasks, the implementation in Eq. (16) is sufficient. But for better fitting with a hierarchical Bayesian model, we further introduce an implementation by taking the kernel as an ARD covariance function [24]. An individual length-scale parameter $\alpha$ is added for each input dimension which determines the relevancy of the input to the task:

GPs defined with Eq. (16) or Eq. (17) can be taken as the mixture of GPs. Previous studies show that the mixture of Gaussian has a universal approximation ability in fitting with continuous distributions [58, 43].

The next goal is to prove that Eq. (16) is PSD regarding its spatial variable. Eq. (16) is the composition of two functions: $e^{-\lambda\left\|v\right\|}$ and $cos(\lambda\left\|v\right\|+\phi)$. It is acknowledged that the exponential function and the cosine function are PSD regarding variable $v$. According to the composition property of a PSD function, the result of positive real function/constant times a PSD function is still PSD [4]. So, Eq. (16) is spatially PSD. Eq. (16) can be taken as a compositional kernel.

Eq. (1) indicates a family of kernels. This opinion is generalized as:

The kernel implemented in Eq. (16) or Eq. (17) is named as the geometry-aware convolutional (GAC) kernel. The GAC kernel satisfies the two characteristics discussed in Introduction, which are summarized as two lemmas:

We validate the proposed GAC kernel in two different studies. The first one is to use it in a regular GPR model. An unsupervised salient point selection algorithm will be introduced. This implementation fully utilizes geometry-awareness property. The second one is to adopt it in a Bayesian network layer. Hierarchical deep Bayesian learning models will be discussed. The feature aggregation will benefit from the nice intra-kernel convolution property.

GP regression (GPR) is widely used in spatial inference known as the “Kriging” method [48, 22]. Being inspired by the landmarking in face recognition [61], we use a GPR-based unsupervised salient point selection to process the irregular inputs. This method is essentially a manifold learning technique [31, 37, 47]. Therefore, applying our method also achieves nonlinear data dimensionality reduction. One advantage of GPR is the availability of uncertainty estimation [48]. We leverage this advantage and define the saliency score as the variance-based uncertainty [65]. By iteratively selecting new salient points and adding previously selected points to the prior, we successively collect a set of salient points. A geometry-aware kernel guarantees the salient points are significant to represent the original massive data. This strategy has been successfully applied in [37, 22, 20].

Suppose a set of $\kappa$ salient points is denoted as $\tilde{v}=\{\tilde{v}^{1},...,\tilde{v}^{\kappa}\}$. Define $\lambda=\sqrt{0.2\pi n},n=[1,...,N_{fre}]$. $\phi$ equals to a $N_{fre}$-length vector by dividing $[0,\frac{\pi}{2}]$ into $N_{fre}$ equal line-spaces. A multi-frequency multi-phase GAC kernel ($MMK$) is defined in a weighted squared form: $MMK=K\times W\times K$ ($K$ is symmetric). The weight $W$ is a diagonal matrix with the sum of absolute values of each row in GAC kernel $K$ as the diagonal entries: $W(v)=\sum\left|K(v,\cdot)\right|$. The saliency score $\Sigma_{\mathcal{M}}$ of $v_{i}$ during selecting the $({\kappa+1})^{th}$ salient point is defined as:

Only the point with the highest uncertainty score is selected as the salient point: $\tilde{v}:=argmax_{v}\Sigma$. The first salient point is the vertex with the maximum variance in $MMK$. From Eq. (18)-(20) we can see that a newly-selected salient point will be added to the prior knowledge and the next saliency score is determined by the previous selections. The whole process follows a GPR framework. The algorithm is summarized in Algorithm 1. Figure 1(b) shows examples of selecting salient points on point clouds. Figure. 4 demonstrates salient points on triangle meshes. Further evaluation results are available in Experiments.

We define a GAC-GP layer with the GAC kernel and follow the framework of DGPs to construct a hierarchical Bayesian learning model by stacking up multiple GP layers. Thanks to the intra-kernel convolution property, the GAC-GP layer has a good feature aggregation ability. In a pure hierarchical Bayesian learning model on manifolds, a computational pipeline is shown in Figure. 3. The first step is to process the input. Assume $\kappa$ salient points are selected by Algorithm 1. The feature on each salient point $f_{\tilde{v}}$ is a vector of length $l$. For shape $i$, we link all features of salient points in the order of their selections: $f_{i}=\{f_{\tilde{v}^{1}}...f_{\tilde{v}^{\kappa}}\}$. Noting that all shapes here belong to the same dataset and all salient points are selected with the same parameter setting in Algorithm 1. Otherwise, point-to-point registration is needed to concatenate features. Suppose $H$ shapes are used, the input $X$ is an $H\times(\kappa\times l)$ matrix. We compose a sequence of layers that map the input $x_{i}$ to its label $y_{i}$ in a hierarchical Bayesian model for classifications:

The output of hidden layer $l$ is a vector of the size $1\times S^{l}$, where $S^{l}$ is the layer size. This is similar to the relationship of the input channel and output channel in a NN. When the batch processing is applied, the output of each hidden layer has dimension $B\times S^{l}$, $B$ is the batch size. A final layer is appended with a softmax multi-class likelihood. The output vector has the dimension $1\times C$, where $C$ is the number of classes. Each entry stands for the probability belonging to a certain class. Arbitrary numbers of GP layers can be added as hidden layers. We use the doubly stochastic variational inference approach to estimate the posterior [49]. The optimization process is to maximize the ELBO in Eq. (11). The K-means method is used to choose inducing points.

Because the input of the Bayesian model is the point-wise features on manifolds, and NNs are strong in feature learning, we are inspired to further explore the potential of NN+Bayesian methods. Such a mixed model can take advantage of both the feature learning ability of NNs and the feature aggregation ability of Bayesian models. We generally follow the pipeline in deep kernel learning [63]. The input format is determined by the NN part. The output of NNs is the shape feature. The feature is then fed into a Bayesian model, and the following processing is the same as the pure Bayesian method. The negative marginal log-likelihood (MLL) is used as the loss function.

In the first experiment, the task is to select salient points on manifolds. The purpose is to evaluate the geometry-awareness of the GAC-GP and its stability in continuous regressions. Additionally, we compare the computational efficiency by recording the average running time. When defining $MMK$, we use the fast marching algorithm to compute geodesic distances and only select $k=200$ nearest neighboring points for each vertex. $N_{fre}$ is set as 5.

Three datasets are used [8]: (i) “Mandibular molars”, or “molar”, contains 116 teeth shapes; (ii) “First metatarsals”, or “metatarsal”, contains 57 shapes; and (iii) Distal radii contains 45 shapes. All shapes are triangle meshes with around 5000 vertices. Examples of each dataset are illustrated in Figure 4(a)-(c). The ROIs of such data are usually the marginal ridges, teeth crowns, and outline contours where the geometric features are rich [8, 12]. Therefore, the meshes in Figure. 4(a)-(c), are rendered by the normalized mean curvature as the ground-truth [27, 26]. The salient points are expected to evenly distribute in yellow regions.

Comparison methods include: (1) RBF kernel GP (RBF, RBF-GP) [48]. RBF is a classical choice; (2) spectral mixture kernel GP (SM, SpectralMixture, SMK-GP) [62]. SMK-GP had good performances in many vision tasks. 10 mixtures are used; (3) Matérn $3/2$ kernel GP (Matérn, Matérn-GP) [55]. Matérn $1/2$ and $5/2$ are excluded because of worse results; (4) mesh saliency (MeshSaliency, MS) [34]. It is a classical saliency detection method; (5) Weighted GP (Weighted, W-GP) in [23]. W-GP is a state-of-the-art GP method on manifolds. The same parameter settings are used.

Figure 4(a)-(c) demonstrate salient point selection results. In (a) and (b), we focus on comparing ours with the W-GP because it is the only GP kernel method that stresses the geometry-awareness, and it outperforms all other comparison methods. We show 20, 80, and 140 salient points. When selecting a small number of salient points, both methods present reasonable results. But the accuracy of W-GP gradually drops with the iterations increasing while GAC-GP shows a much better stability of geometry-awareness. For other comparison methods, we give examples of 50 salient points on one Molar shape in Figure 4(c). Figure. 1(b) shows some examples of salient points on point clouds.

Numerically, we define an Accumulated Curvature (AC) value to measure the selection performance: $AC_{N}=log\sum_{k=1}^{\kappa}(\left|GC_{k}\right|+\left|MC_{k}\right|)$, where $GC$ and $MC$ are normalized Gaussian curvature and mean curvature, and $\kappa$ is the total number of salient points. Drawing the AC values with the increased number of salient points forms an AC curve. Within a dataset, we compute the log average of AC values to plot an average AC curve. This curve reflects the geometry-awareness of different selection methods. When only selecting points with the largest accumulated curvatures at each iteration and draw the AC curve with these maximum values, we get a MaxGC+MC AC curve. This curve stands for the upper bound of saliency selections. The higher and closer to this MaxGC+MC curve an AC curve is, the better the geometry-awareness ability a corresponding method has. The results are plotted in Figure 5(a)-(c). The MaxGC+MC AC curve is drawn in dashes. Generally, AC curves of GAC-GP are above all other methods. Both the empirical visualizations and the numerical measurements verify that the GAC-GP is capable of making geometry-aware inference on manifolds. More importantly, its geometry-awareness is still consistent and stable after continuous regressions.

The computational efficiency is measured by averaging the running time of selecting one salient point, as shown in Figure 5(d). GPR-based methods gradually slow down due to the increased prior knowledge as shown in Eq. (18)-(20). The GAC-GP enjoys strong computational efficiency considering its superior performance because usually only a small number of salient points are needed.

In the second experiment, the task is to classify different human poses modeled by triangle meshes. The first purpose is to further evaluate the salient point selections by fixing the Bayesian learning architecture. The second purpose is to fix the inputs and evaluate different hierarchical Bayesian learning architectures. The pipeline in Figure. 3 is used. We choose the scaled-invariant wave kernel signature (SIWKS) [36, 2] as the feature. When computing the SIWKS, 30 smallest eigenvalues are used. The other parameter settings are the same as those in [2]. Features of 50, 100, and 250 salient points are used. In customizing the Bayesian model, we use the multitask variational strategy and the softmax likelihood in GPytorch. The number of inducing points is 50. We use Adam as the optimizer with an initial learning rate of 0.001. After the first 200 epochs, the learning rate changes to $10^{-4}$. We trained for 2000 epochs. The cost function is the variational ELBO mentioned in Sec 3.3. The comparison methods are the same as the prior experiment.

The SHREC14 non-rigid 3D human model is used [44]. It contains 400 triangle meshes of 40 human subjects with 10 poses. Each mesh has about 15000 vertices. The dataset is randomly split into: 90% for training, 5% for validation, and 5% for testing.

For the first purpose, we fix the Bayesian model to be a one-layer GAC-GP and feed in features from different methods. Table 1 shows the results. Classification with the features of all points has an accuracy of 0.910. Taking this value as a reference, we can draw conclusions that (1) our strategy of selecting salient points works for distinguishing different shapes. When enough salient points are selected, it is possible to use a small subset to represent the original data; (2) the geometry-aware selection of GAC-GP is more distinguishable than other comparison methods.

For the second purpose, we fix the inputs to be GAC-GP salient features and evaluate different Bayesian learning architectures. Here we use GCGP [60] as a comparison method. The results are shown in Table. 2. Noting that the code of GCGP is not available and we use our implementations, so we put a star mark on GCGP’s result. We can see that the accuracy is generally increased after adding a GAC layer, supporting a strong feature aggregation property. Meanwhile, a hierarchical concatenation of GAC layers shows a better accuracy than the single layer structure.

In the third experiment, the task is to classify different point cloud models. Our purpose is to demonstrate the work of integrating NNs with Bayesian learning. Here, we use the hierarchical feature learning architecture in PointNet++ [46] to learn the pointwise feature. We perform multi-class classification on ModelNet40 which contains 12311 3D CAD models of 40 categories. Each point cloud has 10000 points. We use 9843 models for training and 2468 models for testing. In the feature aggregation part, we use one single GAC-GP layer (ten mixtures). 64 inducing points are used. The optimizer is Adam and the initial learning rate is 0.04. The comparison methods include PointNet++, PointNet++ with normal information, and the Pointwise Convolutional NNs (PCNN) [30]. Table 3 shows that (1) the mechanism of NN+Bayesian can be jointly trained for tasks on manifolds; (2) models with Bayesian aggregation layers generally outperforms the classical multiple fully connected layers in our tests. We notice that the performance gain of using single GAC layer shrinks after adding normal information. Our hypothesis is that the features become more complicated, and the inference capability of single GAC layer is not powerful enough to well aggregate the new features. By adding 2&3 GAC layers, the improvements increase to 0.9% and 1.2%, respectively. The overall results show that architectures with GAC layers universally perform better than their original versions, which proves that such a co-design benefits the performance. A reasonable outlook is to investigate more effective architectures that integrate both methods for end-to-end tasks on manifolds.

We provide more salient point selection results on point clouds here. Figure. 7 shows twelve examples of salient point selection on McGill 3D shape benchmark [53]. The shapes have different sampling densities. The original data type is the triangle mesh. We remove their edge connections and only use the vertex coordinates as inputs. The results show that the GAC-GP can learn the geometric property of the inputs in the prior and the selected salient points are representative in distinguishing the shapes.

Figure. 8 illustrates two examples of salient point selection on Modelnet40 dataset. Figure. 8 (a) demonstrates the saliency maps after selecting 200 salient points. Figure. 8 (b) shows salient points on the original point clouds. Figure. 8 (c) is 200 salient points. We can see that 200 salient points can generally depict the original shapes. Noting that we did not use salient point selection algorithm in the point cloud classifications in the main submission. This is mainly because the comparison methods used the whole data as inputs, for fairness and convenience, we also use 10000 vertices as inputs in the experiments.