Subspace-Constrained Mean Shift (SCMS)

Updated 23 February 2026

SCMS is a gradient-based algorithm that detects density ridges—manifold-like structures with local maxima—in point cloud data.
It combines kernel density estimation with subspace-projected mean-shift updates to efficiently extract filamentary features from complex datasets.
The method incorporates bootstrap uncertainty quantification and adapts to diverse geometries, enhancing its application across astronomy, statistics, and signal processing.

The Subspace-Constrained Mean Shift (SCMS) algorithm is a gradient-based method designed to identify density ridges—manifold-like structures where a probability density function exhibits local maxima along lower-dimensional subspaces—within point cloud data. Originally motivated by cosmic web reconstruction, SCMS is now applied across statistics, astronomy, and signal processing to extract high-density filamentary features such as galaxy filaments, tidal streams, and ridges in generic high-dimensional data. SCMS leverages kernel density estimation (KDE), subspace-projection, and iterative fixed-point methods to provide a statistically principled, parameter-efficient, and uncertainty-quantified approach to nonparametric filament detection (Chen et al., 2015, Chen et al., 2015, Hendel et al., 2018).

1. Formal Definition of Density Ridges

Given points $X_1, \ldots, X_n \subset \mathbb{R}^d$ , a kernel density estimator of the underlying density is defined as

$p(x) = \frac{1}{n h^d} \sum_{i=1}^n K\left(\frac{x - X_i}{h}\right)$

where $K$ is a smooth, radial, symmetric kernel (commonly Gaussian), and $h > 0$ is the bandwidth. Denote the gradient $g(x) = \nabla p(x)$ and Hessian $H(x) = \nabla^2 p(x)$ , with eigenvalues $\lambda_1(x) \geq \cdots \geq \lambda_d(x)$ and corresponding orthonormal eigenvectors $v_1(x), \ldots, v_d(x)$ .

The one-dimensional density ridge, representing a filament, is defined as

$R = \left\{ x : v_j(x)^T g(x) = 0 \quad (j = 2, \ldots, d),\; \lambda_2(x) < 0 \right\}$

This expresses that along the ridge, the density gradient is entirely along $v_1(x)$ (the direction of maximal curvature), and the density curves downward in orthogonal directions (Chen et al., 2015, Chen et al., 2015).

2. The SCMS Update Step

The classical mean-shift vector at $x$ is given by

$m(x) = \frac{\sum_{i=1}^n K\left(\frac{x - X_i}{h}\right) (X_i - x)}{\sum_{i=1}^n K\left(\frac{x - X_i}{h}\right)}$

For a Gaussian kernel, $m(x) = \frac{h^2}{p(x)} g(x)$ . The SCMS algorithm constrains movement to the ridge-attracting subspace orthogonal to the leading eigenvector. Defining $V(x) = [v_2(x)\ \cdots\ v_d(x)]$ (a $d \times (d-1)$ matrix of minor eigendirections), the SCMS update is

$x_{\text{new}} = x + V(x)V(x)^T m(x)$

or, using the gradient,

$x_{\text{new}} = x + V(x)V(x)^T \frac{h^2 g(x)}{p(x)}$

Iterating this update causes points to ascend to the density ridge, constrained within the subspace orthogonal to the ridge direction (Chen et al., 2015, Hendel et al., 2018).

3. Complete SCMS Algorithm

The SCMS pipeline, as typically implemented, proceeds in these stages:

Density Estimation: Choose bandwidth $h$ (see below), compute $p(x)$ , $g(x)$ , and $H(x)$ at mesh points over the domain.
Thresholding: Calculate the root-mean-square (RMS) of $p(x)$ . Discard any $x$ with $p(x) < \mathrm{RMS}(p)$ to suppress spurious ridges in low-density regions.
Ridge Ascent: Initialize a grid (or use data points) as seeds. For each $x_k^{(0)}$ , iterate

$x^{(t+1)} = x^{(t)} + V(x^{(t)})V(x^{(t)})^T m(x^{(t)})$

until the projected mean-shift norm falls below a tolerance ( $\|V V^T m(x)\| < \varepsilon$ ) or a maximum step count is reached. The converged points $\{x_k^{(\infty)}\}$ form an approximation to the ridge $R$ (Chen et al., 2015, Chen et al., 2015).

