polykde: Kernel Density
Estimation with Polyspherical Data
Eduardo García-Portugués and Andrea Meilán-Vila
2026-11-13, v1.3.0
The polysphere \(\mathbb{S}^{d_1,\ldots,d_r} := \mathbb{S}^{d_1}
\times \cdots \times \mathbb{S}^{d_r}\) is the Cartesian product
of \(r\) spheres \(\mathbb{S}^{d_j} := \{\boldsymbol{x} \in
\mathbb{R}^{d_j + 1} : \|\boldsymbol{x}\| = 1\}\), and it hosts
multivariate directional data such as sets of angles or unit 3D vectors,
like the spokes of a skeletal shape representation. The
polykde package implements the kernel density estimator on
the polysphere, together with the bandwidth selectors, samplers,
regression estimators, ridge estimators, and homogeneity tests studied
in García-Portugués and Meilán-Vila (2026)
and García-Portugués and Meilán-Vila
(2023). This vignette showcases this functionality.
Polyspherical data
A sample of \(n\) points on \(\mathbb{S}^{d_1,\ldots,d_r}\) is simply
stored as a matrix of size
c(n, sum(d) + r) whose rows are the \(r\) unit-norm blocks concatenated by
columns. The dimensions are given by the vector d of length
r, so a point on \(\mathbb{S}^1
\times \mathbb{S}^2\) occupies \((1 +
1) + (2 + 1) = 5\) columns.
d <- c(1, 2)
r <- length(d)
X <- r_unif_polysph(n = 4, d = d)
X
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.99999969 0.0007828376 0.62165837 -0.1236640 -0.7734650
#> [2,] 0.29855291 0.9543930845 0.01039503 0.5498818 -0.8351778
#> [3,] -0.82064279 -0.5714415247 -0.33046061 0.3427514 0.8793846
#> [4,] 0.06032532 -0.9981787692 0.18790139 0.3420383 -0.9207078
dim(X)
#> [1] 4 5
# The 0-based indexes delimiting the blocks: columns 1:2 and 3:5
comp_ind_dj(d = d)
#> [1] 0 2 5
# Each block has unit norm
rowSums(X[, 1:2]^2)
#> [1] 1 1 1 1
rowSums(X[, 3:5]^2)
#> [1] 1 1 1 1
Cartesian coordinates can be converted to and from angles. The
generic pair is angles_to_polysph() /
polysph_to_angles(); the shortcuts
angles_to_sph() / sph_to_angles() handle a
single sphere and angles_to_torus() /
torus_to_angles() the torus \((\mathbb{S}^1)^r\).
On \(\mathbb{S}^2\) the first angle
is the polar angle measured from \((1, 0,
0)^\top\) and the second is the azimuth, so that \((\cos\theta_1, \sin\theta_1\cos\theta_2,
\sin\theta_1\sin\theta_2)^\top\).
theta <- rbind(c(pi / 2, pi / 4), c(pi / 3, -pi / 2))
x <- angles_to_sph(theta = theta)
x
#> [,1] [,2] [,3]
#> [1,] 6.123234e-17 7.071068e-01 0.7071068
#> [2,] 5.000000e-01 5.302876e-17 -0.8660254
sph_to_angles(x = x)
#> [,1] [,2]
#> [1,] 1.570796 0.7853982
#> [2,] 1.047198 4.7123890
Simulating polyspherical data
The r_*_polysph() samplers draw from the uniform
distribution, a product of von Mises–Fisher (vMF) distributions, and
mixtures of the latter. Their d_* counterparts evaluate the
corresponding densities with respect to the Lebesgue measure on the
polysphere.
set.seed(42)
d <- c(2, 2)
mu <- c(0, 0, 1, 1, 0, 0)
kappa <- c(5, 10)
X <- r_vmf_polysph(n = 500, d = d, mu = mu, kappa = kappa)
# Densities at the modal point
d_vmf_polysph(x = rbind(mu), d = d, mu = mu, kappa = kappa)
#> [,1]
#> [1,] 1.266572
d_unif_polysph(x = rbind(mu), d = d)
#> [1] 0.006332574
Mixtures take a matrix mu of modes (one per row), a
matrix kappa of concentrations, and the vector
prop of mixing proportions.
mu_mix <- rbind(c(0, 0, 1, 1, 0, 0), c(1, 0, 0, 0, 0, 1))
kappa_mix <- rbind(c(10, 10), c(20, 5))
prop <- c(0.7, 0.3)
X_mix <- r_mvmf_polysph(n = 500, d = d, mu = mu_mix, kappa = kappa_mix,
prop = prop)
head(d_mvmf_polysph(x = X_mix, d = d, mu = mu_mix, kappa = kappa_mix,
prop = prop))
#> [,1]
#> [1,] 0.87912608
#> [2,] 0.14626730
#> [3,] 0.53362937
#> [4,] 0.61843443
#> [5,] 0.03989348
#> [6,] 0.26016942
r_path_s1r() and r_path_s2r() generate
samples concentrated about one-dimensional paths on \((\mathbb{S}^1)^r\) and \((\mathbb{S}^2)^r\), which is convenient for
illustrating ridge and regression estimation.
Y <- r_path_s2r(n = 200, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
dim(Y)
#> [1] 200 3
Kernels
Three kernels are available through the kernel argument:
1 for von Mises–Fisher, 2 for Epanechnikov,
and 3 for softplus (with shape parameter k).
Each can be used as a product kernel over the \(r\) blocks (kernel_type = 1 or
"prod") or as a spherically symmetric kernel on
the whole polysphere (kernel_type = 2 or
"sph"). The function L() evaluates the kernel
profile and c_kern() its normalizing constant.
t_grid <- seq(0, 2, l = 200)
plot(t_grid, L(t = t_grid, kernel = "1"), type = "l", ylim = c(0, 1), lwd = 2,
xlab = "t", ylab = "L(t)")
lines(t_grid, L(t = t_grid, kernel = "2"), col = 2, lwd = 2)
lines(t_grid, L(t = t_grid, kernel = "3", k = 10), col = 3, lwd = 2)
legend("topright", legend = c("von Mises-Fisher", "Epanechnikov", "softplus"),
lwd = 2, col = 1:3)
# Normalizing constant of the product vMF kernel on (S^2)^2
c_kern(h = c(0.2, 0.2), d = c(2, 2), kernel = 1, kernel_type = 1)
#> [1] 3.978874 3.978874
The efficiency of a kernel relative to a reference (the spherically
symmetric Epanechnikov kernel by default) follows from the moments
b_d() and v_d().
eff_kern(d = 2, r = 2, kernel = "1", kernel_type = "sph")
#> [1] 0.75
eff_kern(d = 2, r = 2, kernel = "3", k = 10, kernel_type = "sph")
#> [1] 0.9847307
Kernel density estimation
kde_polysph() evaluates \(\hat{f}(\boldsymbol{x};\boldsymbol{h}) =
\sum_{i=1}^n L_{\boldsymbol{h}}(\boldsymbol{x},
\boldsymbol{X}_i)\) at the rows of x. The sample
X, the dimensions d, and the bandwidths
h (one per sphere) are the only required inputs. Setting
log = TRUE returns the log-density and
wrt_unif = TRUE expresses it with respect to the uniform
measure instead of the Lebesgue measure.
set.seed(42)
d <- 2
X <- r_mvmf_polysph(n = 500, d = d, mu = rbind(c(0, 0, 1),
c(1, 1, 0) / sqrt(2)),
kappa = rbind(10, 25), prop = c(0.6, 0.4))
h <- bw_rot_polysph(X = X, d = d)$bw
# Evaluation on a regular grid of spherical angles, with theta_1 the polar
# angle and theta_2 the azimuth (see sph_to_angles())
th_1 <- seq(0, pi, l = 100)
th_2 <- seq(0, 2 * pi, l = 200)
grid <- angles_to_sph(theta = as.matrix(expand.grid(th_1, th_2)))
kde <- matrix(kde_polysph(x = grid, X = X, d = d, h = h), nrow = length(th_1))
image(th_1, th_2, kde, col = hcl.colors(50, "Viridis"), useRaster = TRUE,
xlab = expression(theta[1]), ylab = expression(theta[2]))
contour(th_1, th_2, kde, add = TRUE, col = "white", drawlabels = FALSE)
points(sph_to_angles(x = X), pch = 16, cex = 0.4)
Two companions of the estimator are r_kde_polysph(),
which samples from the fitted density, and
log_cv_kde_polysph(), which returns the
leave-one-out log-density at each sample point. The latter is a
useful notion of depth: large values flag central observations and small
ones flag outliers.
head(r_kde_polysph(n = 3, X = X, d = d, h = h))
#> [,1] [,2] [,3]
#> [1,] -0.30324662 -0.09621508 0.9480423
#> [2,] 0.50173842 -0.01603889 0.8648707
#> [3,] 0.02832985 -0.24133180 0.9700291
log_dens <- log_cv_kde_polysph(X = X, d = d, h = h)
summary(drop(log_dens))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -7.0554 -1.1518 -0.8069 -0.9739 -0.5611 -0.3895
Bandwidth selection
Four selectors are available, in increasing order of computational
cost:
bw_mrot_polysph(): a marginal rule of thumb, one sphere
at a time.bw_rot_polysph(): a rule of thumb minimizing the AMISE
under a product of vMF reference densities.bw_cv_polysph(type = "LCV"): likelihood
cross-validation.bw_cv_polysph(type = "LSCV"): least squares
cross-validation.
set.seed(423432)
d <- c(1, 2)
mu <- c(0, 1, 0, 0, 1)
kappa <- c(5, 10)
X <- r_vmf_polysph(n = 200, d = d, mu = mu, kappa = kappa)
bw_mrot_polysph(X = X, d = d)
#> [1] 0.1644680 0.1356396
bw_rot_polysph(X = X, d = d)$bw
#> [1] 0.2033518 0.1497927
bw_cv_polysph(X = X, d = d, type = "LCV")$bw
#> [1] 0.2309878 0.1874455
The LSCV loss is a Monte Carlo approximation by default, but a closed
form is available for the vMF kernel through
exact_vmf = TRUE. common_h = TRUE forces a
single bandwidth across all spheres, which is the sensible choice when
r is large and the spheres play symmetric roles, and
arcsinh = TRUE stabilizes the loss numerically in high
dimensions.
bw_cv_polysph(X = X, d = d, type = "LSCV", exact_vmf = TRUE)$bw
#> [1] 0.1240745 0.1674424
bw_cv_polysph(X = X, d = d, type = "LSCV", exact_vmf = TRUE,
common_h = TRUE)$bw
#> [1] 0.1566476
Kernel regression
kre_polysph() regresses a polyspherical response
Y on a real covariate X by local constant
(p = 0, Nadaraya–Watson) or local linear
(p = 1) fitting, projecting the fit back onto the
polysphere. bw_cv_kre_polysph() selects the bandwidth by
cross-validation, returning both the minimizer h_min and
the more conservative one-standard-error rule h_1se.
set.seed(423432)
n <- 100
x <- seq(0, 1, l = n)
Y <- r_path_s2r(n = n, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
h_kre <- bw_cv_kre_polysph(X = x, Y = Y, d = 2, p = 0, plot_cv = FALSE)$h_1se
h_kre
#> [1] 0.09800103
x_grid <- seq(-0.25, 1.25, l = 200)
sc3 <- scatterplot3d::scatterplot3d(Y, pch = 16, xlim = c(-1, 1),
ylim = c(-1, 1), zlim = c(-1, 1),
xlab = "", ylab = "", zlab = "",
color = gray(0.5))
sc3$points3d(kre_polysph(x = x_grid, X = x, Y = Y, d = 2, h = h_kre, p = 0),
type = "l", col = 2, lwd = 2)
sc3$points3d(kre_polysph(x = x_grid, X = x, Y = Y, d = 2, h = h_kre, p = 1),
type = "l", col = 3, lwd = 2, lty = 2)
legend("topleft", legend = c("Local constant", "Local linear"), lwd = 2,
lty = 1:2, col = 2:3)
Density ridges
The density ridge is the one-dimensional set of points where the
density is at a local maximum in the direction of least curvature (Ozertem and Erdogmus 2011; Genovese et al.
2014). euler_ridge() runs the Euler algorithm of
García-Portugués and Meilán-Vila (2023)
from a set of starting points x, using the ridge bandwidth
h and the step size h_euler. The gradient and
Hessian it relies upon are exported as
grad_hess_kde_polysph() and
proj_grad_kde_polysph().
set.seed(987204452)
n <- 100
X <- r_path_s2r(n = n, r = 1, spiral = FALSE, Theta = cbind(c(1, 0, 0)),
sigma = 0.35)[, , 1]
rid <- euler_ridge(x = X, X = X, d = 2, h = 0.5, h_euler = 0.25, N = 30,
keep_paths = TRUE, show_prog = FALSE)
# The ridge points, the iterations spent, and the convergence flags
head(rid$ridge_y)
#> [,1] [,2] [,3]
#> [1,] 0.004171400 0.9838102 0.17916485
#> [2,] 0.025532579 0.9986988 0.04414606
#> [3,] -0.001888754 0.9739162 0.22689998
#> [4,] -0.002713002 0.9722556 0.23390526
#> [5,] -0.005584929 0.9658791 0.25893307
#> [6,] -0.010970432 0.9504782 0.31059765
table(rid$iter)
#>
#> 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
#> 1 2 3 5 8 7 12 9 15 7 4 5 9 2 2 5 2 1 1
table(rid$conv)
#>
#> 1
#> 100
# Starting points, ridge points, and the paths joining them
sc3 <- scatterplot3d::scatterplot3d(X, pch = 16, xlim = c(-1, 1),
ylim = c(-1, 1), zlim = c(-1, 1),
xlab = "x", ylab = "y", zlab = "z",
color = gray(0.5))
for (k in seq_len(nrow(rid$paths))) {
sc3$points3d(t(rid$paths[k, , ]), type = "l", col = gray(0.75))
}
sc3$points3d(rid$ridge_y, pch = 16, cex = 0.75, col = 2)
clean_euler_ridge() removes non-convergent or spurious
endpoints, and index_ridge() parameterizes the resulting
curve, delivering scores that play the role of a first principal
component along the ridge.
rid_clean <- clean_euler_ridge(e = rid, X = X)
ind <- index_ridge(endpoints = rid_clean$ridge_y, X = X, d = 2,
probs_scores = seq(0, 1, l = 50))
summary(ind$scores_X)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -1.8121 -1.8121 -0.1706 0.0000 1.9234 1.9234
For a polysphere with many spheres,
parallel_euler_ridge() distributes the starting points
across cores and block_euler_ridge() runs the Euler
algorithm marginally on groups of spheres.
Homogeneity tests
hom_test_polysph() performs permutation tests for the
equality of the distributions of \(k\)
polyspherical samples. The omnibus statistic is the Jensen–Shannon
distance between the \(k\) kernel
density estimators (type = "jsd", the only one supporting
\(k > 2\)); "hd" uses
the Hellinger distance, while "mean" and
"scatter" compare only the group means and scatter
matrices.
set.seed(423432)
n <- c(100, 100)
X1 <- r_vmf_polysph(n = n[1], d = 2, mu = c(0, 0, 1), kappa = 5)
X2 <- r_vmf_polysph(n = n[2], d = 2, mu = c(0, 0, 1), kappa = 5)
# Homogeneity holds
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
type = "jsd", h = 0.5, B = 500, show_prog = FALSE)
#>
#> Permutation-based Jensen--Shannon distance test of homogeneity
#>
#> data: rbind(X1, X2)
#> Tn = -0.0041028, p-value = 0.6307
#> alternative hypothesis: any alternative to homogeneity
# It does not: the second sample is more concentrated
X2 <- r_vmf_polysph(n = n[2], d = 2, mu = c(0, 0, 1), kappa = 20)
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
type = "jsd", h = 0.5, B = 500, show_prog = FALSE)
#>
#> Permutation-based Jensen--Shannon distance test of homogeneity
#>
#> data: rbind(X1, X2)
#> Tn = 0.036579, p-value = 0.001996
#> alternative hypothesis: any alternative to homogeneity
# The mean-based test is blind to this alternative
hom_test_polysph(X = rbind(X1, X2), d = 2, labels = rep(1:2, times = n),
type = "mean", B = 500, show_prog = FALSE)
#>
#> Permutation-based test of means equality
#>
#> data: rbind(X1, X2)
#> Tn = 0.00058435, p-value = 0.6986
#> alternative hypothesis: both group means are different
Hippocampus shapes
The hippocampus dataset contains skeletal
representations (s-reps) of the hippocampi of \(177\) 6-month-old infants, \(34\) of whom later developed autism (Liu et al. 2021). Each s-rep has \(168\) spokes, so the spoke
directions of a subject form a point on \((\mathbb{S}^2)^{168}\). Converting the
dirs array into the c(n, sum(d) + r) layout is
a one-liner.
data("hippocampus")
dirs <- hippocampus$dirs
r <- ncol(dirs)
d <- rep(2, r)
X <- do.call(cbind, args = lapply(1:r, function(i) dirs[, i, ]))
dim(X)
#> [1] 177 504
Following García-Portugués and Meilán-Vila
(2026), we use the spherically symmetric softplus kernel and a
rule-of-thumb bandwidth, and then rank the subjects by their
leave-one-out log-density. The extremes are the most representative and
the most outlying hippocampi.
h_mrot <- bw_mrot_polysph(X = X, d = d, kernel = 3, k = 100)
h_rot <- bw_rot_polysph(X = X, d = d, bw0 = 2^(-1:2) %o% h_mrot,
kernel = 3, kernel_type = 2, k = 100)$bw
log_dens <- log_cv_kde_polysph(X = X, d = d, h = h_rot, wrt_unif = TRUE,
kernel = 3, kernel_type = 2, k = 100)
deep <- which.max(log_dens)
outly <- which.min(log_dens)
c("deepest" = deep, "most outlying" = outly)
#> deepest most outlying
#> 156 137
view_srep() displays an s-rep, statically via
scatterplot3d or interactively via rgl.
par(mfrow = c(1, 2))
titles <- c("Deepest", "Most outlying")
for (j in 1:2) {
i <- c(deep, outly)[j]
view_srep(base = hippocampus$base[i, , ], dirs = dirs[i, , ],
bdry = hippocampus$bdry[i, , ], static = TRUE, main = titles[j])
}
Finally, a homogeneity test compares the spoke directions of the two
diagnostic groups. A small number of permutations is used here to keep
the vignette light; García-Portugués and
Meilán-Vila (2026) uses B = 5e3.
set.seed(423432)
hom_test_polysph(X = X, d = d, labels = hippocampus$ids_labs, type = "jsd",
h = 4 * h_rot, kernel = 3, kernel_type = 2, k = 100,
B = 100, show_prog = FALSE)
#>
#> Permutation-based Jensen--Shannon distance test of homogeneity
#>
#> data: X
#> Tn = 0.0015581, p-value = 0.08911
#> alternative hypothesis: any alternative to homogeneity
The scripts in the /replication folder of the GitHub repository
reproduce the complete numerical experiments and data applications of
García-Portugués and Meilán-Vila (2026)
and its Supplementary Material.
References
García-Portugués, E., and A. Meilán-Vila. 2023.
“Hippocampus Shape
Analysis via Skeletal Models and Kernel Smoothing.” In
Statistical Methods at the Forefront of Biomedical Advances,
edited by Y. Larriba. Springer.
https://doi.org/10.1007/978-3-031-32729-2_4.
García-Portugués, E., and A. Meilán-Vila. 2026.
“Kernel Density
Estimation with Polyspherical Data and Its Applications.” J.
Amer. Statist. Assoc. 121 (553): 427–39.
https://doi.org/10.1080/01621459.2025.2521898.
Genovese, C. R., M. Perone-Pacifico, I. Verdinelli, and L. Wasserman.
2014.
“Nonparametric Ridge Estimation.” Ann.
Statist. 42 (4): 1511–45.
https://doi.org/10.1214/14-AOS1218.
Liu, Z., J. Hong, J. Vicory, J. N. Damon, and S. M. Pizer. 2021.
“Fitting Unbranching Skeletal Structures to Objects.”
Med. Image Anal. 70: 102020.
https://doi.org/10.1016/j.media.2021.102020.
Ozertem, U., and D. Erdogmus. 2011. “Locally Defined Principal
Curves and Surfaces.” J. Mach. Learn. Res. 12: 1249–86.