The parametrix kernel K^prx previously used in the literature is just a Gaussian RBF function with θ=arccosx^·y^ as the radial distance [1]:

Definition 3. The parametrix kernel is a non-negative function

defined for t > 0 and attaining global maximum 1 at θ = 0.

Note that this kernel is assumed to be restricted to the positive orthant. The normalization factor (4πt)-n-12 is numerically unstable as t ↓ 0 and complicates hyperparameter tuning; as a global scaling factor of the kernel can be absorbed into the misclassification C-parameter in SVM, this overall normalization term is ignored in this paper. Importantly, the parametrix kernel K^prx is merely the Gaussian multiplicative factor without any asymptotic expansion terms in the full parametrix expansion G^prx of the heat kernel on a hypersphere [1, 12], as described below. The leading order term in the full parametrix expansion was derived in Berger et al. [12] and quoted in Lafferty and Lebanon [1]. For the sake of completeness, we include the details here and also extend the calculation to higher order terms.

The Laplace operator on manifold M equiped with a Riemannian metric g_μν acts on a function f that depends only on the geodesic distance r from a fixed point as

where g ≡ det(g_μν) and ′ denotes the radial derivative. Due to the nonvanishing metric derivative in Equation (1), the canonical diffusion function

does not satisfy the heat equation; that is, (Δ − ∂_t)G(r, t) ≠ 0 (Supplementary Material, section 2). For sufficiently small time t and geodesic distance r, the parametrix expansion of the heat kernel on a full hypersphere proposes an approximate solution

where the functions u_i should be found such that K_p satisfies the heat equation to order t^p−d/2, which is small for t ≪ 1 and p > d/2; more precisely, we seek u_i such that

Taking the time derivative of K_p yields

while the Laplacian of K_p is

One can easily compute

and

The left-hand side of Equation (3) is thus equal to G multiplied by

and we need to solve for u_i such that all the coefficients of t^q in this expression, for q < p, vanish.

For q = −1, we need to solve

or equivalently,

Integrating with respect to r yields

where we implicitly take only the radial part of logg. Thus, we get

as the zeroth-order term in the parametrix expansion. Using this expression of u₀, the remaining terms become

and we obtain the recursion relation

Algebraic manipulations show that

from which we get

Integrating this equation and rearranging terms, we finally get

Setting k = 0 in this recursion equation, we find the second correction term to be

From our previously obtained solution for u₀, we find

and

Substituting these expressions into the recursion relation for u₁ yields

For the hypersphere S^d, where d ≡ n − 1 and g = const. × sin^{2(d − 1)}r, we have

and

Thus,

Notice that u₁(r) = 0 when d = 1 and u₁(r) = u₀(r) when d = 3. For d = 2, u₁/u₀ is an increasing function in r and diverges to ∞ at r = π. By contrast, for d > 3, u₁/u₀ is a decreasing function in r and diverges to −∞ at r = π; u₁/u₀ is relatively constant for r < π and starts to decrease rapidly only near π. Therefore, the first order correction is not able to remove the unphysical behavior near r = 0 in high dimensions where, according to the first order parametrix kernel, the surrounding area is hotter than the heat source.

Next, we apply Equation (4) again to obtain u₂ as

After some cumbersome algebraic manipulations, we find

Again, d = 1 and d = 3 are special dimensions, where u₂(r) = 0 for d = 1, and u₂(r) = u₀/2 for d = 3; for other dimensions, u₂(r) is singular at both r = 0 and π. Note that on S¹, the metric in geodesic polar coordinate is g₁₁ = 1, so all parametrix expansion coefficients u_k(r) must vanish identically, as we have explicitly shown above.

Thus, the full G^prx defined on a hypersphere, where the geodesic distance r is just the arc length θ, suffers from numerous problems. The zeroth order correction term u0=(sinθ/θ)-n-22 diverges at θ = π; this behavior is not a major problem if θ is restricted to the range [0,π2]. Moreover, G^prx is also unphysical as θ ↓ 0 when (n − 2)t > 3; this condition on dimension and time is obtained by expanding e-θ2/4t=1-θ24t+O(θ4) and (sinθ/θ)-n-22=1+θ212(n-2)+O(θ3), and noting that the leading order θ² term in the product of the two factors is a non-decreasing function of distance θ when n-212≥14t, corresponding to the unphysical situation of nearby points being hotter than the heat source itself. As the feature dimension n is typically very large, the restriction (n − 2)t < 3 implies that we need to take the diffusion time to be very small, thus making the similarity measure captured by G^prx decay too fast away from each data point for use in clustering applications. In this work, we further computed the first and second order correction terms, denoted u₁ and u₂ in Equations (5), (6), respectively. In high dimensions, the divergence of u₁/u₀ and u₂/u₀ at θ = π is not a major problem, as we expect the expansion to be valid only in the vicinity θ ↓ 0; however, the divergence of u₂/u₀ at θ = 0 (to −∞ in high dimensions) is pathological, and thus, we truncate our approximation to O(t2). Since u₁(θ) is not able to correct the unphysical behavior of the parametrix kernel near θ = 0 in high dimensions, we conclude that the parametrix approximation fails in high dimensions. Hence, the only remaining part of G^prx still applicable to SVM classification is the Gaussian factor, which is clearly not a heat kernel on the hypersphere. The failure of this perturbative expansion using the Euclidean heat kernel as a starting point suggests that diffusion in ℝ^d and S^d are fundamentally different and that the exact hyperspherical heat kernel derived from a non-perturbative approach will likely yield better insights into the diffusion process.