Extensions to spheres and product manifolds involve adapting the KDE, gradients, and projection operators to non-Euclidean geometry, allowing SCMS to be applied on domains such as $\mathbb{S}^2$ and $\mathbb{S}^2 \times \mathbb{R}$ (Zhang et al., 2021, Zhang et al., 2022, Zhang et al., 2021).

4. Choice of Smoothing Bandwidth and Density Estimation

Bandwith selection critically affects filament geometry. The standard reference rule ("Silverman's rule") is

$h \propto n^{-1/(d+4)} \widehat{\sigma}$

with $\widehat{\sigma}$ the empirical standard deviation of the data. For cosmological applications, the increasing redshift reduces galaxy density, so $h$ is adapted per slice, ranging, for example, from $\sim 5^\circ$ at low $z$ to $\sim 15^\circ$ at high $z$ (Chen et al., 2015). On directional or mixed domains, analogous rules-of-thumb based on estimated concentration or marginal variance are used (Zhang et al., 2021, Zhang et al., 2022).

5. Uncertainty Quantification via Bootstrap

Uncertainty in the detected ridge is quantified by resampling:

Bootstrap Sampling: Draw $B$ resamples of the data with replacement.
Ridge Estimation: Apply the complete SCMS procedure separately to each bootstrap sample, yielding $R^{(1)}, \ldots, R^{(B)}$ .
Projection Distance: For each original ridge point $x \in R$ , compute distances $d_b(x) = \min_{y \in R^{(b)}} \|x - y\|$ over $b = 1,\ldots, B$ .
Summarization: Report uncertainty at each $x$ as the mean, quantiles, or RMS of $\{d_b(x)\}$ . Typical $B$ is $100$–$1000$ (Chen et al., 2015, Chen et al., 2015).

This yields a pointwise, data-driven uncertainty measure and enables the construction of geometry-adaptive uncertainty bands around each filament.

6. Implementation Parameters and Practical Considerations

Typical SCMS settings and steps for large-scale cosmic web mapping include:

Redshift Slicing: Data is partitioned into thin redshift bins (e.g., $\Delta z = 0.005$ ), with galaxies projected onto 2D angular coordinates per slice.
Spatial Window: The working area is restricted (example: RA $\in [150^\circ, 200^\circ]$ , Dec $\in [5^\circ, 30^\circ]$ ).
Thresholding: A density RMS threshold is enforced on all candidate points in each slice.
Seed Grid: The initial mesh is a uniform lattice, typically spaced at about $h/2$ .
Convergence: Iterations stop when $\|VV^T m(x)\| < 10^{-4}$ or a maximum iteration count (e.g., 200) is reached.
Intersection/Junction Detection: Each filament point is tested for intersection status by clustering neighboring points within an annulus; $x$ is flagged as a junction if at least three clusters are identified in its neighborhood (Chen et al., 2015).
Computational Considerations: Each Hessian computation is $O(n d^2)$ , eigen-decompositions $O(d^3)$ , so acceleration via spatial data structures and parallelization is common for large $n$ and $d$ (Hendel et al., 2018).
Parameter Sensitivity: Ridge extraction quality depends on $h$ , density threshold, and convergence tolerance. Robustness to these is a practical requirement for large astrophysical catalogues (Hendel et al., 2018).

7. Theoretical Properties, Convergence, and Extensions

While early SCMS lacked formal convergence proofs, recent developments establish SCMS as a specific instance of subspace-constrained gradient ascent (SCGA) with locally adaptive step sizes (Zhang et al., 2021, Zhang et al., 2021). Under mild regularity conditions (smoothness, eigengap, and path smoothness), SCMS exhibits local linear convergence: $\|x^t - x^*\| \leq \Upsilon^t \|x^0 - x^*\| + O(h^2) + O_P\left((n h^{d+4})^{-1/2}\right)$ for iterates $x^t$ initialized sufficiently close to a true ridge $x^*$ , with contraction rate $\Upsilon < 1$ (Zhang et al., 2021). Generalizations extend these guarantees to directional and product spaces, e.g., the sphere $\mathbb{S}^2$ and mixtures such as $\mathbb{S}^2 \times \mathbb{R}$ (Zhang et al., 2021, Zhang et al., 2022).

The ridge definition is stable to small perturbations in the density estimate, and consistency of the filament estimator in Hausdorff distance is achievable at the nonparametric minimax rate provided $n h^{d+4} \to \infty$ (Qiao et al., 2021, Zhang et al., 2021).