Unsupervised Learning of Density Estimates with Topological Optimization

Kernel density estimation is a key component of a wide variety of algorithms in machine learning, Bayesian inference, stochastic dynamics and signal processing. However, the unsupervised density estimation technique requires tuning a crucial hyperparameter: the kernel bandwidth. The choice of bandwidth is critical as it controls the bias-variance trade-off by over- or under-smoothing the topological features. Topological data analysis provides methods to mathematically quantify topological characteristics, such as connected components, loops, voids et cetera, even in high dimensions where visualization of density estimates is impossible. In this paper, we propose an unsupervised learning approach using a topology-based loss function for the automated and unsupervised selection of the optimal bandwidth and benchmark it against classical techniques – demonstrating its potential across different dimensions.

💡 Research Summary

**
The paper addresses a long‑standing problem in kernel density estimation (KDE): how to choose the bandwidth h without supervision. Classical approaches (rules of thumb, plug‑in estimators, cross‑validation) optimize statistical criteria such as mean integrated squared error (MISE) or likelihood. While these criteria guarantee good global accuracy, they ignore the geometric and topological shape of the estimated density. Small changes in h can cause level‑sets to merge or split, dramatically altering the number of connected components, loops, or higher‑dimensional voids. Recent theoretical work has shown that an MISE‑optimal bandwidth may completely misrepresent the underlying topology.

To fill this gap, the authors propose a fully unsupervised method that selects h by minimizing a loss function derived from persistent homology, the main tool of topological data analysis (TDA). The procedure is as follows:

1. KDE Construction – Given data X⊂ℝᵈ, compute the Gaussian KDE \hat f_h(x)=1/(n hᵈ)∑K((x−x_i)/h). The kernel can be any smooth function, but the experiments use a Gaussian.

2. Cubical Super‑Level Filtration – Normalize \hat f_h, then for each threshold a construct the super‑level set X_a={x | \hat f_h(x)≥a}. This set is represented as a cubical complex, which is computationally convenient for image‑like data and for high‑dimensional grids.

3. Persistent Homology – Compute the persistence diagram (PD) of the filtration. Each point (b,d) in the PD corresponds to a topological feature (connected component, loop, void, etc.) that appears at birth b and disappears at death d. The lifetime ℓ=d−b measures its significance.

4. Topological Loss – Two quantities are extracted from the PD:

   • Feature Count = ∑σ(ℓ_i), where σ is a sigmoid (σ(t)=1/(1+e^{−t})). This smooth surrogate approximates the discrete Betti numbers, rewarding a moderate number of persistent features.
   • Total Persistence = ∑ℓ_i, which penalizes many short‑lived features that typically arise from over‑fitting. The loss is defined as
     L(h)=α_count·count(h) − α_TP·TP(h).
     In all experiments α_count=α_TP=1, so the loss simply balances simplicity (few features) against complexity (many short‑lived features).

5. Optimization – Because the loss is differentiable (thanks to the sigmoid), automatic differentiation can be applied. The authors use stochastic gradient descent (SGD) to minimize L(h) with respect to h. Persistent homology is stable under small perturbations, so L(h) often exhibits flat plateaus where many h values give identical topological structure. SGD drifts within these plateaus toward a value that also respects the smoothness of the loss surface.

6. Experiments – The method (referred to as “TDA‑based”) is benchmarked against a wide suite of classical bandwidth selectors:

   • Rules of thumb: Scott, Silverman, Normal‑Reference‑Rule (NRR)
   • Plug‑in: Improved Sheather‑Jones (ISJ, 1‑D only), Botev’s FFT‑based estimator, diagonal plug‑in
   • Cross‑validation: maximum‑likelihood CV (ML‑CV), least‑squares CV (LSCV), biased CV (BCV)
   • Adaptive/other: BotevProj, PluginDiag Synthetic data include 1‑D bimodal Gaussian, a mixed Gaussian‑Cauchy distribution, 2‑D clustered mixtures, and elongated elliptical mixtures, as well as a 4‑D synthetic example and MNIST digit images (784‑dimensional). Performance is measured by Kullback‑Leibler divergence (KLD) and Earth Mover’s Distance (EMD).

   Results – Across 500 repetitions per setting, the TDA‑based selector consistently yields low KLD and competitive EMD. In 1‑D, it achieves the lowest KLD on the bimodal case and comparable EMD to the best cross‑validation methods, while showing far less variance than LSCV or BCV. For the more complex mixed distribution, it outperforms most rule‑of‑thumb and plug‑in methods, though ISJ (a 1‑D specialist) still edges it in KLD. In 2‑D, TDA‑based bandwidths capture cluster topology better than overly smooth selectors (Silverman, Scott) and are far more stable than LSCV/BCV, which often diverge. On MNIST, the method produces sensible bandwidths that preserve digit‑wise cluster structure without manual tuning.

7. Discussion – The authors highlight three main advantages:

   • Topology‑aware smoothing – By directly penalizing undesirable topological changes, the estimator avoids both over‑smoothing (loss of genuine modes) and over‑fitting (spurious loops).
   • No hyper‑parameter tuning – The loss uses a single sigmoid with fixed scale; α‑weights are set to 1, eliminating the need for significance thresholds or bootstrapping.
   • Fully unsupervised – The approach requires only the data and a kernel; no labeled validation set or external criteria are needed.

   Limitations include reliance on Gaussian kernels in the experiments, computational cost of building cubical complexes in very high dimensions, and the fixed α‑weights that may not be optimal for all data regimes.

8. Future Work – Suggested extensions are: (i) supporting other kernels and anisotropic bandwidth matrices, (ii) developing more scalable persistent homology pipelines for dimensions > 4, (iii) enriching the loss with additional topological statistics (e.g., persistence entropy, moments of lifetimes), and (iv) applying the method to downstream tasks such as non‑parametric Bayesian inference, clustering, or anomaly detection.

Overall Assessment – The paper introduces a novel, mathematically grounded framework that bridges KDE and topological data analysis. By turning persistent homology into a differentiable loss, it provides a principled way to select bandwidths that respect the shape of the underlying distribution. Empirical results demonstrate that the method is competitive with, and often superior to, a broad set of classical selectors, especially in multimodal or heterogeneous settings where preserving topological features is crucial. The work opens a promising research direction at the intersection of statistics, geometry, and machine learning.