By definition, the exact heat kernel Gext(x^,y^;t) is the fundamental solution to heat equation ∂tu+L^2u=0 where -L^2 is the hyperspherical Laplacian [13, 14, 19, 20]. In the language of operator theory, Gext(x^,y^;t) is an integral kernel, or Green's function, for the operator exp{-L^2t} and has an associated eigenfunction expansion. Because L^2 and exp{-L^2t} share the same eigenfunctions, obtaining the eigenfunction expansion of Gext(x^,y^;t) amounts to solving for the complete basis of eigenfunctions of L^2. The spectral decomposition of the Laplacian is in turn facilitated by embedding S^{n − 1} in ℝⁿ and utilizing the global rotational symmetry of Sⁿ⁻¹ in ℝⁿ. The Euclidean space harmonic functions, which are the solutions to the Laplace equation ∇²u = 0 in ℝⁿ, can be projected to the unit hypersphere Sⁿ⁻¹ through the usual separation of radial and angular variables [17, 18]. In this formalism, the hyperspherical Laplacian -L^2 on Sⁿ⁻¹ naturally arises as the angular part of the Euclidean Laplacian on ℝⁿ, and L^2 can be interpreted as the squared angular momentum operator in ℝⁿ [18].

The resulting eigenfunctions of L^2 are known as the hyperspherical harmonics and generalize the usual spherical harmonics in ℝ³ to higher dimensions. Each hyperspherical harmonic is equipped with a triplet of parameters or “quantum numbers” (ℓ, {m_i}, α): the degree ℓ, magnetic quantum numbers {m_i} and α=n2-1. In the eigenfunction expansion of exp{-L^2t}, we use the addition theorem of hyperspherical harmonics to sum over the magnetic quantum number {m_i} and obtain the following result:

Theorem 1. The exact hyperspherical heat kernel Gext(x^,y^;t) can be expanded as a uniformly and absolutely convergent power series

in the interval x^·y^∈[-1,1] and for t > 0, where Cℓα(w) are the Gegenbauer polynomials and ASn-1=2πn2Γ(n2) is the surface area of Sⁿ⁻¹. Since the kernel depends on x^ and ŷ only through x^·y^, we will write Gext(x^,y^;t)=Gext(x^·y^;t).

A proof of similar expansions can be found in standard textbooks on orthogonal polynomials and geometric analysis, e.g., [17]. For the sake of completeness, we include a proof in Supplementary Material section 2.5.3.

Note that the exact kernel G^ext is a Mercer kernel re-expressed by summing over the degenerate eigenstates indexed by {m}. As before, we will rescale the kernel by self-similarity and define:

Definition 4. The exact kernel Kext(x^,y^;t) is the exact heat kernel normalized by self-similarity:

which is defined for t > 0, is non-negative, and attains global maximum 1 at x^·y^=1.

Note that unlike Kprx(x^,y^;t), Kext(x^,y^;t) explicitly depends on the feature dimension n. In general, SVM kernel hyperparameter tuning can be computationally costly for a data set with both high feature dimension and large sample size. In particular, choosing an appropriate diffusion time scale is an important challenge. On the one hand, choosing a very large value of t will make the series converge rapidly; but, then, all points will become uniformly similar, and the kernel will not be very useful. On the other hand, a too small value of t will make most data pairs too dissimilar, again limiting the applicability of the kernel. In practice, we thus need a guideline for a finite time scale at which the degree of “self-relatedness” is not singular, but still larger than the “relatedness” averaged over the whole hypersphere. Examining the asymptotic behavior of the exact heat kernel in high feature dimension n shows that an appropriate time scale is t~O(logn/n); in this regime the numerical sum in Theorem 1 satisfies a stopping condition at low orders in ℓ and the sample points are in moderate diffusion proximity to each other so that they can be accurately classified (Supplementary Material, section 2.5.4).

Figure 1A illustrates the diffusion process captured by our exact kernel Kext(x^,y^;t) in three feature dimensions at time t = t^*log3/3, for t^* = 0.5, 1.0, 2.0. In Figure 1B, we systematically compared the behavior of (1) dimension-independent parametrix kernel K^prx at time t = 0.5, 1.0, 2.0 and (2) exact kernel K^ext on Sⁿ⁻¹ at t = t^*logn/n for t^* = 0.5, 1.0, 2.0 and n = 3, 100, 200. By symmetry, the slope of K^ext vanished at the south pole θ = π for any time t and dimension n. In sharp contrast, K^prx had a negative slope at θ = π, again highlighting a singular behavior of the parametrix kernel. The “relatedness” measured by K^ext at the sweet spot t = logn/n was finite over the whole hypersphere with sufficient contrast between nearby and far away points. Moreover, the characteristic behavior of K^ext at t = logn/n did not change significantly for different values of the feature dimension n, confirming that the optimal t for many classification applications will likely reside near the “sweet spot” t = logn/n.

Linear SVM seeks a separating hyperplane that maximizes the margin, i.e., the distance to the nearest data point. The primal formulation of SVM attempts to minimize the norm of the weight vector w that is normal to the separating hyperplane, subject to either hard or soft margin constraints. In the so-called Lagrange dual formulation of SVM, one applies the Representer Theorem to rewrite the weight as a linear combination of data points; in this set-up, the dot products of data points naturally appear, and kernel SVM replaces the dot product operation with a chosen kernel evaluation. The ultimate hope is that the data points will become linearly separable in the new feature space implicitly defined by the kernel.

We evaluated the performance of kernel SVM using the

on two independent data sets: (1) WebKB data of websites from four universities (WebKB-4-University) [21], and (2) glioblastoma multiforme (GBM) mutation data from The Cancer Genome Atlas (TCGA) with 5-fold cross-validations (CV) (Supplementary Material, section 1). The WebKB-4-University data contained 4,199 documents in total comprising four classes: student (1,641), faculty (1,124), course (930), and project (504); in our analysis, however, we selected an equal number of representative samples from each class, so that the training and testing sets had balanced classes. Table 1 shows the average optimal prediction accuracy scores of the five kernels for a varying number of representative samples, using 393 most frequent word features (Supplementary Material, section 1). The exact kernel outperformed the Gaussian RBF and parametrix kernel, reducing the error by 41 ~ 45% and by 1 ~ 7%, respectively. Changing the feature dimension did not affect the performance much (Table 2).

In the TCGA-GBM data, there were 497 samples, and we aimed to impute the mutation status of one gene—i.e., mutant or wild-type—from the mutation counts of other genes. For each imputation target, we first counted the number m_r of mutant samples and then selected an equal number of wild-type samples for 5-fold CV. Imputation tests were performed for top 102 imputable genes (Supplementary Material, section 1). Table 3 shows the average prediction accuracy scores for 5 biologically interesting genes known to be important for cancer [22]:

For the remaining genes, the exact kernel generally outperformed the linear, cosine and parametrix kernels (Figure 2). However, even though the exact kernel dramatically outperformed the Gaussian RBF in the WebKB-4-University classification problem, the advantage of the exact kernel in this mutation analysis was not evident (Figure 2). It is possible that the radial degree of freedom ∑i=1nxi in this case, corresponding to the genome-wide mutation load in each sample, contained important covariate information not captured by the hyperspherical heat kernel. The difference in accuracy between the hyperspherical kernels (cos, prx, and ext) and the Euclidean kernels (lin and rbf) also hinted some weak dependence on class size m_r (Figure 2), or equivalently the sample size m = 2m_r. In fact, the level of accuracy showed much stronger correlation with the “effective sample size” m~ related to the empirical Vapnik-Chervonenkis (VC) dimension [4, 7, 36–38] of a kernel SVM classifier (Figures 3A–E); moreover, the advantage of the exact kernel over the Guassian RBF kernel grew with the effective sample size ratio m~cos/m~lin (Figure 3F, Supplementary Material, section 2.5.5).

By construction, our definition of the hyperspherical map exploits only the positive portion of the whole hypersphere, where the parametrix and exact heat kernels seem to have similar performances. However, if we allow the data set to assume negative values, i.e., the feature space is the usual ℝⁿ\{0} instead of ℝ≥0n\{0}, then we may apply the usual projective map, where each vector in the Euclidean space is normalized by its L²-norm. As shown in Figure 1B, the parametrix kernel is singular at θ = π and qualitatively deviates from the exact kernel for large values of θ. Thus, when data points populate the whole hypersphere, we expect to find more significant differences in performance between the exact and parametrix kernels. For example, Table 4 shows the kernel SVM classifications of 91 S&P500 Financials stocks against 64 Information Technology stocks (m = 155) using their log-return instances between January 5, 2015 and November 18, 2016 as features. As long as the number of features was greater than sample size, n > m, the exact kernel outperformed all other kernels and reduced the error of Gaussian RBF by 29 ~ 51% and that of parametrix kernel by 17 ~ 51%.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.