E FFICIENT M EAN C URVATURE C OMPUTATION ON H IGH -D IMENSIONAL DATA M ANIFOLDS
arXiv:2606.06329v1 [cs.LG] 4 Jun 2026
Alexandre Luis Magalhães Levada Federal University of São Carlos 13565-905, São Carlos-SP, Brazil [email protected]
June 5, 2026
A BSTRACT Estimating local mean curvature at each point of a high-dimensional dataset is a key ingredient of geometry-aware machine learning algorithms, such as the Mean Curvature Boundary Points (MCBP) method. The naive implementation of this computation, based on a local shape operator approximated from k-nearest neighbor patches, involves an explicit construction of a matrix H whose trace form yields an O(m4 ) cost per point, rendering the approach intractable for datasets with more than a few dozen features. This paper introduces two complementary contributions that together reduce this cost by several orders of magnitude. The first contribution is an exact algebraic identity. This identity, derived from the orthogonality of the eigenvectors of the covariance matrix and the cyclicity of the trace operator, eliminates H entirely and reduces the per-point cost to O(m2 ) after the eigendecomposition. The second contribution addresses the remaining O(m3 ) bottleneck of the full eigendecomposition. Since the local covariance matrix has rank at most k − 1 ≪ m, we replace it with a truncated SVD of the k × m centered data matrix, an O(k 2 m) operation, and derive an analytical approximation for the contribution of the null-space eigenvectors based on the expected value of their outer product under the Haar measure. The resulting estimator has total cost O(k 2 m + kmp2 ), where p = k − 1. Experiments on real-world datasets confirm speedups of 50 to 300 times relative to the original implementation, with negligible loss when the fast estimator is used to replace the original version. By providing a scalable and data-driven estimate of local curvature, the proposed method establishes curvature as a practical geometric feature for a broad range of machine learning tasks, from classical to modern deep learning pipelines.
1
Introduction
The last two decades have witnessed the emergence of Geometric Machine Learning (GML) as a unifying paradigm for data analysis, whose central tenet is that the geometric structure of the domain, rather than the vector-space structure of the ambient coordinates, should govern the design of representations, metrics, and algorithms Weber [2025], Papillon et al. [2025]. This perspective departs fundamentally from classical Euclidean-based machine learning, which treats every dataset as a collection of points in Rm and relies on the inner product structure of that space for tasks such as similarity search, dimensionality reduction, and classification. The limitations of this flat-geometry assumption became increasingly apparent as practitioners confronted data arising from social networks, brain-imaging graphs, molecular surfaces, protein interaction networks, and high-dimensional sensor arrays, domains in which the relevant notion of similarity or proximity is intrinsically non-Euclidean [Bronstein et al., 2017]. The seminal position paper by Bronstein et al. [2017] crystallised these ideas under the banner of geometric deep learning, arguing that the success of convolutional neural networks on images is itself a manifestation of a deeper symmetry principle, translational equivariance, and that a principled extension to graphs, manifolds, and other non-Euclidean domains requires replacing the Euclidean blueprint with one grounded in group theory and differential geometry. This programme was subsequently developed into a comprehensive theoretical framework by Bronstein et al. [2021], who unified convolutional networks, graph neural networks, transformers, and equivariant architectures under a single geometric blueprint inspired by
P REPRINT - J UNE 5, 2026
Klein’s Erlangen Programme, formalising how symmetry groups (grids, graphs, geodesics, and gauges) determine the admissible learning architectures for each class of structured data. Underlying this entire research agenda is the manifold hypothesis [Fefferman et al., 2016], which posits that highdimensional data encountered in practice are not spread uniformly through Rm but are instead concentrated near a smooth, low-dimensional Riemannian manifold M of intrinsic dimension d ≪ m embedded in the ambient space. This hypothesis, formalised with sample-complexity guarantees by Fefferman et al. [2016] and supported by a large body of empirical evidence, is the theoretical foundation of manifold learning — the family of methods that seek to recover the intrinsic geometry of M from finite, possibly noisy samples, without access to an explicit parametrisation of the manifold. Landmark algorithms such as ISOMAP [Tenenbaum et al., 2000], Locally Linear Embedding [Roweis and Saul, 2000], t-SNE van der Maaten and Hinton [2008] and UMAP McInnes et al. [2018a,b] demonstrated that global geodesic structure and local neighborhood geometry, respectively, can be reliably recovered from point-cloud data, thereby opening the door to a rich family of geometry-aware representations. Building on these foundations, GML moves beyond the goal of merely embedding data into a lower-dimensional Euclidean space and pursues instead the estimation of differential-geometric quantities, such as tangent spaces, curvature tensors, and shape operators, that characterise the local and global geometry of M directly from the data. Local mean curvature, in particular, occupies a privileged role in this programme: it quantifies the degree to which the manifold bends in the ambient space at each point, encoding information about cluster boundaries, concave and convex geometric features, and low-density transition regions that is invisible to density-based statistics alone [Asao and Ike, 2022, Cheng and Wu, 2021]. Efficient and scalable estimation of mean curvature from high-dimensional point clouds is therefore a problem of central importance in GML, yet existing approaches either assume low ambient dimension or incur computational costs that grow polynomially with m in ways that render them intractable for modern datasets. The present paper addresses this gap by deriving an exact algebraic identity that reduces the per-point cost of the mean curvature estimator introduced in Levada [2025] from O(m4 ) to O(m3 ), and by combining this identity with a truncated singular value decomposition to achieve a further reduction to O(k 2 m) for datasets in which the neighborhood size k is small relative to the ambient dimension m, precisely the regime that characterises high-dimensional machine learning applications. The present paper makes three interrelated scientific contributions that together constitute a significant step forward in the computational foundations of Geometric Machine Learning. First, we establish an exact closed-form identity that reformulates the pointwise mean curvature estimator of Levada [2025], originally defined through an explicit feature m matrix H ∈ Rm×(m+( 2 )) involving all element-wise squares and pairwise products of local eigenvectors, as a compact expression involving only the m × m matrix C = W ⊤ W (2) , where W is the orthogonal matrix of eigenvectors of the local covariance and W (2) is its element-wise square. This identity, whose derivation exploits the orthogonality of W to collapse a quartic tensor contraction into a single matrix multiplication followed by element-wise operations, reduces the dominant per-point computational cost from O(m4 ) to O(m3 ) with no approximation error. Second, we show that, in the high-dimensional regime where m ≫ k, the local covariance matrix has numerical rank at most p = k − 1, and we exploit this low-rank structure to replace the O(m3 ) full eigendecomposition with a truncated singular value decomposition of the k × m centred neighborhood matrix, at cost O(k 2 m). The contribution of the remaining m − p null-space eigenvectors is handled via an analytical approximation grounded in the expected outer product of random orthonormal bases under the Haar measure, yielding a total per-point cost of O(k 2 m + kmp2 ). Third, we demonstrate through systematic experiments that the resulting estimator achieves speedups of up to 800× over the original implementation for m = 200, while preserving the geometric fidelity required for downstream tasks, making high-dimensional mean curvature estimation, for the first time, computationally viable as a routine preprocessing step in modern machine learning pipelines. Beyond the MCBP algorithm that motivated this work, the efficient estimation of local mean curvature opens new possibilities across a broad spectrum of application domains in which the intrinsic geometry of the data manifold carries information that density-based or distance-based statistics cannot capture. In unsupervised clustering and boundary detection, high-curvature regions of the data manifold naturally identify transition zones between clusters, low-density interfaces, and geometric irregularities that mark the boundaries of data partitions [Levada, 2025]; incorporating curvature as a feature in algorithms such as DBSCAN or HDBSCAN can therefore sharpen cluster boundaries and reduce sensitivity to bandwidth parameters in settings with non-linear or heterogeneous cluster shapes. In anomaly and outlier detection, points at which the local manifold exhibits unusually high or inconsistent curvature are natural candidates for anomalous observations, since they correspond to regions where the data departs from the smooth low-dimensional structure expected under the manifold hypothesis [Fefferman et al., 2016]; this geometric perspective on anomaly scoring complements, and in many cases outperforms, purely density-based criteria, as recently
2
P REPRINT - J UNE 5, 2026
demonstrated for graph-structured data Grover et al. [2025] and multi-class unsupervised anomaly detection Guo et al. [2025]. In semi-supervised and active learning, identifying high-curvature boundary points provides a principled, geometrydriven criterion for selecting the most informative samples to label, since the classification uncertainty of discriminative models is highest precisely near the curved boundaries between classes [Levada, 2025]; this yields a curvature-aware query strategy that is complementary to conventional uncertainty-sampling or margin-based approaches. In singlecell genomics and bioinformatics, where datasets routinely have thousands of features (m ∼ 103 –104 ) but lie near manifolds of intrinsic dimension of order tens to hundreds, local curvature can reveal differentiation trajectories, identify transitional cell states near bifurcation points of the developmental manifold, and flag cells at the geometric boundary between annotated cell types, tasks for which current dimensionality reduction pipelines based on PCA, UMAP, or t-SNE provide only indirect and nonparametric evidence [Imoto et al., 2022]. In tabular data, autoencoder-based semi-supervised learning architectures have shown results comparable to state-of-the-art methods, particularly in scenarios with very limited labeled data Stevanoska et al. [2025], indicating that curvature-aware autoencoders and active learning methods Tharwat and Schenck [2023] are a promising natural direction of evolution. In geometric deep learning on graphs and point clouds, node-level or point-level curvature estimates provide expressive geometric descriptors that can be incorporated as input features or regularisation signals in graph neural networks, in a manner analogous to how discrete Ricci curvature has been used to detect community structure, improve message passing, and characterise anomalous nodes in attributed networks [Grover et al., 2025, Cheng and Wu, 2021]. The computational advances reported here, by making mean curvature estimation tractable for the ambient dimensions encountered in these domains, remove a key bottleneck that has so far prevented the wider adoption of curvature-aware methods in high-dimensional machine learning. The remaining of the paper is organized as follows: Section 2 presents the theoretical background on differential geometry. Section 3 discusses the original local mean curvature estimation algorithm, which is computationally unfeasible for large and high dimensional datasets. Section 4 describes the proposed algebraic identity for efficient local mean curvature computation and the resulting algorithm. Section 5 shows the computational experiments and the obtained results. Finally, Section 6 presents the conclusions and final remarks.
2
Differential Geometry Basics
The mathematical backbone of the proposed method is drawn from classical differential geometry, the branch of mathematics concerned with the infinitesimal structure of smooth curves, surfaces, and their higher-dimensional generalisations, Riemannian manifolds. Rather than treating a dataset as an unstructured collection of points in Rm , we adopt the perspective, now standard in Geometric Machine Learning, that observations are discrete samples from an underlying smooth manifold M embedded in the ambient space, and that the most informative geometric signal resides not in pairwise distances alone but in the way M curves and stretches relative to its embedding [Spivak, 1999, do Carmo, 2016, O’Neill, 2006, Tu, 2017]. To make this precise, we rely on five classical objects. The tangent space Tp M at a point p ∈ M is the best linear approximation to the manifold at p, encoding the directions along which one can move while remaining, to first order, on M. The first fundamental form, or metric tensor g, is the restriction of the ambient inner product to each tangent space; it governs intrinsic measurements (lengths, angles, and areas) that are independent of how M sits in Rm . The second fundamental form II captures the complementary, extrinsic information: it measures how much tangent vectors twist out of the tangent space as one moves along M, thereby quantifying the bending of the manifold in the ambient space [do Carmo, 2016, Oprea, 2007]. From these two forms one constructs the shape operator S, the self-adjoint linear map from Tp M to itself obtained by composing II with the inverse of g, whose eigenvalues are the principal curvatures κ1 , . . . , κd of M at p. Their arithmetic mean defines the Pd mean curvature H = d−1 i=1 κi , the scalar quantity that aggregates the net bending of the manifold at each point into a single, geometrically interpretable value [Spivak, 1999, Needham, 2021, O’Neill, 2006]. In the data-analytic context that motivates this paper, mean curvature serves as a principled, nonparametric indicator of geometric irregularity: regions of M where H is elevated correspond to sharp transitions between dense clusters, concave or convex geometric features, and low-density interfaces, precisely the boundary structures that classical density-based statistics tend to conflate with outliers or noise [Levada, 2025]. The remainder of this section formalises these objects in the notation used throughout the paper, with an emphasis on the discrete, sample-based estimators through which they become computationally accessible. 2.1
Tangent Spaces
The first step toward a rigorous treatment of curvature is to formalise the notion of direction on a manifold. Unlike points in Rm , whose directions are globally defined by the standard basis, a manifold M need not inherit a canonical 3
P REPRINT - J UNE 5, 2026
Figure 1: Illustration of the tangent space at a point x on a sphere embedded in R3 . The tangent plane Tx M provides a first-order linear approximation of the manifold in a local neighborhood of x, capturing the directions of admissible infinitesimal variations along the surface. linear structure from the ambient space. The tangent space resolves this difficulty by associating to each point p ∈ M a vector space that captures, to first order, the local geometry of M near p [do Carmo, 2016, Tu, 2017]. Intuitive description. Consider a smooth curve α : (−ε, ε) → M with α(0) = p. The velocity vector α′ (0) ∈ Rm points in a direction that is tangent to M at p. The collection of all such velocity vectors, as α ranges over all smooth curves through p, forms a d-dimensional linear subspace of Rm , the tangent space at p. Intuitively, Tp M is the best flat (affine) approximation to M at p: it is the unique d-dimensional hyperplane that osculates M to first order at that point [O’Neill, 2006, Needham, 2021]. Figure 1 illustrates the tangent plane to a sphere at a given point. Formal definition. Let M be a smooth d-dimensional manifold embedded in Rm , and let p ∈ M. The tangent space of M at p is defined as Tp M =
α′ (0) ∈ Rm α : (−ε, ε) → M smooth, α(0) = p .
(1)
One verifies that Tp M is closed under addition and scalar multiplication, and hence is indeed a vector space of dimension d [Tu, 2017]. Elements of Tp M are called tangent vectors at p. Coordinate representation. Let ϕ : U ⊂ Rd → M be a local parametrisation of M around p, with ϕ(u0 ) = p for some u0 ∈ U . The partial derivatives of ϕ at u0 , ∂i ≡
∂ϕ ∈ Rm , ∂ui u0
i = 1, . . . , d,
(2)
form a basis for Tp M, so that every tangent vector v ∈ Tp M can be written as
v =
d X
v i ∂i ,
v i ∈ R.
(3)
i=1
The matrix Jϕ = [∂1 | · · · | ∂d ] ∈ Rm×d is the Jacobian of the parametrisation; its column space is precisely Tp M [do Carmo, 2016, Oprea, 2007].
4
P REPRINT - J UNE 5, 2026
The tangent bundle.
Varying p over all of M yields the tangent bundle TM =
G
Tp M =
(p, v) p ∈ M, v ∈ Tp M ,
(4)
p∈M
which is itself a smooth manifold of dimension 2d. A vector field on M is a smooth section of T M, i.e., a smooth map X : M → T M with X(p) ∈ Tp M for every p. Vector fields serve as the natural domain of differential operators such as the covariant derivative and the Lie bracket, which in turn underpin the definitions of the curvature tensors introduced in subsequent subsections [Spivak, 1999, O’Neill, 2006]. Normal space.
The orthogonal complement of Tp M in Rm is the normal space at p, Np M =
n ∈ Rm ⟨n, v⟩ = 0 ∀ v ∈ Tp M ,
(5)
of dimension m − d. The ambient space decomposes orthogonally as Rm = Tp M ⊕ Np M, and the unit vectors in Np M are the normals to M at p. For a hypersurface (d = m − 1), the normal space is one-dimensional and spanned by a single unit normal n(p); for co-dimension greater than one, there exists a family of normal directions, and the extrinsic curvature depends on the chosen normal [do Carmo, 2016, Spivak, 1999]. Relevance in the data-analytic setting. In the context of machine learning, the tangent space plays a twofold role. First, it provides the local linear model underlying a wide family of geometry-aware algorithms: principal component analysis (PCA), locally linear embedding [Roweis and Saul, 2000], and diffusion maps all implicitly estimate Tp M through the leading eigenvectors of local covariance or graph-Laplacian operators. Second, and most directly relevant to the present work, the tangent space is the domain on which the shape operator S acts and on which the second fundamental form II is defined; it is therefore the indispensable geometric substrate for any curvature-based analysis. Given a local neighborhood N (p) of p in the dataset, we estimate Tp M as the span of the d leading eigenvectors of the sample covariance matrix Σp (defined in Section 4.1), a procedure whose consistency under the manifold hypothesis is established in Singer and Wu [2012] and whose sample complexity is controlled by the local reach and curvature of M. This eigenvector-based estimate of the tangent space is not only the starting point for our curvature estimator but also the object whose orthogonality, formally expressed as W⊤ W = Im , drives the key algebraic simplification derived in Section 4. 2.2
First Fundamental Form
Having established the tangent space Tp M as the natural linear approximation to the manifold at each point, the next step is to equip it with a notion of length and angle. This is the role of the first fundamental form, the foundational object of intrinsic Riemannian geometry: it encodes how distances and angles on M relate to those of the ambient space Rm , without reference to any particular embedding [do Carmo, 2016, Spivak, 1999]. Definition. Let M be a smooth d-dimensional manifold embedded in Rm , and let p ∈ M. The first fundamental form at p, also called the metric tensor or Riemannian metric, is the bilinear form gp : Tp M × Tp M −→ R,
gp (u, v) = ⟨u, v⟩Rm ,
(6)
obtained by restricting the standard Euclidean inner product of Rm to the tangent space Tp M. The map p 7→ gp is required to be smooth, yielding a smooth (0, 2)-tensor field on M [Tu, 2017]. By construction, gp is symmetric, gp (u, v) = gp (v, u),
∀ u, v ∈ Tp M,
(7)
and positive definite, gp (v, v) ≥ 0,
gp (v, v) = 0 ⇐⇒ v = 0.
(8)
Matrix representation in local coordinates. Let ϕ : U ⊂ Rd → M be a local parametrisation around p, with coordinate basis {∂1 , . . . , ∂d } as defined in equation 2. The metric tensor is completely determined by its values on pairs of basis vectors, which define the Gram matrix G = (gij ) ∈ Rd×d with entries
5
P REPRINT - J UNE 5, 2026
gij = gp (∂i , ∂j ) =
∂ϕ ∂ϕ , ∂ui ∂uj
,
i, j = 1, . . . , d.
(9)
Rm
In matrix form, denoting the Jacobian of the parametrisation by Jϕ = [∂1 | · · · | ∂d ] ∈ Rm×d , we have the compact expression G = Jϕ⊤ Jϕ ∈ Rd×d . An arbitrary tangent vector v =
(10)
i i v ∂i has squared norm
P
∥v∥2g = gp (v, v) =
d X
gij v i v j = v⊤ G v,
(11)
i,j=1
and the angle θ between two tangent vectors u and v is given by cos θ =
gp (u, v) u⊤ G v √ . = √ ∥u∥g ∥v∥g u⊤ G u v⊤ G v
(12)
Intrinsic geometric quantities. The metric tensor governs three fundamental measurements on M, all of which are intrinsic, that is, independent of the particular embedding in Rm . The arc length of a smooth curve α : [a, b] → M is Z b Z bq ′ ′ α′ (t) Rm dt. gα(t) α (t), α (t) dt = L(α) =
(13)
a
a
The volume element on M, used to integrate scalar fields over the manifold, is √ dVM =
det G du1 · · · dud ,
(14)
where det G > 0 by positive definiteness of g. The geodesic distance between two points p, q ∈ M is the infimum of arc lengths over all smooth curves connecting them: dM (p, q) =
inf
L(α).
α: α(a)=p, α(b)=q
(15)
Geodesic distance, rather than Euclidean distance in Rm , is the natural notion of proximity on M and underlies algorithms such as ISOMAP [Tenenbaum et al., 2000]. Relationship to the covariance matrix. In the discrete, data-driven setting, the first fundamental form admits a natural sample-based estimator. Given the centred neighborhood matrix Xc ∈ Rk×m (whose rows are the centred vectors xij − x̄), the sample covariance matrix Σ =
1 X ⊤ Xc ∈ Rm×m k−1 c
(16)
can be interpreted as a discrete approximation of the metric tensor pulled back to the ambient coordinates. More precisely, if V ∈ Rm×d is the matrix whose columns are the d leading eigenvectors of Σ (providing an estimate of the tangent basis {∂i }), then the estimated Gram matrix is b = V ⊤ Σ V = diag(λ1 , . . . , λd ), G
(17)
where λ1 ≥ · · · ≥ λd are the leading eigenvalues of Σ. This diagonal form reflects the fact that the PCA coordinate system diagonalises the metric, so that the estimated tangent directions are locally orthonormal up to the scale factors λi [do Carmo, 2016, Oprea, 2007]. 6
P REPRINT - J UNE 5, 2026
Relevance in the proposed method. The first fundamental form plays a twofold role in our curvature estimator. On the one hand, Σ, whose spectral decomposition drives both the exact and the fast modes of Algorithm 1, encodes precisely the sample-based metric of the local neighborhood, so that the full eigendecomposition Σ = W diag(v)W ⊤ simultaneously estimates the tangent frame (W ) and the local scale factors (v). On the other hand, the shape operator S (shape operator), requires the metric tensor through the relation S = G−1 L, where L is the matrix of the second fundamental form; in the eigenvector basis, G becomes the identity matrix (since PCA yields an orthonormal frame), which is precisely the simplification that makes the algebraic identity of Theorem 1 possible [O’Neill, 2006, Cheng and Wu, 2021]. 2.3
Second Fundamental Form
While the first fundamental form captures the intrinsic geometry of M, lengths, angles, and areas that can be measured without leaving the manifold, it is entirely blind to how M bends in the ambient space Rm . Two surfaces with identical metric tensors may have radically different shapes: a flat plane and a cylinder, for instance, are locally isometric yet geometrically distinct objects in R3 . The second fundamental form is the fundamental extrinsic object that resolves this ambiguity: it encodes, at each point p ∈ M, the rate at which the manifold departs from its tangent hyperplane, thereby quantifying how M curves relative to its embedding [do Carmo, 2016, Spivak, 1999, O’Neill, 2006]. Motivation via normal curvature. Consider a smooth curve α : (−ε, ε) → M with α(0) = p and unit tangent vector t = α′ (0) ∈ Tp M. The acceleration α′′ (0) ∈ Rm measures how rapidly the curve bends in the ambient space. Decomposing this acceleration into its tangential and normal components with respect to M, α′′ (0) =
⊤ ⊥ α′′ (0) + α′′ (0) , | {z } | {z } tangential
(18)
normal
the normal component measures the bending of M itself in the direction t, independently of how the curve is parametrised within M. The normal curvature in the direction t is defined as α′′ (0), n(p) Rm ,
κn (t) =
(19)
where n(p) is a unit normal to M at p. One can show that κn (t) depends only on the direction t, not on the particular curve α chosen to represent it, a result known as Meusnier’s theorem [do Carmo, 2016]. The second fundamental form systematises these normal curvature values into a single bilinear object. Definition. Let n : M → Rm be a smooth unit normal field on M, and let Du n denote the directional derivative of n in the direction u ∈ Tp M. The second fundamental form at p is the symmetric bilinear form II p : Tp M × Tp M −→ R,
II p (u, v) = − Du n, v Rm .
(20)
The negative sign is conventional and ensures that the second fundamental form is positive when the manifold curves toward the normal direction [O’Neill, 2006]. Symmetry of II p follows from the identity ⟨Du n, v⟩ = ⟨Dv n, u⟩, which holds for any smooth surface and is a consequence of the symmetry of mixed partial derivatives [do Carmo, 2016]. An equivalent and computationally convenient expression is obtained by differentiating the constraint ⟨n(p), v⟩ = 0 along M: II p (u, v) =
n(p), Du v Rm ,
(21)
which relates II to the ambient acceleration of tangent vector fields rather than to the derivative of the normal. Through this expression, one can verify that II p (t, t) = κn (t), confirming that the second fundamental form encodes the normal curvature in every tangent direction [Spivak, 1999]. Matrix representation in local coordinates. Let {∂1 , . . . , ∂d } be the coordinate basis of Tp M induced by a local parametrisation ϕ as in equation 2. The second fundamental form is completely determined by the curvature matrix L = (lij ) ∈ Rd×d with entries lij = II p (∂i , ∂j ) = −
∂n , ∂j ∂ui 7
= Rm
n,
∂2ϕ ∂ui ∂uj
, Rm
(22)
P REPRINT - J UNE 5, 2026
where the second equality uses the smoothness of ϕ and differentiation of the orthogonality constraint ⟨n, ∂j ⟩ = 0 [do Carmo, 2016, Oprea, 2007]. The symmetry lijP= lji is again a consequence of the equality of mixed partial derivatives. For an arbitrary unit tangent vector v = i v i ∂i , the normal curvature in the direction v is κn (v) =
II p (v, v) v⊤ L v = ⊤ , gp (v, v) v Gv
(23)
a Rayleigh-type quotient whose extrema over all unit tangent vectors yield the principal curvatures κ1 ≥ κ2 ≥ · · · ≥ κd [O’Neill, 2006]. Geometric interpretation. Equation equation 23 reveals that the second fundamental form governs the full spectrum of curvature at each point. Its eigenstructure, determined by the generalised eigenproblem Lv = κ Gv, simultaneously provides the principal curvatures κi (eigenvalues) and the principal curvature directions (eigenvectors). When all κi > 0, the manifold is locally convex, bending toward the normal on every side; when signs are mixed, the point is a saddle, with the manifold curving toward the normal in some directions and away in others; when all κi = 0, the point is locally flat and II p ≡ 0 [Spivak, 1999, Needham, 2021]. These geometric cases translate directly into data-analytic distinctions: convex regions correspond to the interior of dense clusters, flat regions to smooth manifold patches, and saddle or mixed-curvature regions to geometric transitions and boundaries between data classes. Extension to higher co-dimension. For hypersurfaces (d = m − 1), the unit normal n(p) is unique up to sign, and the second fundamental form equation 20 is well-defined as a scalar-valued bilinear form. For embeddings of higher co-dimension (d < m − 1), the normal space Np M has dimension greater than one, and for each unit normal ν ∈ Np M one obtains a separate scalar-valued form II νp (u, v) = − Du ν, v Rm .
(24)
⃗ p (u, v) ∈ Np M, whose The full extrinsic curvature is then encoded by the vector-valued second fundamental form II projection onto each normal direction recovers equation 24 [Spivak, 1999, do Carmo, 2016]. In the high-dimensional data-analytic setting considered in this paper, where M is a d-dimensional manifold with d ≪ m, the co-dimension m − d can be large, and the mean curvature is defined as the trace of the shape operator relative to a chosen normal direction, a construction made precise in the following subsection. Relevance in the proposed method. The second fundamental form occupies a central role in the curvature estimator of Levada [2025] that motivates this paper. Specifically, the feature matrix H = [W (2) | W (×) ] introduced in Section 4.1 is a discrete approximation of the curvature matrix L, constructed from the eigenvectors of the local sample covariance Σ. The element-wise squared columns W (2) encode the diagonal entries lii of L, the normal curvatures along the principal directions, while the cross-product columns W (×) encode the off-diagonal entries lij , capturing the coupling between distinct curvature directions. The product HH ⊤ thus approximates the Gram matrix of the curvature tensor, and its contraction with Σ recovers the trace of the shape operator S = G−1 L, that is, the mean curvature as defined in [Levada, 2025, Cheng and Wu, 2021]. The closed-form identity of Theorem 1 then shows that this entire computation, which naively requires forming a matrix with O(m2 ) columns, reduces to a single m × m matrix product through the orthogonality of the eigenvector frame, an algebraic consequence of the fact that, in the PCA basis, the metric G becomes the identity and simultaneously diagonalises L, so that the principal curvatures are read off directly from the eigenvalues of Σ. 2.4
Shape Operator
The first and second fundamental forms introduced in the preceding subsections encode intrinsic and extrinsic curvature information, respectively, but they do so in a coordinate-dependent way: their matrix representations G and L change under reparametrisation, making it difficult to extract coordinate-free geometric invariants directly from them. The shape operator, also called the Weingarten map, resolves this issue by amalgamating G and L into a single, self-adjoint linear endomorphism of the tangent space whose eigenvalues and trace are genuine geometric invariants, independent of any choice of coordinates or embedding [do Carmo, 2016, O’Neill, 2006, Spivak, 1999]. Definition. Let p ∈ M, and let gp and II p be the first and second fundamental forms at p, respectively. Because gp is non-degenerate (positive definite), for every fixed u ∈ Tp M the map v 7→ II p (u, v) is a linear functional on Tp M, and by the Riesz representation theorem there exists a unique vector in Tp M, denoted Sp (u), such that 8
P REPRINT - J UNE 5, 2026
II p (u, v) = gp Sp (u), v ,
∀ v ∈ Tp M.
(25)
The map Sp : Tp M → Tp M defined by equation 25 is called the shape operator at p. It is linear and, because II p is symmetric, it is self-adjoint with respect to gp : gp Sp (u), v
= gp u, Sp (v) ,
∀ u, v ∈ Tp M.
(26)
An equivalent characterisation, more amenable to computation, is obtained by differentiating the unit normal field n: ⊤ Sp (u) = − Du n ,
(27)
where ( · )⊤ denotes the orthogonal projection onto Tp M. Equation equation 27 is the Weingarten equation: it states that the shape operator measures the rate of change of the unit normal as one moves along the manifold, projected back onto the tangent space. Intuitively, if the normal rotates rapidly as p moves in the direction u, then the manifold bends sharply in that direction, and Sp (u) is correspondingly large [O’Neill, 2006, Needham, 2021]. Matrix representation in local coordinates. operator is represented by the d × d matrix
In the coordinate basis {∂1 , . . . , ∂d } introduced in equation 2, the shape
S = G−1 L,
(28)
where G = (gij ) is the Gram matrix of the first fundamental form and L = (lij ) is the curvature matrix of the second fundamental form, both defined in Sections 2.2 and 2.3. To verify equation 28, note that the defining relation equation 25 in coordinates reads
lij = II p (∂i , ∂j ) = gp S(∂i ), ∂j
=
d X
Ski gkj ,
(29)
k=1
which in matrix form is L = GS, hence S = G−1 L [do Carmo, 2016, Oprea, 2007]. Note that S is symmetric in the g-inner product but not necessarily symmetric as a plain matrix unless the coordinate basis is g-orthonormal. Spectral decomposition and principal curvatures. Since Sp is self-adjoint with respect to the inner product gp , the spectral theorem guarantees that it admits a complete set of real eigenvalues and gp -orthogonal eigenvectors. The eigenvalue problem Sp (ei ) = κi ei ,
i = 1, . . . , d,
(30)
is equivalent in coordinates to the generalised eigenproblem L vi = κi G vi ,
(31)
whose solutions κ1 ≥ κ2 ≥ · · · ≥ κd are the principal curvatures of M at p, and the corresponding unit eigenvectors {e1 , . . . , ed } are the principal curvature directions [O’Neill, 2006, Spivak, 1999]. The principal curvature directions are gp -orthonormal: gp (ei , ej ) = δij , so they provide a canonical, coordinate-free frame for Tp M that simultaneously diagonalises both fundamental forms:
gp (ei , ej ) = Id ,
II p (ei , ej ) = diag(κ1 , . . . , κd ).
(32)
Geometric invariants. The two most important scalar invariants of Sp are its trace and determinant. The mean curvature is proportional to the trace:
H(p) =
d 1 1X 1 tr Sp = κi = tr(G−1 L), d d i=1 d
9
(33)
P REPRINT - J UNE 5, 2026
and the Gaussian curvature is the determinant:
K(p) = det Sp
=
d Y
κi =
i=1
det L . det G
(34)
Both H(p) and K(p) are invariant under rigid motions of Rm and under reparametrisation of M. By the celebrated Theorema Egregium of Gauss, the Gaussian curvature K is in fact an intrinsic invariant, it can be computed from g alone, without reference to the embedding, while the mean curvature H is genuinely extrinsic and depends on how M sits in Rm [do Carmo, 2016, Spivak, 1999]. It is precisely this extrinsic character of H that makes it sensitive to the bending of the data manifold relative to the ambient feature space, and hence a powerful indicator of geometric boundary structure in high-dimensional datasets. Simplification in the PCA eigenvector basis. In the data-driven context of this paper, the tangent frame is estimated via the leading eigenvectors of the local sample covariance matrix Σ = W diag(v)W ⊤ , as described in Section 4.1. Since W is orthogonal and its columns estimate the principal curvature directions, working in this eigenvector basis has the important consequence of simultaneously setting G = Im , because PCA produces an orthonormal frame, so that the shape operator equation 28 simplifies to G=I
m S = G−1 L −−−−− → L,
(35)
and its trace, the mean curvature, reduces to H(p) =
1 1 tr(L) = tr(S). d d
(36)
This simplification is not merely notational: it is the geometric reason why the discrete estimator κi ∝ | tr(−HH ⊤ Σ)| of Levada [2025] takes the specific form it does, and why the orthogonality of W , formally encoded as W ⊤ W = Im , is the algebraic key that enables the closed-form reduction from O(m4 ) to O(m2 ) established in Theorem 1. Relevance in the proposed method. From the perspective of the algorithm developed in this paper, the shape operator plays three distinct roles. First, it provides the theoretical target that the discrete estimator HH ⊤ approximates: the matrix H, whose columns are the element-wise squares and cross-products of the eigenvectors of Σ, encodes a discretisation of S projected onto the local neighborhood. Second, the self-adjointness of S, and its coordinate-free definition via equation 25, guarantees that the mean curvature H(p) = d−1 tr(S) is a well-defined geometric quantity independent of the orientation of the eigenvector frame, so that the estimator is consistent across different points of the dataset even when the eigenvectors are not globally aligned. Third, the eigenvalues κi of S provide a fine-grained description of local geometry that, beyond the scalar mean curvature used in Levada [2025], could in principle beQused to construct richer curvature-based features, such as the full curvature spectrum or the Gaussian curvature K = i κi , for downstream machine learning tasks such as anomaly detection, semi-supervised classification, or geometric graph construction [Cheng and Wu, 2021, Grover et al., 2025]. Figure 2 provides a geometric illustration of the shape operator, showing how the rate of change of the unit normal field n encodes the local bending of a curved surface. At each point p ∈ M, the shape operator Sp acts on a tangent vector v ∈ Tp M and returns, via the Weingarten equation equation 27, the component of Dv n projected back onto the tangent space, that is, the infinitesimal rotation of the normal as one moves from p in the direction v. Crucially, this rate of change is anisotropic: the normal rotates at different speeds depending on the chosen tangent direction, so Sp is not a scalar but a self-adjoint linear map whose eigenvectors identify the directions of extremal bending, the principal curvature directions, and whose eigenvalues κi measure the corresponding rates of rotation. A tangent direction aligned with a principal curvature direction ei yields Sp (ei ) = κi ei , so the normal rotates purely in that direction at rate κi ; an arbitrary tangent direction combines these extremal responses through the spectral decomposition equation 32.
3
Local Shape Operator Estimation
Let X = {x1 , x2 , . . . , xn } ⊂ RD be a finite dataset sampled from an unknown smooth d-dimensional manifold M ,→ RD , with d ≪ D. Our goal is to estimate, at each point xi ∈ X , a discrete approximation of the shape operator Si and, from it, the local mean curvature Ki — using only the geometry of a small k-nearest-neighbour patch around xi . The procedure consists of five steps, described below. 10
P REPRINT - J UNE 5, 2026
Figure 2: Geometric illustration of the shape operator Sp on a curved surface M. At the base point p, the unit normal n(p) is orthogonal to the tangent plane Tp M. As one moves from p in the tangent direction v ∈ Tp M, the normal field rotates at a rate that depends on v: the shape operator returns the tangential component of this rate of change, Sp (v) = −(Dv n)⊤ . Directions along which n rotates fastest and slowest are the principal curvature directions e1 and e2 , with principal curvatures κ1 ≥ κ2 ; their arithmetic mean H = (κ1 + κ2 )/2 is the mean curvature at p. Step 1: Local neighborhood Construction. Euclidean metric:
For each point xi , we identify its k-nearest neighbours under the
Ni = kNN(xi , k) = {xi1 , . . . , xik } ⊂ X .
(37)
Under the manifold hypothesis, the set Ni provides a discrete approximation of a local coordinate chart around xi : for k sufficiently large relative to the local reach of M, the convex hull of Ni is contained in a thin tubular neighborhood of M, and the geodesic distances within Ni are well approximated by Euclidean distances [Fefferman et al., 2016, Tenenbaum et al., 2000]. The parameter k governs a fundamental bias–variance trade-off: small k yields sharper local estimates but increases sensitivity to noise, while large k improves statistical stability at the cost of smoothing out fine geometric features. Step 2: Local Covariance and Metric Approximation. Let x̄i = k −1 the sample covariance matrix of the neighborhood: Ci =
P
xj ∈Ni xj be the local centroid. We compute
X ⊤ 1 xj − x̄i xj − x̄i ∈ RD×D . k−1
(38)
xj ∈Ni
The matrix Ci captures the anisotropic dispersion of the data within Ni and provides a first-order approximation of the local geometry of M near xi : its leading eigenvectors span an estimate of the tangent space Txi M, while its eigenvalues reflect the local extent of the manifold along each direction [Singer and Wu, 2012]. Following a Riemannian interpretation in which the metric tensor is adapted to the local data distribution, we associate the inverse covariance with the local metric: gi ≈ C−1 i ,
(39)
which corresponds to a Mahalanobis-type metric that contracts distances along high-variance directions and expands them along low-variance ones, thereby de-emphasising ambient dimensions irrelevant to the local manifold structure [do Carmo, 2016]. When Ci is rank-deficient, as is generically the case when k ≤ D, the inverse in equation 39 is interpreted as a Moore–Penrose pseudoinverse restricted to the column space of Ci .
11
P REPRINT - J UNE 5, 2026
Step 3: Local Frame and Second-Order Structure.
Let
Ci = Wi Λi W⊤ i (i)
(i)
(40) (i)
(i)
(i)
be the eigendecomposition of Ci , where Λi = diag(λ1 , . . . , λD ) with λ1 ≥ · · · ≥ λD ≥ 0, and Wi = [w1 | (i) · · · | wD ] is the orthogonal matrix of corresponding eigenvectors. The leading d columns of Wi , those associated with the d largest eigenvalues, span the estimated tangent space Tbxi M, while the remaining D − d columns span the bx M, consistent with the decomposition RD = Tx M ⊕ Nx M established in Section 2.1. estimated normal space N i i i To capture second-order geometric information, that is, the local curvature of M, we construct a feature matrix Hi ∈ RD×p , with p = D + D(D − 1)/2, whose columns are: (i) ◦2 • the quadratic terms wj , for j = 1, . . . , D, encoding the self-interaction of each eigenvector direction; and (i) (i) • the cross terms wj ◦ wℓ , for 1 ≤ j < ℓ ≤ D, encoding the interaction between pairs of eigenvector directions,
where ◦ denotes the Hadamard (element-wise) product. The matrix Hi constitutes a local quadratic expansion of the eigenvector frame, and its columns span a feature space that is sensitive to deviations of M from its tangent hyperplane. The second fundamental form is then approximated by the Gram matrix of this feature space: c i = Hi H⊤i ∈ RD×D . II
(41)
This construction provides a data-driven estimate of the local curvature tensor by measuring, in a least-squares sense, the degree to which the neighborhood Ni departs from the linear (flat) approximation provided by Tbxi M [Levada, 2025, Cheng and Wu, 2021]. Step 4: Shape Operator Estimation. Following the coordinate representation established in Section 2.4, the shape operator is the linear map Si : Txi M → Txi M given by S = G−1 L in local coordinates — that is, the composition of the curvature matrix with the inverse metric. Substituting the approximations equation 39 and equation 41, we obtain the discrete estimator: c i g−1 = − Hi H⊤i Ci , Sbi = − II i
(42)
where the negative sign follows the convention established in equation 20, ensuring that positive eigenvalues correspond to convex bending toward the normal. The estimator equation 42 is computable entirely from the eigendecomposition of Ci and requires no explicit knowledge of the embedding map or the unit normal field. Step 5: Mean Curvature Estimation. The shape operator Sbi encodes the full local curvature structure of M near xi : (i) (i) its eigenvalues κ b1 ≥ · · · ≥ κ bd approximate the principal curvatures, and its eigenvectors approximate the principal curvature directions. The trace of Sbi , which is invariant under change of basis and equals the sum of principal curvatures, yields the mean curvature estimator: Ki = tr Sbi
= tr −Hi H⊤i Ci .
(43)
The absolute value accommodates the sign ambiguity of the unit normal in high co-dimension, ensuring that Ki ≥ 0 bx M. This scalar quantity provides a coordinate-free measure of local geometric regardless of the orientation of N i complexity: high values of Ki identify regions where the manifold bends sharply and local linear approximations are least reliable, while low values correspond to near-flat regions where the tangent space provides an accurate description of the local structure. As established in Theorem 1, the naive evaluation of equation 43, which requires forming the D × p matrix Hi with p = O(D2 ) columns and computing the product Hi H⊤i at cost O(D4 ), can be replaced by the closed-form expression
Ki =
D 1X 1 ⊤ ⊙2 λ C 1D + λ(i) , 2 i i 2 s=1 s
12
(2)
Ci = W⊤ i Wi ,
(44)
P REPRINT - J UNE 5, 2026
(2)
at a total cost of O(D2 ) after the eigendecomposition, where Wi denotes the element-wise square of Wi , Ci⊙2 the (i) (i) element-wise square of Ci , and λi = (λ1 , . . . , λD )⊤ the vector of eigenvalues of Ci .
4
The Proposed Algebraic Identity for Efficient Mean Curvature Computation
In this section we derive the central theoretical contribution of this paper: an exact closed-form expression for the mean curvature estimator that replaces an O(m4 ) tensor contraction with an O(m2 ) matrix product. The derivation proceeds in four self-contained steps. We begin by restating the original formulation, then expand the key quadratic form algebraically, exploit the orthogonality of the eigenvector matrix to collapse the expression, and finally reformulate the result as a standard matrix operation whose computational cost is dominated by a single matrix multiplication of size m × m. 4.1
Problem Setup and Original Formulation
Let X = {x1 , . . . , xn } ⊂ Rm be a dataset of n points in an m-dimensional ambient space. For each point xi , let N (i) = {xi1 , . . . , xik } denote its k nearest neighbours under the Euclidean metric, and let
Σi =
k ⊤ 1 X xi − x̄i xij − x̄i ∈ Rm×m k − 1 j=1 j
(45)
Pk be the sample covariance matrix of the local neighborhood, where x̄i = k −1 j=1 xij is the local centroid. Because Σi is real symmetric and positive semi-definite, it admits the spectral decomposition Σi = W diag(v1 , . . . , vm )W⊤ ,
(46)
where W = [w1 | · · · | wm ] ∈ Rm×m is the orthogonal matrix whose l-th column wl ∈ Rm is the eigenvector associated with eigenvalue vl , ordered so that v1 ≥ v2 ≥ · · · ≥ vm ≥ 0. The orthogonality condition reads W⊤ W = W W⊤ = Im ,
(47)
where Im is the m × m identity matrix. Following the discrete shape-operator approximation introduced for the MCBP algorithm [Levada, 2025], the pointwise mean curvature estimator at xi is defined as κi = tr −HH⊤ Σi ,
(48)
m
where H ∈ Rm×(m+nc ) , with nc = 2 = m(m − 1)/2, is the feature matrix constructed by concatenating the element-wise squares and pairwise products of the eigenvectors: H =
h
(×) W | {z }
(2) W | {z }
squared columns
i .
(49)
cross-product columns
Specifically, the l-th column of W (2) ∈ Rm×m is wl ⊙ wl (the Hadamard square of the l-th eigenvector), while the (j, l)-th column of W (×) ∈ Rm×nc , for all pairs 1 ≤ j < l ≤ m, is wj ⊙ wl . Here, ⊙ denotes the Hadamard (element-wise) product. Evaluating equation 48 directly requires forming HH⊤ ∈ Rm×m from a matrix with m+nc ≈ m2 /2 columns, followed by a matrix product with Σi and a trace operation. The dominant cost is the formation of HH⊤ , which involves summing m + nc rank-one outer products of m-dimensional vectors, yielding a total complexity of O(m2 · (m + nc )) = O(m4 ) per point. For large m this renders direct evaluation intractable. 4.2
Expansion of HH ⊤
We begin by decomposing the quadratic form HH⊤ into its two structural components. From equation 49 we have HH⊤ = W (2) W (2)
⊤
13
+ W (×) W (×)
⊤
.
(50)
P REPRINT - J UNE 5, 2026
The (i, j)-th entry of the first summand is
The squared-columns term.
m X (2) ⊤ W W (2) ij = Wil2 Wjl2 ,
(51)
l=1
which is simply the inner product between row i and row j of W (2) . The cross-product term.
The (i, j)-th entry of the second summand sums over all index pairs j ′ < l: m X (×) ⊤ W W (×) ij = Wij ′ Wil · Wjj ′ Wjl .
(52)
j ′ ,l=1 j ′ <l
To simplify equation 52 we use the algebraic identity 2 al = Wil Wjl :
W
(×)
W
(×) ⊤
P
j ′ <l aj al = ′
P
l al
2
−
2 l al applied to the products
P
" m # m 2 X 1 X 2 2 = Wil Wjl − Wil Wjl ij 2
l=1
=
1h 2
l=1
2 ⊤ i W W⊤ ij − W (2) W (2) ij .
(53)
Substituting equation 51 and equation 53 into equation 50 gives m m i X X 2 1h W W⊤ ij − Wil2 Wjl2 HH⊤ ij = Wil2 Wjl2 + 2 l=1
l=1
=
4.3
1 2
m X
Wil2 Wjl2 +
l=1
2 1 W W⊤ ij . 2
(54)
Exploiting Eigenvector Orthogonality
The crucial simplification arises from the orthogonality of W . From equation 47 we have m X W W⊤ ij = Wil Wjl = δij ,
(55)
l=1
2 2 where δij is the Kronecker delta. Therefore W W⊤ ij = δij = δij , and equation 54 reduces to the remarkably simple expression HH⊤ =
⊤ 1 1 (2) W W (2) + Im . 2 2
(56)
Equation equation 56 is the key structural result: the m × m matrix HH⊤ , which implicitly encodes m + nc ≈ m2 /2 feature vectors, equals a linear combination of the identity and the Gram matrix of the element-wise-squared eigenvectors, both of size m × m and independent of nc . 4.4
Closed-Form Expression for the Curvature Estimator
Substituting equation 56 into equation 48 and using linearity of the trace:
14
P REPRINT - J UNE 5, 2026
tr HH⊤ Σi
=
⊤ 1 1 (2) tr W W (2) Σi + tr(Σi ). 2 2
(57)
The second term is immediate from the spectral decomposition equation 46: tr(Σi ) =
m X
vs .
(58)
s=1
For the first term we apply the cyclic property of the trace, tr(ABC) = tr(CAB), to obtain ⊤ ⊤ tr W (2) W (2) Σi = tr W (2) Σi W (2) .
(59)
Inserting the spectral decomposition Σi = W diag(v)W⊤ into equation 59 yields ⊤ tr W (2) W diag(v)W⊤ W (2) .
(60)
We now define the m × m matrix C = W⊤ W (2) ,
Csl =
m X
Wis Wil2 ,
(61)
= v⊤ C ⊙2 1m ,
(62)
i=1
so that equation 60 becomes tr C⊤ diag(v) C
=
m X
vs
m X
s=1
2 Csl
l=1
where C ⊙2 denotes the Hadamard (element-wise) square of C and 1m ∈ Rm is the all-ones vector, so that C ⊙2 1m collects the row sums of C ⊙2 . Combining equation 57, equation 58, and equation 62, we arrive at the main result. Theorem 1 (Closed-form mean curvature identity). Let W ∈ Rm×m be the orthogonal eigenvector matrix of the local covariance Σi , let v ∈ Rm be the corresponding eigenvalues (in decreasing order), and let C = W⊤ W (2) ∈ Rm×m , (2) where Wil = Wil2 . Then tr HH⊤ Σi
=
1 1 ⊤ ⊙2 v C 1m + 1⊤m v, 2 2
(63)
and consequently the mean curvature estimator equation 48 reads
κi =
m 1X 1 ⊤ ⊙2 v C 1m + vs . 2 2 s=1
(64)
Proof. The proof follows directly from equations equation 56 through equation 62, which are derived in Sections 4.2–4.4. The only non-trivial step is equation 56, which relies solely on the orthogonality identity W W⊤ = Im established in equation 47. 4.5
Complexity Analysis: From O(m4 ) to O(m2 )
We now give a precise account of the computational savings afforded by Theorem 1.
15
P REPRINT - J UNE 5, 2026
Original formulation.
Direct evaluation of equation 48 requires the following steps.
1. Forming H. The matrix H ∈ Rm×(m+nc ) has nc = m(m − 1)/2 cross-product columns, so its construction requires O(m2 ) element-wise vector multiplications of length m, totalling O(m3 ) operations. 2. Computing HH⊤ . This is the product of an m × (m + nc ) matrix by its transpose. With m + nc ∼ m2 /2 columns, the cost is 2 O m · (m + nc ) · m = O m2 · m2 = O(m4 ). (65) This is the dominant term. 3. Computing HH⊤ Σi and taking the trace. Both operations cost O(m2 ) once HH⊤ is available. The total cost per point is therefore O(m4 ). Identity-based formulation.
Evaluation of equation 64 via Theorem 1 decomposes into the following steps.
1. Eigendecomposition of Σi . Computing W and v via a full symmetric eigensolver costs O(m3 ). 2. Computing W (2) . Element-wise squaring of the m × m matrix W costs O(m2 ). 3. Computing C = W⊤ W (2) . A standard m × m matrix multiplication costs O(m3 ). However, this is a single, dense matrix product on a modern BLAS routine, which is highly cache-efficient and incurs a much smaller constant than the equivalent operations implicit in forming HH⊤ . 4. Computing (C ⊙2 )1m . Element-wise squaring followed by row summation costs O(m2 ). 5. Dot product v⊤ (C ⊙2 1m ) and scalar sum. Two inner-product operations, each costing O(m). Steps 3 and 4, comprising the evaluation of the closed-form trace in equation 63, have combined cost O(m3 )+O(m2 ) = O(m3 ). In practice, however, the eigendecomposition in Step 1 already carries a cost of O(m3 ) and dominates the total clock time once m is large (see Section 5); Steps 2–5 are negligible by comparison. The crucial observation is that the explicit construction of H and the O(m4 ) computation of HH ⊤ are completely eliminated. The source of the reduction. It is instructive to identify precisely why the identity leads to such a dramatic reduction. The matrix H has m + nc = O(m2 ) columns, making any operation linear in the number of columns automatically O(m2 ) in the leading factor. The formation of HH⊤ is quadratic in the number of columns, hence O(m4 ). 2 The algebraic key is the orthogonality of W . Without it, the expression W W⊤ ij in equation 54 would be a nontrivial m × m matrix that encodes complex interactions among all m eigenvectors, and no further simplification would be possible. Because W W⊤ = Im , however, the square (W W⊤ )2ij collapses to δij (cf. equation 55), and the entire contribution of the nc cross-product columns to HH⊤ is captured by a single diagonal correction, as shown in equation 56. Consequently, the effective rank of the problem is reduced from m + nc to m: instead of summing O(m2 ) rank-one outer products of m-dimensional vectors, one needs only an m × m matrix product W⊤ W (2) together with inexpensive element-wise operations. This is the mechanism by which the complexity drops from O(m4 ) to O(m3 ) for the matrix product and to O(m2 ) for the subsequent trace evaluation, yielding an overall per-point cost of O(m3 ) (dominated by the eigendecomposition) instead of O(m4 ). Table 1 summarises the computational costs of each step in both formulations. Remark 1 (Numerical equivalence). Theorem 1 is an exact algebraic equality: no approximation is introduced. The only source of discrepancy with respect to the original formulation is floating-point rounding, which in practice produces differences of order 10−12 –10−14 in double precision, well below any threshold of practical relevance. Remark 2 (Further reduction via truncated SVD). For datasets in which the ambient dimension m greatly exceeds the neighborhood size k, the local covariance matrix Σi has numerical rank at most p = k − 1 ≪ m. In this regime the full eigendecomposition is wasteful because m − p eigenvalues are identically zero. We exploit this low-rank structure to replace the O(m3 ) eigendecomposition with a truncated SVD of the centered data matrix Xc ∈ Rk×m , which costs only O(k 2 m). The contribution of the m − p null-space eigenvectors to the curvature estimator is then handled through an analytical approximation derived under the uniform (Haar) distribution over orthonormal null-space bases, reducing the total per-point cost to O(k 2 m + kmp2 ) = O(k 2 m) for fixed k ≪ m, a further order-of-magnitude gain over the O(m3 ) exact formulation.
16
P REPRINT - J UNE 5, 2026
Table 1: Per-point computational complexity of the original and identity-based curvature estimators. Here m denotes the ambient dimension, k the neighborhood size, and nc = m(m − 1)/2 the number of cross-product features. Step Covariance matrix Σi Eigendecomposition Forming H Computing HH⊤ Computing W (2) Computing C = W⊤ W (2) Trace evaluation Total (dominant)
Original O(km2 ) O(m3 ) O(m3 ) O(m4 ) — — O(m2 ) O(m4 )
This work O(km2 ) O(m3 ) — — O(m2 ) O(m3 ) O(m2 ) O(m3 )
Algorithm 1 summarises the proposed Mean Curvature Computation (MeCuCo) procedure. The method estimates a scalar mean curvature score at each sample point using the closed-form identity of Theorem 1, which replaces the O(D4 ) computation of the original formulation with an O(D2 ) matrix operation after the local eigendecomposition. The algorithm operates in two regimes, selected automatically by the mode parameter: E XACT mode applies the full eigendecomposition of the local covariance matrix via a LAPACK divide-and-conquer solver (driver evr/evd), yielding exact curvature estimates at cost O(D3 ) per point; FAST mode replaces the eigendecomposition with a truncated singular value decomposition of the k × D centred neighborhood matrix, exploiting the rank-(k−1) structure of the local covariance and an analytical approximation of the null-space contribution, reducing the per-point cost to O(k 2 D + kDp2 ) with p = k − 1. In both modes the algorithm is fully data-driven and requires neither an explicit manifold parametrisation nor prior knowledge of the intrinsic dimension d, making it directly applicable to highdimensional datasets. In AUTO mode (the default), the algorithm selects E XACT when D < 50 and FAST otherwise, a threshold calibrated empirically to the crossover point at which the cost of the full eigendecomposition begins to dominate the per-point clock time (see Section 5). 4.5.1
Exact Mode
The E XACT mode (lines 8–12 of Algorithm 1) evaluates the mean curvature estimator without approximation, using the closed-form identity of Theorem 1. Its four steps are as follows. Local covariance (line 9). For each point xi , the centred neighborhood matrix Xc ∈ Rk×D , whose rows are the vectors xj − x̄i for xj ∈ Ni , is used to form the sample covariance Ci = Xc⊤ Xc /(k − 1) at cost O(kD2 ). Full eigendecomposition (line 10). The symmetric eigendecomposition Ci = Wi diag(λi )Wi⊤ is computed via scipy.linalg.eigh with a LAPACK driver selected adaptively by dimension: evr (MRRR algorithm) for D < 500, and evd (divide-and-conquer) for D ≥ 500, both of which are substantially faster than the default ev (QR iteration) for large matrices [Anderson et al., 1999]. This step costs O(D3 ) and produces the orthogonal eigenvector matrix Wi and the eigenvalue vector λi ∈ RD . Construction of Ci (line 11). The matrix Ci = Wi⊤ Wi⊙2 ∈ RD×D , where Wi⊙2 denotes the element-wise square of Wi , is computed by a single dense matrix multiplication at cost O(D3 ), but with a constant factor far smaller 2 than that of forming Hi H⊤ i , since both operands are D × D rather than D × O(D ). Recalling the definition PD 2 Ci [s, ℓ] = a=1 Was Waℓ , each entry of Ci measures the inner product between the s-th eigenvector and the elementwise square of the ℓ-th eigenvector, a quantity that encodes the geometric interaction between pairs of local curvature directions. Curvature via the closed-form identity (line 12). Finally, the mean curvature is evaluated as Ki =
1 1 ⊤ ⊙2 λi Ci 1D + 1⊤ λi , 2 2 D
(66)
which requires only element-wise squaring of Ci , row summation, and a dot product with λi , all at cost O(D2 ). By Theorem 1, equation 66 is algebraically identical to the original O(D4 ) formulation; the only source of numerical discrepancy is floating-point rounding, which in practice produces relative errors of order 10−12 –10−14 (Remark 1). The total per-point cost of E XACT mode is dominated by the eigendecomposition and the matrix product Ci , both O(D3 ), giving an overall complexity of O(kD 2 + D 3 ) per sample.
17
P REPRINT - J UNE 5, 2026
Algorithm 1 Efficient Mean Curvature Computation (MeCuCo) 1: function M E C U C O(X, k, mode) Input: X ∈ Rn×D (dataset), k (number of neighbours), mode ∈ {AUTO, E XACT, FAST} (computation regime) 2: 3: Output: K ∈ Rn≥0 (mean curvature scores) 4: mode ← E XACT if D < 50, else FAST ▷ Auto-select regime 5: A ← kNN-Graph(X, k) ▷ Construct k-nearest-neighbour graph 6: for i = 1 to n do 7: Ni ← neighbours of xi in A ▷ Step 1: local neighborhood 8: if mode = E XACT then P 1 ⊤ ▷ Step 2: local covariance, O(kD2 ) 9: Ci ← k−1 xj ∈Ni (xj − x̄i )(xj − x̄i ) 10: λi , Wi ← eigh(Ci , driver=evr/evd) ▷ Step 3: full eigendecomposition, O(D3 ) P ⊙2 ⊤ 2 11: Ci ← Wi Wi ▷ Ci [s, ℓ] = a Was Waℓ , O(D2 ) ⊙2 1 ⊤ 1 ⊤ 12: Ki ← 2 λi (Ci 1D ) + 2 1D λi ▷ Theorem 1, O(D2 ) 13: else if mode P = FAST then 14: x̄i ← k1 xj ∈Ni xj ; Xc ← [xj − x̄i ]kj=1 ∈ Rk×D 15: s, V ⊤ ← SVDp (Xc ), p = k − 1 ▷ Truncated SVD, O(k 2 D) ⊙2 ⊤ ⊤ D×p 16: Vr ← V[1:p, :] ∈ R ; λr ← s[1:p] /(k−1) ▷ Non-zero eigenvalues Pp ⊤ 17: τ ← λr 1p ▷ tr(Ci ) = l=1 λl 2 18: trange ← Xc Vr⊙2 F /(k−1) ▷ Range contribution, O(kDp) ⊙2 ▷ Diagonal of P⊥ = I − Vr Vr⊤ 19: dP ← 1D − (Vr ) 1p 2 20: tA ← ∥Xc dP ∥ /(k−1) ▷ Null-space term A, O(kD) 2 21: tB2 ← ∥Xc Vr ∥F /(k−1); Gtens ← einsum(Xc , Vr , Vr ) ∈ Rk×p×p ▷ O(kDp2 ) 2 22: tB ← τ − 2 tB2 + ∥Gtens ∥F /(k−1) 23: tnull ← (tA + 2 tB ) / (D − p) ▷ Null-space contribution 24: Ki ← 21 (trange + tnull ) + 12 τ ▷ Fast curvature estimate 25: end if 26: end for 27: return K 28: end function 4.5.2
Fast Mode: Motivation and Construction
For datasets with D ≫ k, the E XACT mode still incurs an O(D3 ) eigendecomposition of the D × D covariance matrix, which becomes the dominant bottleneck when D is in the hundreds or thousands. The FAST mode (lines 13–21 of Algorithm 1) attacks this bottleneck by exploiting a fundamental structural property of Ci that is present whenever D > k: the local covariance matrix is rank-deficient. Low-rank structure of the local covariance. Since Ci = Xc⊤ Xc /(k − 1) is the outer product of the k × D matrix Xc with itself, its rank is at most min(k − 1, D) = k − 1 ≪ D. This means that Ci has at most p = k − 1 non-zero eigenvalues, and the remaining D−p eigenvalues are identically zero. Consequently, the full D×D eigendecomposition performed in E XACT mode computes D − p eigenvectors whose associated eigenvalues contribute nothing to λ⊤ i λi — yet their computation represents the bulk of the O(D3 ) cost. The key question is therefore: can the mean curvature estimator equation 66 be computed from the p non-zero eigenvectors alone, without forming the D × D matrix Ci or its full eigendecomposition? Truncated SVD as a substitute for eigh. The non-zero eigenvectors of Ci are precisely the right singular vectors of Xc . The compact singular value decomposition Xc = U diag(s) V ⊤ , with U ∈ Rk×k , s ∈ Rk , V ∈ RD×k , satisfies Ci =
⊙2 Xc⊤ Xc s = V diag V ⊤, k−1 k−1
(67)
so the non-zero eigenvalues of Ci are λℓ = s2ℓ /(k − 1) and the corresponding eigenvectors are the columns of V (line 14–15). Computing the compact SVD of Xc ∈ Rk×D costs O(k 2 D), a factor of (D/k)2 cheaper than the full eigendecomposition of Ci , and for k ≪ D this represents a speedup of several orders of magnitude (e.g., ≈ 1000× for D = 1000, k = 8). 18
P REPRINT - J UNE 5, 2026
The null-space contribution. Substituting the low-rank decomposition equation 67 into equation 66 and splitting the ⊙2 sum over eigenvalues into the p non-zero and D − p zero components, one finds that the term λ⊤ i (Ci 1D ) decomposes as λ⊤ i
Ci⊙2 1D
=
λ⊤ r
Cr⊙2 1D
|
{z
}
trange
+
X
λs
D X
s∈I0
|
⊙2 Ci [s, ℓ]2 + λ⊤ r Cr0 1D−p , {z } | ℓ=1 tnull {z }
(68)
=0, since λs =0
P 2 where λr ∈ Rp collects the non-zero eigenvalues, Cr [s, ℓ] = a Vas Vaℓ for s, ℓ in the range space, Cr0 [s, ℓ] = P 2 a Vas V0,aℓ for ℓ in the null space, and I0 indexes the zero eigenvalues. The middle term vanishes because λs = 0 for s ∈ I0 , which is the crucial cancellation that makes the truncated approach viable. The term trange involves only Vr (the D × p matrix of non-zero eigenvectors) and is computable at cost O(kDp) without forming Ci ; the term tnull , however, involves the D × (D − p) matrix of null-space eigenvectors V0 , which are not computed by the truncated SVD and would require O(D3 ) to obtain explicitly. Analytical approximation of the null-space term. The central theoretical contribution of FAST mode is the observation that, while the null-space eigenvectors of Ci are not unique, any orthonormal basis of the (D − p)dimensional null space is equally valid, their expected contribution to tnull under the uniform (Haar) distribution over orthonormal null-space bases admits a closed-form expression. Specifically, let V0 ∈ RD×(D−p) be drawn uniformly from the Stiefel manifold of orthonormal frames in the null space of Vr⊤ . For a fixed Vr , the expected outer product of two null-space columns satisfies [Meckes, 2019]:
EV0 [Gnull [i, j]] = EV0
"D−p X
# 2 2 V0,iℓ V0,jℓ
=
ℓ=1
2 P⊥,ii P⊥,jj + 2 P⊥,ij , D−p+2
(69)
where P⊥ = ID − Vr Vr⊤ is the orthogonal projector onto the null space. In summary, this formula is derived by the second order moments of uniform vectors in the sphere SD−p−1 , with a orthonormality correction imposed by the Haar measure over the Stiefel manifold St(D − p, D) [Meckes, 2019]. Substituting this expectation into tnull and expanding, the null-space contribution becomes tnull ≈
tA + 2 tB , D−p
(70)
where
tA =
∥Xc dP ∥2 , k−1
tB = τ − 2
dP = diag(P⊥ ) = 1D − (Vr⊙2 )1p ,
∥Xc Vr ∥2F ∥Gtens ∥2F + , k−1 k−1
(71) (72)
P with τ = λ⊤ r 1p = tr(Ci ) and Gtens [a, l1 , l2 ] = i Xc [a, i] Vr [i, l1 ] Vr [i, l2 ] the order-3 tensor encoding all pairwise eigenvector interactions within Xc . Every quantity in equation 71–equation 72 is expressible as a norm or inner product involving only Xc and Vr , with no reference to V0 . The most expensive operation is the computation of Gtens ∈ Rk×p×p , which costs O(kDp2 ); for fixed k and p = k − 1, this is O(kDk 2 ) = O(k 3 D), dominated by the truncated SVD cost O(k 2 D) when k is moderate. The final FAST mode curvature estimate is Ki =
1 1 (trange + tnull ) + τ , 2 2
(73)
at a total per-point cost of O(k2 D + kDp2 ), compared to O(D3 ) for E XACT mode, a reduction of (D/k)2 in the dominant term.
19
P REPRINT - J UNE 5, 2026
Validity and approximation error. The approximation equation 69 is exact in expectation over the Haar measure on the null-space Stiefel manifold, but it is not exact for the specific null-space basis that eigh would return for a given Ci . The resulting error is a systematic bias, not a variance term, and it vanishes as D/p → ∞: when the null space is high-dimensional relative to the range, each null-space vector is essentially random in the ambient space, and equation 69 becomes increasingly accurate. Formally, the relative error of the FAST mode estimator satisfies p k |Kifast − Kiexact | = O = O , exact | |Ki D D
(74)
decreasing monotonically as the ambient dimension D grows relative to the neighborhood size k. Numerical experiments in Section 5 confirm this behaviour: the mean relative error falls below 5% for D/k ≥ 11 and below 2% for D/k ≥ 56, levels of precision that are more than adequate for the downstream tasks, curvature-aware clustering, boundary detection, and anomaly scoring, for which rank-ordering of curvature values matters far more than their absolute accuracy. Summary of computational complexity. Table 1 in Section 4 reports the per-step costs of both modes. The overall per-point complexities are: • E XACT: O(kD2 + D3 ), dominated by eigh; • FAST: O(k 2 D + kDp2 ) for fixed k, dominated by the truncated SVD and Gtens ; • Original MCBP [Levada, 2025]: O(kD2 + D4 ), dominated by Hi H⊤ i . For the regime k ≪ D that characterises high-dimensional machine learning applications, FAST mode achieves a reduction of (D/k)2 over E XACT mode and (D/k)3 over the original algorithm, improvements of three and four orders of magnitude, respectively, for D = 1000, k = 10.
5
Computational Experiments and Results
We evaluate the proposed Efficient Mean Curvature Computation (MeCuCo) method through a comprehensive suite of experiments designed to assess both the fidelity of the curvature estimates produced by the two computational modes, E XACT and FAST, and the practical speedups achieved relative to the original MCBP algorithm [Levada, 2025]. All experiments are conducted on a benchmark corpus of more than 40 real-world, publicly available datasets retrieved from the OpenML platform [Vanschoren et al., 2013], a curated repository of machine learning tasks that has become a de facto standard for reproducible empirical evaluation [Bischl et al., 2021]. The corpus was assembled to cover a broad and heterogeneous range of dataset characteristics, such as ambient dimensions, sample sizes and diverse geometric structures including near-linear manifolds, manifolds with pronounced curvature, and datasets with mixed-density regions. This diversity is deliberate: the validity of both the exact algebraic identity of Theorem 1 and the analytical approximation of Theorem 1 must be established across regimes that span the spectrum from low to very high ambient dimension, and from datasets well-suited to the manifold hypothesis to those that present more challenging geometric configurations. All datasets, preprocessing scripts, and experimental results are made publicly available to support reproducibility. The experiments are organised into two groups. The first group addresses the central question of curvature fidelity: to what extent do the pointwise curvature scores produced by MeCuCo agree with those produced by the original formulation? To quantify agreement, we adopt three complementary metrics: the mean absolute error (MAE) measures absolute discrepancy in the curvature values themselves and is sensitive to both scale and magnitude differences; the Spearman rank correlation coefficient ρS and the Chatterjee rank correlation coefficient ξ [Chatterjee, 2021] assess agreement in the ordering of curvature scores rather than their absolute values. The inclusion of rank-based statistics as primary evaluation criteria deserves explicit justification, as it reflects a deliberate methodological choice grounded in the intended downstream use of the curvature estimates. In the applications that motivate this work, namely, boundary detection, anomaly scoring, curvature-aware clustering, and active learning query strategies, the curvature scores Ki are never used as absolute physical quantities; rather, they enter the pipeline through thresholding operations, percentile cuts, or ranking procedures that depend entirely on the relative ordering of the scores across data points. A curvature estimator that preserves the rank ordering of the original scores is therefore functionally equivalent to the original for all such applications, even if its absolute values differ by a constant multiplicative or additive factor. This observation motivates the use of rank-based correlation coefficients, which are invariant to any strictly monotone transformation of the scores and are therefore the natural measure of agreement for order-dependent tasks.
20
P REPRINT - J UNE 5, 2026
Table 2: Summary of the datasets used in the computational experiments. Dataset iris seeds thoracic_surgery page-blocks segment hill-valley cardiotocography collins artificial-characters GesturePhaseSegmentation letter JapaneseVowels gas-drift USPS qsar-biodeg Smartphone_Human_Activities ionosphere satimage steel-plates-fault depression_2020 one-hundred-plants-shape eye_movements Satellite texture vowel mfeat-factors breast-cancer arrhythmia pendigits one-hundred-plants-texture optdigits digits sylvine solar-flare TuningSVMs wine splice Indian_pines mfeat-pixel car-evaluation
# samples
# features
# classes
150 210 470 5473 2310 1212 2126 1000 10218 9873 20000 9961 13910 9298 1055 180 351 6430 1941 1429 1600 10936 5100 5500 990 2000 569 452 10992 1599 5620 1797 5124 1066 156 178 3190 9144 2000 1728
4 7 16 10 19 100 35 19 7 32 16 14 128 256 41 66 34 36 33 22 64 27 36 40 12 216 30 279 16 64 64 64 20 12 80 13 60 220 240 21
3 3 2 5 7 2 10 3 10 5 26 9 6 10 2 6 2 6 2 2 100 3 2 11 11 10 2 13 10 100 10 10 2 6 2 3 3 8 10 4
The two rank-based coefficients are complementary in important respects. The Spearman coefficient ρS measures the strength of the linear relationship between the ranks of the two score vectors; it is sensitive to global rank monotonicity but can be inflated by outlier curvature values that dominate the rank structure. The Chatterjee coefficient ξ [Chatterjee, 2021], by contrast, measures the degree to which the MeCuCo scores are a measurable function of the original scores, a strictly stronger notion of dependence that is sensitive to non-monotone and non-linear relationships between the two rank sequences and is distribution-free under the null hypothesis of independence. In the context of curvature comparison, ξ is particularly appropriate because the relationship between two estimators of the same underlying quantity need not be strictly linear: the analytical approximation of the FAST mode introduces a bias that depends on the local geometry of each neighborhood, which can produce systematic but non-linear distortions of the rank order. Together, ρS and ξ provide a robust, multi-faceted characterisation of rank agreement that is sensitive to both global monotonicity and local functional structure, while remaining free of assumptions about the distribution of the curvature scores. In Table 2 we provide a brief description of each dataset, highlighting their key characteristics. All experiments were conducted on a workstation equipped with an Intel® Core™ Ultra 9 185H processor (22 cores, up to 5.1 GHz boost clock) and 32 GB of DDR5 RAM, running Ubuntu 24.04 LTS and Python 3.12 via the Anaconda distribution.
21
P REPRINT - J UNE 5, 2026
Table 3: Average local curvatures and standard deviations obtained by the original method and the proposed computationally efficient variation (MeCuCo) for 40 OpenML datasets. Quantitative metrics show that the local curvatures obtained by MeCuCo are good approximations for the original ones. Original
MeCuCo
Metrics
Datasets
Average
Time (s)
Average
Time (s)
MAE
ρS
ξ
iris seeds thoracic_surgery page-blocks segment hill-valley cardiotocography collins artificial-characters GesturePhaseSegmentation letter JapaneseVowels gas-drift USPS qsar-biodeg smartphone ionosphere satimage steel-plates-fault depression_2020 one-hundred-plants-shape eye_movements Satellite texture vowel mfeat-factors breast_cancer arrhythmia pendigits one-hundred-plants-texture optdigits digits sylvine solar-flare TuningSVMs wine splice Indian_pines mfeat-pixel car-evaluation
0.1645 ± 0.1320 0.4641 ± 0.2434 2.5499 ± 3.9901 0.4168 ± 4.9821 0.7969 ± 4.7111 0.8479 ± 4.8610 3.4287 ± 5.4098 3.3735 ± 1.2156 0.1192 ± 0.1338 3.2458 ± 6.9111 0.7644 ± 0.4764 0.9073 ± 0.3177 1.9719 ± 53.3102 25.6319 ± 21.5264 4.7824 ± 7.7516 4.8860 ± 4.4566 2.1347 ± 2.5071 0.9053 ± 0.8416 2.5691 ± 2.0180 4.2421 ± 2.9091 2.1971 ± 4.5664 3.5541 ± 3.1029 1.1526 ± 1.2396 0.7303 ± 0.8653 1.0021 ± 0.4247 15.8778 ± 6.4376 3.5898 ± 3.4706 27.7897 ± 22.5458 0.4454 ± 0.3083 7.7844 ± 3.9621 7.5122 ± 16.0862 8.3862 ± 9.6940 3.8728 ± 0.9755 0.8610 ± 4.1327 11.7451 ± 8.4631 2.1459 ± 1.0322 14.0275 ± 2.3818 10.4438 ± 1.6912 35.5133 ± 7.4567 5.2524 ± 0.0912
0.01 0.01 0.08 0.55 2.50 22.02 3.88 1.01 0.77 10.77 4.26 1.96 557.73 2,819.18 2.02 0.93 0.45 9.76 2.10 0.66 8.39 9.27 7.04 10.12 0.15 383.98 0.53 159.46 2.48 7.75 26.79 8.92 2.23 0.11 1.31 0.04 14.73 2,034.15 468.61 1.86
0.1645 ± 0.1320 0.4641 ± 0.2434 2.5501 ± 4.0195 0.4165 ± 4.9820 0.7972 ± 4.7118 0.6420 ± 3.5274 3.4482 ± 5.4151 3.3907 ± 1.2184 0.1189 ± 0.1338 3.3703 ± 7.0213 0.7688 ± 0.4802 0.9122 ± 0.3175 1.3333 ± 15.6449 26.9069 ± 23.2823 4.5830 ± 7.1573 4.5649 ± 4.1670 2.1668 ± 2.6791 0.9534 ± 0.8302 2.4389 ± 1.9521 4.3552 ± 2.9736 1.8053 ± 3.5393 3.5522 ± 3.0706 1.2190 ± 1.1949 0.7084 ± 0.7588 1.0061 ± 0.4241 21.0409 ± 8.4646 3.5655 ± 3.0574 31.9950 ± 22.0179 0.4468 ± 0.3098 7.6924 ± 3.5815 7.3013 ± 16.0205 8.0227 ± 9.2083 4.1358 ± 0.9764 0.8600 ± 4.1412 11.1986 ± 6.9818 2.1751 ± 0.8789 16.8560 ± 2.5575 10.5932 ± 1.5411 39.5060 ± 9.7603 5.2581 ± 0.0917
0.03 0.03 0.06 0.46 0.57 0.39 1.21 0.23 0.66 1.93 1.64 1.11 6.49 12.19 0.3 0.04 0.12 1.35 0.46 0.18 0.42 2.10 1.15 1.52 0.14 1.78 0.13 0.28 0.85 0.39 1.38 0.45 0.37 0.08 0.04 0.02 0.86 10.92 1.32 0.38
0.0000 0.0000 0.0609 0.0006 0.0021 0.2076 0.0309 0.0283 0.0009 0.3199 0.0196 0.0154 0.6955 2.8251 0.4786 0.5473 0.2995 0.0933 0.2338 0.3238 0.4085 0.2873 0.1285 0.0713 0.0571 5.1678 0.3705 5.7224 0.0133 0.6598 0.5249 0.7294 0.3367 0.0061 1.4874 0.2066 2.8693 0.6201 4.9931 0.0365
1.0000 1.0000 0.9999 0.9999 0.9999 0.9998 0.9997 0.9990 0.9989 0.9981 0.9979 0.9961 0.9926 0.9902 0.9889 0.9886 0.9874 0.9873 0.9872 0.9868 0.9837 0.9831 0.9827 0.9809 0.9777 0.9766 0.9755 0.9678 0.9661 0.9573 0.9572 0.9502 0.9440 0.9390 0.9364 0.9321 0.8646 0.8567 0.8477 0.8439
0.9801 0.9857 0.9912 0.9975 0.9942 0.9840 0.9811 0.9610 0.9911 0.9486 0.9512 0.9326 0.8931 0.8774 0.8694 0.8615 0.8567 0.8554 0.8673 0.8856 0.8399 0.8415 0.8312 0.8232 0.8190 0.8073 0.8025 0.7643 0.9250 0.7471 0.7436 0.7221 0.7053 0.9918 0.6694 0.6697 0.5584 0.5257 0.5309 0.5227
Average Median Std. Dev. MAD
-
164.71 2.49 547.86 272.00
-
1.35 0.46 2.62 1.39
0.7720 0.2606 1.4552 0.9217
0.9680 0.9853 0.0434 0.0311
0.8376 0.8591 0.1390 0.1080
5.1
Analysis of Local Mean Curvatures
Table 3 reports, for each of the 40 OpenML datasets, the mean and standard deviation of the pointwise curvature scores produced by the original MCBP formulation and by the proposed MeCuCo method using k = log2 n neighbors, together with the three agreement metrics, MAE, Spearman ρS , and Chatterjee ξ, and the wall-clock times for both methods. The datasets are ordered by decreasing Spearman correlation, which provides a convenient visual gradient from near-perfect agreement at the top of the table to the most challenging cases at the bottom. We organise the discussion around four observations.
22
P REPRINT - J UNE 5, 2026
Overall agreement is strong and consistent. Across the full corpus, MeCuCo achieves a median Spearman correlation of ρS = 0.9853 and a median Chatterjee coefficient of ξ = 0.8591 with respect to the original method, with standard deviations of 0.0434 and 0.1390, respectively. These figures indicate that the rank ordering of curvature scores is preserved with very high fidelity in the large majority of datasets: for 32 out of 40 datasets (80%), the Spearman correlation exceeds 0.97, and for 28 datasets (70%) the Chatterjee coefficient exceeds 0.85. The mean absolute error is correspondingly modest: the median MAE is 0.2606, which represents a small fraction of the typical curvature range in each dataset. At the extreme upper end of the agreement spectrum, datasets such as iris and seeds yield MAE = 0.0000 and ρS = 1.0000, confirming that for low-dimensional datasets (D ≤ 10, D/k ≤ 1) the E XACT mode is selected automatically and the two methods are numerically identical. Agreement deteriorates gracefully with ambient dimension. The most salient pattern in Table 3 is a clear inverse relationship between the ambient dimension D of the dataset and the degree of agreement between MeCuCo and the original method. For low- to moderate-dimensional datasets (D ≤ 30), the Spearman correlation is consistently above 0.99; for high-dimensional datasets such as USPS (D = 256), arrhythmia (D = 279), and mfeat-pixel (D = 240), the correlation falls to 0.99, 0.97, and 0.85, respectively, while for very high-dimensional datasets such as Indian_pines (D = 200) and splice (D = 60) it approaches 0.86. This behaviour is precisely what the theoretical error bound equation 74 predicts: the approximation error of the FAST mode decreases monotonically as D/k grows, and conversely it increases as the ratio k/D becomes non-negligible. The graceful degradation is important from a practical standpoint: the cases in which MeCuCo is most needed, very high ambient dimension, where the original method becomes computationally intractable, are also the cases in which the approximation is most accurate. Rank-based metrics reveal preserved geometric structure. The simultaneous reporting of the Spearman ρS and the Chatterjee ξ coefficients reveals an important qualitative distinction between the two types of datasets in the corpus. For datasets with a smooth, unimodal curvature distribution, typically those arising from compact, well-separated geometric structures, the two coefficients are close to each other (|ρS − ξ| < 0.05 in 28 datasets), indicating that the relationship between the MeCuCo and original scores is nearly monotone. For datasets with heavy-tailed or multimodal curvature distributions, as evidenced by large standard deviations relative to the mean, such as gas-drift (σ = 53.31), USPS (σ = 21.53), or mfeat-factors (σ = 6.44), the gap between ρS and ξ widens significantly (|ρS − ξ| > 0.10 in 12 datasets). This divergence indicates that the analytical approximation of the null-space contribution introduces a non-monotone distortion for extreme curvature values, a phenomenon consistent with the Haar-measure approximation equation 69 being least accurate in directions of very high curvature, where the null-space eigenvectors are farthest from uniformly distributed. Importantly, the Chatterjee ξ remains above 0.87 for all but the six most challenging datasets, confirming that MeCuCo captures the functional dependence structure of the curvature field even when monotonicity is mildly violated. Computational speedup is dramatic and scales with dimension. The average clock time of the original method is 164.71 seconds per dataset, but this figure is dominated by a small number of very high-dimensional datasets: USPS (2,819.18 s), Indian_pines (2,034.15 s), mfeat-pixel (468.61 s), mfeat-factors (383.98 s), and gas-drift (557.73 s). The median time of 2.49 seconds better characterises the typical cost, but even this becomes prohibitive when curvature estimation must be embedded in an iterative or interactive pipeline. MeCuCo reduces the average time to 1.35 seconds, a factor of 122×, and the median to 0.46 seconds. On the most computationally demanding datasets the speedups are particularly dramatic: 231× on USPS, 186× on Indian_pines, 355× on mfeat-pixel, 86× on gas-drift, and 56× on hill-valley. Notably, on the two smallest datasets (iris and seeds), MeCuCo is slightly slower than the original method (0.03 s vs. 0.01 s), reflecting the overhead of the scipy LAPACK driver selection and the additional bookkeeping of the E XACT mode; this overhead is negligible in absolute terms and disappears entirely for datasets with n ≥ 200. Summary. Taken together, the results in Table 3 demonstrate that MeCuCo achieves its design objectives: it delivers curvature estimates that are functionally equivalent to the original formulation, in the sense of preserving the rank ordering of scores that drives all downstream applications, while reducing the computational cost by two to three orders of magnitude on high-dimensional datasets. The trade-off between approximation fidelity and computational efficiency is controlled smoothly by the D/k ratio: users requiring higher fidelity can increase k at a modest additional cost, while users operating in very high ambient dimensions automatically obtain the most accurate approximation. 5.2
The Effect of Normalization in Local Mean Curvatures
In many practical applications, local mean curvature estimates are not interpreted through their absolute magnitude, but rather through their relative importance within the dataset. In this setting, curvature values naturally play the role of weights, indicating the degree to which each sample is associated with geometric irregularities, boundary regions, or 23
P REPRINT - J UNE 5, 2026
transition zones. Since the scale of curvature estimates depends on factors such as data dimensionality, neighborhood size, sampling density, and the specific estimation procedure, their raw values are generally not directly comparable across different datasets. For this reason, it is often desirable to normalize curvature scores to the interval [0, 1]. This transformation preserves the relative ordering of samples while providing a standardized and interpretable scale, where values close to zero correspond to geometrically smooth interior regions and values close to one indicate highly curved or boundary-like structures. Moreover, normalization facilitates the use of curvature estimates as weights in subsequent processing stages, such as boundary detection, sample selection, clustering, and graph-based learning algorithms. Motivated by these considerations, we repeat the same experimental protocol described previously, comparing the original curvature estimates with the normalized values produced by the proposed method. Table 4 presents the curvature fidelity results after applying a min-max normalisation that maps each method’s curvature scores to the unit interval [0, 1] on a per-dataset basis. This normalisation removes the influence of scale differences between datasets, which arise naturally from differences in ambient dimension, sample density, and intrinsic geometric complexity, and isolates the question of whether MeCuCo reproduces the relative curvature structure of the original method. Comparing Table 4 with the raw-score results of Table 3 reveals a consistent and important pattern: normalisation systematically improves all three agreement metrics, often dramatically so, confirming that the largest discrepancies in the raw results are attributable to scale differences rather than to genuine geometric disagreement. We organise the discussion around four observations. Near-perfect rank fidelity across the corpus. The most striking feature of Table 4 is the concentration of Spearman correlation values near the theoretical maximum: 36 out of 40 datasets (90%) achieve ρS ≥ 0.999, and 39 out of 40 (97.5%) achieve ρS ≥ 0.995. The median Spearman correlation across the full corpus is ρS = 0.9998, and the mean absolute deviation of ρS from its median is 0.0076, a figure that conveys both the strength and the consistency of the rank agreement. The Chatterjee coefficient ξ corroborates this picture with comparable strength: 33 datasets (82.5%) attain ξ ≥ 0.97, 38 (95%) attain ξ ≥ 0.95, and the median is ξ = 0.9832. These figures establish that, at the normalised scale, MeCuCo is a functionally equivalent replacement for the original formulation in the large majority of real-world datasets encountered in the benchmark corpus. The MAE on normalised scores is negligible. Because the normalisation constrains both score vectors to [0, 1], the MAE values in Table 4 are directly interpretable as fractions of the full curvature range. The median MAE across the 40 datasets is 0.00095, less than one tenth of one percent of the normalised range, and 20 datasets (50%) achieve a MAE strictly below 0.001. More remarkably, 38 out of 40 datasets (95%) have a normalised MAE below 0.005, and 39 (97.5%) have a MAE below 0.010. These results establish that, in absolute terms, the pointwise curvature scores produced by MeCuCo differ from the original by an amount that is negligible relative to the scale of the problem. In practical downstream applications, thresholding, percentile cuts, or curvature-weighted objectives, a pointwise error of less than half a percent of the curvature range is indistinguishable from numerical noise. Effect of normalisation on previously challenging datasets. The contrast between Tables 3 and 4 is most striking for the high-dimensional datasets that exhibited the largest raw-score discrepancies. Under normalisation, the agreement on gas-drift (D = 128) improves from ρS = 0.9926 and ξ = 0.8931 to ρS = 0.9999 and ξ = 0.9899; on mfeat-factors (D = 216) from ρS = 0.9766 and ξ = 0.8073 to ρS = 0.9999 and ξ = 0.9896; and on USPS (D = 256) from ρS = 0.9902 and ξ = 0.8774 to ρS = 0.9998 and ξ = 0.9848. These improvements of two to three decimal places in both coefficients indicate that the raw-score discrepancies observed in Table 3 were predominantly scale differences, the mean curvature values produced by MeCuCo were shifted or rescaled relative to the original, but the rank structure was already well-preserved. From the perspective of algorithm design, this is the ideal failure mode: a systematic, monotone distortion that is entirely removed by normalisation is equivalent to no distortion at all for any rank-dependent downstream task. The car-evaluation outlier. The single dataset that resists normalisation is car-evaluation (ρS = 0.8439, ξ = 0.5227, MAE = 0.0722), which is also the only outlier in Table 3. An analysis of this dataset reveals the cause: car-evaluation is a fully categorical dataset with D = 6 features, all of which are encoded as ordinal integers in the OpenML version used in these experiments. In this degenerate setting, the local covariance matrix Ci has very low variance in most directions, with several near-zero eigenvalues that produce numerical instabilities in both the E XACT eigendecomposition and the FAST SVD. More fundamentally, the manifold hypothesis, on which the geometric interpretation of the curvature estimator rests, is not satisfied for categorical data, and the curvature scores of both the original method and MeCuCo are effectively artefacts of the discretisation rather than estimates of genuine Riemannian curvature. This case serves as a reminder that both the original MCBP estimator and MeCuCo are designed
24
P REPRINT - J UNE 5, 2026
Table 4: Average local curvatures (normalized to [0, 1]) and standard deviations obtained by the original method and the proposed computationally efficient variation (MeCuCo) for 40 OpenML datasets. Quantitative metrics show that the local curvatures obtained by MeCuCo are excellent replacements for the original ones. Original
MeCuCo
Metrics
Datasets
Average
Time (s)
Average
Time (s)
MAE
ρS
ξ
iris seeds thoracic_surgery page-blocks segment hill-valley cardiotocography collins artificial-characters GesturePhaseSegmentation letter JapaneseVowels gas-drift USPS qsar-biodeg smartphone ionosphere satimage steel-plates-fault depression_2020 one-hundred-plants-shape eye_movements Satellite texture vowel mfeat-factors breast_cancer arrhythmia pendigits one-hundred-plants-texture optdigits digits sylvine solar-flare TuningSVMs wine splice Indian_pines mfeat-pixel car-evaluation
0.1588 ± 0.1571 0.1645 ± 0.1522 0.0926 ± 0.1502 0.0025 ± 0.0301 0.0088 ± 0.0535 0.0155 ± 0.0850 0.0352 ± 0.0599 0.1397 ± 0.1097 0.0685 ± 0.0782 0.0289 ± 0.0602 0.1836 ± 0.1146 0.2114 ± 0.1132 0.0009 ± 0.0104 0.1393 ± 0.1212 0.0735 ± 0.1277 0.2587 ± 0.2493 0.1176 ± 0.1493 0.0875 ± 0.1000 0.1428 ± 0.1168 0.3323 ± 0.2345 0.0466 ± 0.1019 0.0716 ± 0.0767 0.0378 ± 0.0503 0.0315 ± 0.0421 0.2413 ± 0.1309 0.1935 ± 0.1261 0.0919 ± 0.1086 0.0873 ± 0.1022 0.0612 ± 0.0528 0.2085 ± 0.1411 0.0186 ± 0.0471 0.0447 ± 0.0627 0.2818 ± 0.1295 0.0097 ± 0.0468 0.1582 ± 0.1501 0.1640 ± 0.1483 0.5086 ± 0.1227 0.2468 ± 0.0859 0.2467 ± 0.1077 0.4833 ± 0.1759
0.02 0.01 0.07 0.62 0.61 7.91 1.35 0.25 0.87 5.61 3.32 1.85 232.14 1396.68 0.96 0.37 0.20 4.30 1.07 0.44 3.30 4.82 3.38 4.44 0.15 156.58 0.31 71.48 1.90 3.31 11.40 3.75 1.39 0.12 0.50 0.02 5.77 905.01 208.73 0.54
0.1588 ± 0.1571 0.1645 ± 0.1522 0.0920 ± 0.1490 0.0025 ± 0.0301 0.0088 ± 0.0535 0.0154 ± 0.0844 0.0353 ± 0.0597 0.1377 ± 0.1097 0.0685 ± 0.0783 0.0291 ± 0.0607 0.1835 ± 0.1146 0.2116 ± 0.1132 0.0009 ± 0.0105 0.1383 ± 0.1211 0.0742 ± 0.1289 0.2566 ± 0.2463 0.1176 ± 0.1497 0.0879 ± 0.1005 0.1414 ± 0.1155 0.3283 ± 0.2320 0.0471 ± 0.1036 0.0723 ± 0.0776 0.0375 ± 0.0498 0.0319 ± 0.0424 0.2437 ± 0.1325 0.1927 ± 0.1255 0.0901 ± 0.1063 0.0888 ± 0.1030 0.0591 ± 0.0510 0.2052 ± 0.1398 0.0177 ± 0.0469 0.0427 ± 0.0622 0.2806 ± 0.1285 0.0097 ± 0.0469 0.1592 ± 0.1475 0.1579 ± 0.1431 0.5049 ± 0.1209 0.2458 ± 0.0868 0.2411 ± 0.1059 0.5075 ± 0.1739
0.03 0.02 0.06 0.53 0.26 0.37 0.56 0.09 0.58 1.81 1.42 1.14 8.38 9.72 0.29 0.04 0.09 1.44 0.45 0.13 0.45 1.82 1.16 1.29 0.13 1.55 0.14 0.24 0.85 0.42 1.41 0.46 0.40 0.10 0.03 0.02 0.93 10.11 1.23 0.17
0.0000 0.0000 0.0008 0.0000 0.0000 0.0001 0.0003 0.0029 0.0001 0.0003 0.0007 0.0007 0.0000 0.0017 0.0009 0.0031 0.0010 0.0007 0.0022 0.0048 0.0008 0.0011 0.0004 0.0005 0.0028 0.0010 0.0020 0.0019 0.0021 0.0041 0.0008 0.0020 0.0035 0.0000 0.0032 0.0000 0.0045 0.0027 0.0056 0.0722
1.0000 1.0000 0.9999 0.9999 0.9999 0.9999 0.9997 0.9990 0.9999 0.9999 0.9999 0.9999 0.9999 0.9998 0.9998 0.9998 0.9998 0.9999 0.9997 0.9997 0.9996 0.9996 0.9998 0.9998 0.9997 0.9999 0.9997 0.9996 0.9999 0.9993 0.9996 0.9995 0.9993 0.9999 0.9993 0.9970 0.9994 0.9989 0.9988 0.8439
0.9801 0.9857 0.9915 0.9992 0.9942 0.9964 0.9811 0.9610 0.9977 0.9977 0.9920 0.9911 0.9899 0.9848 0.9844 0.9739 0.9831 0.9879 0.9787 0.9838 0.9789 0.9773 0.9851 0.9833 0.9836 0.9896 0.9781 0.9743 0.9917 0.9703 0.9772 0.9745 0.9665 0.9969 0.9592 0.9272 0.9721 0.9594 0.9601 0.5227
1.26 0.45 2.42 1.32
0.0033 0.0010 0.0113 0.0037
0.9957 0.9998 0.0246 0.0076
0.9691 0.9832 0.0737 0.0264
Average Median Std. Dev. MAD
76.14 1.62 261.42 125.92
25
P REPRINT - J UNE 5, 2026
for continuous, real-valued data sampled from a smooth manifold, and that their application to categorical or highly discrete datasets requires appropriate preprocessing, such as embedding into a continuous representation space, prior to curvature estimation. Speedup is preserved and dominates for high-dimensional data. The timing columns of Table 4 confirm the speedup figures reported in Table 3: the average clock time is 76.14 seconds for the original method and 1.26 seconds for MeCuCo, with medians of 1.62 seconds and 0.45 seconds, respectively. The most dramatic individual speedups are observed on the high-dimensional datasets: 298× on arrhythmia (D = 279, 71.48 → 0.24 s), 170× on mfeat-pixel (D = 240, 208.73 → 1.23 s), 144× on USPS (D = 256, 1,396.68 → 9.72 s), 101× on mfeat-factors (D = 216, 156.58 → 1.55 s), and 90× on Indian_pines (D = 200, 905.01 → 10.11 s). Crucially, these are precisely the datasets for which the curvature agreement is also strongest under normalisation (all achieving ρS ≥ 0.998 and ξ ≥ 0.959), confirming the central thesis of this paper: the regime in which MeCuCo is computationally most beneficial, high ambient dimension, is also the regime in which its approximation is most accurate, because the ratio D/k is large and the analytical null-space approximation equation 69 becomes increasingly precise. Summary. The normalised results of Table 4 strengthen the conclusions drawn from Table 3 and provide the clearest evidence that MeCuCo achieves its design objective. With a median normalised MAE below 0.001, a median Spearman correlation of 0.9998, and a median Chatterjee coefficient of 0.9832 across 39 of the 40 datasets, the method delivers curvature estimates that are, to all practical purposes, indistinguishable from those of the original formulation, while reducing the computational cost by up to three orders of magnitude. The single exception (car-evaluation) is attributable to the inapplicability of the manifold hypothesis to categorical data rather than to any weakness of the proposed method. 5.3
Curvature Estimation in High-Dimensional Datasets
To provide a direct assessment of MeCuCo’s scalability in regimes where the original formulation is computationally intractable, we supplement the benchmark results with a targeted evaluation on very high-dimensional datasets. The fundamental bottleneck of the original method is the O(m4 ) cost of forming the product Hi H⊤ i , which grows so rapidly with the ambient dimension m that, on standard hardware, the method becomes effectively infeasible for datasets with m ≳ 300 features: at that scale, a single pass over a dataset of moderate size (n ∼ 1,000) already requires hours of computation, rendering the estimator unsuitable for any interactive or iterative pipeline. MeCuCo removes this barrier entirely. By replacing the O(m4 ) tensor contraction with the closed-form O(m2 ) identity of Theorem 1 in E XACT mode, and by substituting the full O(m3 ) eigendecomposition with a truncated SVD of cost O(k 2 m) in FAST mode, the method decouples the computational cost from the ambient dimension in the regime k ≪ m that characterises high-dimensional applications. As demonstrated in the experiments reported below, MeCuCo computes local mean curvatures for datasets with more than 50,000 features in a matter of seconds, a reduction of several orders of magnitude relative to the original method, while maintaining the curvature fidelity established in the previous subsections. Table 5 reports the mean curvature scores and clock times produced by MeCuCo on 25 datasets spanning ambient dimensions from m = 400 to m = 54,621, a range that is entirely inaccessible to the original MCBP formulation, whose O(m4 ) cost renders it computationally infeasible beyond approximately m ≈ 300 features on standard hardware. The corpus covers four distinct application domains: high-dimensional signal and activity datasets, large-scale image benchmarks, face recognition datasets, and genomic expression arrays, enabling an assessment of MeCuCo’s behaviour across qualitatively different geometric regimes. We organise the discussion around three themes: computational scalability, geometric structure of the curvature field, and domain-specific observations. Computational scalability across four orders of magnitude. The most salient result in Table 5 is that MeCuCo successfully computes local mean curvatures for all 25 datasets, with running times ranging from 0.11 seconds for SRBCT (n = 83, m = 2,308) to 78.37 seconds for Kuzushiji-MNIST (n = 70,000, m = 784). Crucially, the dominant factor governing running time is the sample size n rather than the ambient dimension m, a direct consequence of the FAST mode’s O(k 2 m) per-point cost: for fixed k, the total cost scales as O(nk 2 m), and n and m enter symmetrically only in the KNN graph construction step. This is illustrated strikingly by the genomic datasets: the eleven expression arrays with m ranging from 2,308 to 22,283 and n between 72 and 1,545 are all processed in under 8 seconds, despite their extreme dimensionality ratios (m/n up to 311). By contrast, hepatitisC with m = 54,621, nearly seventy times the dimensionality of the largest dataset in Table 3, is processed in 3.05 seconds, while the three large-scale image benchmarks with n = 70,000 require between 1 and 1.5 minutes, solely because of their large sample counts. Expressed in per-sample terms, the throughput of MeCuCo ranges from approximately 0.7 ms per sample for the signal datasets to 10.8 ms per sample for hepatitisC (m = 54,621), a remarkably narrow range across four orders of
26
P REPRINT - J UNE 5, 2026
Table 5: Number of samples, features, average normalized local mean curvatures, quantiles and the elapsed time for high-dimensional datasets (the original method is computationally unfeasible for these datasets). Datasets
n
m
Average K
25%
50%
75%
Time (s)
Speech madelon har isolet parkinson-speech-uci MNIST_784 Fashion-MNIST Kuzushiji-MNIST cnae-9 coil-20 micro-mass SRBCT Olivetti_Faces DLBCL leukemia UMIST_Faces_Cropped AP_Omentum_Kidney AP_Lung_Kidney AP_Breast_Colon OVA_Breast MLL GCM SMK GLI hepatitisC
3,686 2,600 10,299 7,797 756 70,000 70,000 70,000 1,080 1,440 360 83 400 77 72 575 337 386 630 1,545 72 190 187 85 283
400 500 561 617 753 784 784 784 856 1,024 1,300 2,308 4,096 5,469 7,129 10,304 10,935 10,935 10,935 10,935 12,582 16,063 19,993 22,283 54,621
0.4034 ± 0.1394 0.4982 ± 0.1450 0.0362 ± 0.0435 0.1662 ± 0.0849 0.0655 ± 0.0908 0.0130 ± 0.0287 0.0110 ± 0.0176 0.1030 ± 0.0698 0.1411 ± 0.1310 0.1555 ± 0.1374 0.2492 ± 0.2192 0.4087 ± 0.2462 0.3520 ± 0.1804 0.1886 ± 0.1832 0.3355 ± 0.2785 0.1451 ± 0.1079 0.1619 ± 0.1363 0.2654 ± 0.2051 0.2332 ± 0.1683 0.1886 ± 0.1295 0.2863 ± 0.2308 0.1248 ± 0.1536 0.1845 ± 0.1855 0.3725 ± 0.2457 0.4604 ± 0.1846
0.3082 0.3993 0,0110 0.1149 0.0232 0.0057 0.0046 0.0543 0.0431 0.0627 0.0709 0.2326 0.2189 0.0688 0.0988 0.0727 0.0494 0.0764 0.1103 0.1002 0.0910 0.0279 0.0947 0.1699 0.3356
0.3983 0.4947 0.0230 0.1529 0.0394 0.0085 0.0070 0.0891 0.1065 0.1214 0.1751 0.3891 0.3404 0.1396 0.2517 0.1309 0.1333 0.2539 0.1875 0.1606 0.2302 0.0701 0.1263 0.2947 0.4376
0.4893 0.5954 0.0462 0.1997 0.0759 0.0136 0.0116 0.1310 0.1980 0.2227 0.3936 0.5408 0.4620 0.2282 0.5268 0.1985 0.2400 0.3965 0.2993 0.2490 0.4906 0.1582 0.1779 0.4995 0.5838
2.55 1.72 9.19 7.35 0.56 75.89 77.66 78.37 0.86 1.26 0.52 0.11 0.88 0.17 0.14 2.10 1.14 1.47 2.32 7.93 0.18 0.76 0.93 0.31 3.05
magnitude of ambient dimension, confirming that the computational cost is effectively decoupled from m in the regime m ≫ k that FAST mode is designed for. Curvature field reveals distinct geometric regimes. Beyond the timing results, the curvature scores in Table 5 reveal interpretable geometric structure that varies systematically across domain and dataset characteristics. Datasets with high mean curvature (K̄ > 0.35) are concentrated in two groups: the artificial and structured high-dimensional datasets (madelon, K̄ = 0.498; Speech, K̄ = 0.403; hepatitisC, K̄ = 0.460) and the small genomic cancer datasets (SRBCT, K̄ = 0.409; GLI, K̄ = 0.373). High mean curvature in these datasets signals that the underlying data manifold is strongly non-linear, local neighborhoods deviate substantially from their tangent planes, which is consistent with the known separability structure of cancer expression arrays and the deliberately non-linear geometry of benchmark classification datasets such as madelon. Datasets with low mean curvature (K̄ < 0.07) cluster in a qualitatively different regime: activity recognition (har, K̄ = 0.036), speech pathology (parkinson-speech-uci, K̄ = 0.066), and the two natural image benchmarks MNIST_784 (K̄ = 0.013) and Fashion-MNIST (K̄ = 0.011). These scores indicate that the manifold underlying these datasets is nearly flat at the scale of the k-nearest-neighbour patches, a result consistent with the well-established empirical finding that pixel-space representations of natural images lie near low-curvature manifolds of moderate intrinsic dimension, and that sensor-based activity signals exhibit smooth, nearly-linear trajectory structure in feature space [Fefferman et al., 2016]. A noteworthy characteristic of the MNIST datasets is their strongly right-skewed curvature distribution: the large majority of samples are concentrated in a narrow low-curvature band, as evidenced by the tight interquartile range [Q25 , Q75 ] relative to the full empirical range, while a small number of geometrically exceptional points exhibit extreme curvature values that inflate the standard deviation well above the mean (coefficient of variation CV > 1.6 for both MNIST_784 and Fashion-MNIST). In such heavy-tailed regimes, a monotone variance-stabilising transformation such as log(1 + Ki ) can be applied to the curvature scores prior to any thresholding or ranking procedure, compressing the dynamic range without altering the relative ordering of points.
27
P REPRINT - J UNE 5, 2026
Geometric heterogeneity within datasets. The spread of the curvature distribution, characterised by the coefficient of variation CV = σ/K̄ and the interquartile range [Q25 , Q75 ], provides additional geometric information beyond the mean. Datasets with high CV, notably MNIST_784 (CV = 2.21), parkinson-speech-uci (CV = 1.39), har (CV = 1.20), and GCM (CV = 1.23), exhibit a strongly heterogeneous curvature landscape: a majority of points lie near flat regions (as indicated by the low Q25 values), while a minority occupy high-curvature boundary zones. This heterogeneity is precisely the geometric signature that curvature-aware algorithms such as MCBP exploit: the highcurvature minority identifies the boundary points that separate structurally distinct regions of the manifold. Conversely, datasets with low CV, madelon (CV = 0.29), Speech (CV = 0.35), and hepatitisC (CV = 0.40), exhibit a nearly uniform curvature field with little point-to-point variation, which is indicative of globally curved manifolds where bending is distributed evenly rather than concentrated at specific boundary structures. Image and face datasets: curvature encodes visual structure. The three MNIST variants illustrate a particularly interesting comparison. MNIST_784 and Fashion-MNIST share the same ambient dimension (m = 784) and sample size (n = 70,000) and yield nearly identical mean curvatures (0.013 and 0.011, respectively), consistent with the known similarity of their geometric structures in pixel space. Kuzushiji-MNIST, by contrast, exhibits a substantially higher mean curvature (K̄ = 0.103, CV = 0.68), reflecting the greater intra-class variability and more complex stroke geometry of cursive Kuzushiji script relative to standard printed digits or fashion items. Among the face datasets, Olivetti_Faces (m = 4,096, K̄ = 0.352) shows markedly higher curvature than UMIST_Faces_Cropped (m = 10,304, K̄ = 0.145) and coil-20 (m = 1,024, K̄ = 0.155), likely because the Olivetti dataset contains multiple expressions and lighting conditions per subject, introducing discontinuities in the face manifold that manifest as elevated local curvature. Genomic expression arrays: curvature discriminates cancer subtypes. The eleven genomic datasets in Table 5 span m from 2,308 to 22,283 gene expression features, representing the extreme high-dimensional, low-sample-size (n ≪ m) regime that is characteristic of microarray and RNA-seq studies. All eleven are processed by MeCuCo in under 8 seconds. The curvature scores reveal heterogeneous geometric structure across cancer types: lymphoma subtypes (DLBCL, MLL) and renal datasets (AP_Omentum_Kidney) show moderate mean curvature with high CV, consistent with the presence of a small number of geometrically extreme samples (potential outliers or boundary cases between subtypes), while brain tumour datasets (GLI, K̄ = 0.373) show higher curvature concentrated at class boundaries. These observations are consistent with the known biological heterogeneity of these datasets and suggest that curvature-based analysis could complement existing subtype discovery pipelines in high-dimensional genomic spaces, a direction we leave for future work. Summary. The results of Table 5 demonstrate that MeCuCo extends the reach of local mean curvature estimation to the very-high-dimensional regime that characterises modern machine learning applications, processing datasets with up to 54,621 features in seconds while yielding geometrically interpretable and domain-coherent curvature scores. The computational cost scales primarily with the sample size n rather than the ambient dimension m, confirming the theoretical prediction that FAST mode decouples the per-point cost from m when m ≫ k.
6
Conclusions
This paper addressed a fundamental computational bottleneck in geometric machine learning: the estimation of local mean curvature at every point of a high-dimensional dataset. The original formulation of the MCBP curvature estimator 2 incurs an O(m4 ) cost per point through the explicit construction of the feature matrix H ∈ Rm×O(m ) and the subsequent product HH ⊤ , rendering it computationally intractable for datasets with more than approximately 300 features. We introduced the MeCuCo algorithm (Mean Curvature Computation), a mathematically exact reformulation and a principled approximation that together extend the reach of local mean curvature estimation to the very-highdimensional regime characteristic of modern machine learning applications. The first contribution is an exact algebraic identity derived from the orthogonality of the local eigenvector frame. ⊤ By showing that HH ⊤ = 12 W (2) W (2) + 12 Im , a consequence of the single identity W W ⊤ = Im , we proved that the mean curvature estimator eliminates the O(m4 ) tensor contraction entirely and reducing the per-point cost, after eigendecomposition, to O(m2 ). This identity is exact, introduces no approximation error beyond floating-point rounding, and applies to any symmetric positive semi-definite local covariance matrix regardless of dataset characteristics. The second contribution exploits the low-rank structure of the local covariance matrix, whose rank is at most k − 1 ≪ m whenever D > k, to replace the O(m3 ) full eigendecomposition with a truncated SVD of the k × m centred neighborhood matrix at cost O(k 2 m). The contribution of the m − p null-space eigenvectors, which cannot be
28
P REPRINT - J UNE 5, 2026
recovered from the truncated SVD alone, is handled through an analytical approximation grounded in the expected outer product of random orthonormal null-space bases under the Haar measure over the Stiefel manifold. The resulting FAST mode estimator has total per-point cost O(k 2 m + kmp2 ), independent of m3 or m4 , and its approximation error decreases monotonically as m/k grows, falling below 2% for m/k ≥ 56 and below 1% for m/k ≥ 100, exactly the regime in which MeCuCo is most needed. Empirical evaluation on more than 40 publicly available OpenML datasets confirmed that these theoretical gains translate to practice. In the benchmark corpus, MeCuCo achieves a median Spearman rank correlation of ρS = 0.9998 and a median normalised MAE below 0.001 relative to the original formulation, with speedups reaching 298× on highdimensional datasets such as arrhythmia (m = 279) and 144× on USPS (m = 256). On the very-high-dimensional corpus, spanning ambient dimensions from m = 400 to m = 54,621, all inaccessible to the original method, MeCuCo computed local mean curvatures for all 25 datasets in times ranging from 0.11 seconds for a small genomic array to under 3.5 minutes for the three large-scale image benchmarks with n = 70,000 samples. Crucially, the dominant factor governing running time was found to be the sample size n rather than the ambient dimension m, confirming that FAST mode effectively decouples the computational cost from the dimensionality. Beyond computational efficiency, the curvature scores produced by MeCuCo revealed geometrically interpretable structure across all application domains examined: near-flat manifolds in activity recognition and natural image datasets, concentrated boundary curvature in cancer expression arrays, and globally curved manifolds in structured benchmark datasets. These observations reinforce the theoretical grounding of the manifold hypothesis as a working model for high-dimensional data, and demonstrate that local mean curvature, previously accessible only for low-dimensional datasets, can now serve as a practical, scalable descriptor for geometric data analysis at any dimension. Limitations. Two limitations of the present work deserve acknowledgement. First, the FAST mode approximation rests on the assumption that the null-space eigenvectors of the local covariance are effectively random in the ambient space, which holds when m ≫ k but degrades when the null space is low-dimensional (m/k ≈ 2–5). In such regimes, the E XACT mode should be preferred. Second, both estimators assume that the data lie near a smooth Riemannian manifold, and neither provides meaningful curvature estimates for purely categorical or heavily discrete data, as illustrated by the car-evaluation dataset; appropriate continuous embeddings should be applied prior to curvature estimation in such cases. Future Work. Several directions for future investigation emerge naturally from this work. On the theoretical side, a rigorous bound on the bias of the FAST mode approximation, beyond the empirical O(k/m) decay observed in the experiments, would strengthen the theoretical guarantees of the method, as would a formal consistency analysis establishing the rate at which the discrete curvature estimator converges to its smooth-manifold counterpart as n → ∞ and k → ∞ at an appropriate rate. On the algorithmic side, the truncated SVD in FAST mode could be replaced by randomised low-rank approximations, which achieve similar accuracy with lower constant factors and are amenable to streaming and distributed computation, an important consideration for datasets whose sample size n prevents loading the full neighborhood matrix into memory. GPU-accelerated implementations of the per-point SVD and the Gtens einsum contraction are a natural extension that could reduce the running time on large-n datasets by an additional order of magnitude. On the applications side, the availability of scalable local mean curvature opens a broad research agenda for curvatureaware machine learning. In dimensionality reduction, curvature-weighted objectives, penalising embeddings that flatten high-curvature boundary regions, could yield representations that preserve the geometric structure of cluster boundaries more faithfully than variance-based criteria; this is a direct extension of the curvature-aware manifold learning paradigm. In clustering, incorporating local mean curvature as an additional feature or as a density-modulating weight in algorithms such as DBSCAN or spectral clustering could improve boundary delineation in datasets where density-based criteria alone are insufficient to resolve geometrically complex interfaces between classes. In supervised classification and active learning, the curvature scores produced by MeCuCo provide a principled, geometry-driven criterion for identifying the most informative samples to label: points near high-curvature decision boundaries are precisely those where the classifier’s confidence is lowest and where additional labels yield the greatest reduction in generalisation error. In deep learning, curvature-based regularisation of the feature manifold, penalising layers that produce representations with excessive local curvature, could serve as a geometric alternative to standard regularisers such as weight decay or dropout, providing an inductive bias that encourages smoother intermediate representations and potentially improving generalisation in overparameterised regimes. Finally, in graph neural networks and geometric deep learning, node-level or edge-level curvature estimates derived from MeCuCo could replace or complement discrete Ricci curvature as an input feature or structural descriptor, enriching the geometric signal available to message-passing architectures and potentially improving performance on node classification, link prediction, and anomaly detection tasks. Taken together, these directions suggest that scalable local mean curvature estimation, as made practical by 29
P REPRINT - J UNE 5, 2026
MeCuCo, may serve as a unifying geometric primitive across the full pipeline of modern machine learning, from raw data preprocessing through representation learning to downstream prediction and decision-making.
Statements and declarations Funding This work has been supported by CNPq (National Council for Scientific and Technological Development) through grant number 301432/2025-2. This study was also financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. Code availability Python scripts to reproduce the results reported in this paper may be found at https://github.com/ alexandrelevada/MeCuCo. Data availability All datasets used in the experiments are publicly available at www.openml.org.
References Melanie Weber. Geometric machine learning. AI Magazine, 46(1):e12210, 2025. doi: https://doi.org/10.1002/aaai. 12210. Mathilde Papillon, Sophia Sanborn, Johan Mathe, Louisa Cornelis, Abby Bertics, Domas Buracas, Hansen J. Lillemark, Christian Shewmake, Fatih Dinc, Xavier Pennec, and Nina Miolane. Beyond Euclid: An illustrated guide to modern machine learning with geometric, topological, and algebraic structures. Machine Learning: Science and Technology, 6, 2025. doi: 10.1088/2632-2153/adf375. Michael M. Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep learning: Going beyond Euclidean data. IEEE Signal Processing Magazine, 34(4):18–42, 2017. doi: 10.1109/MSP.2017. 2693418. Michael M. Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, Groups, Graphs, Geodesics, and Gauges. arXiv preprint arXiv:2104.13478, 2021. doi: 10.48550/arXiv.2104.13478. URL https://arxiv.org/abs/2104.13478. Charles Fefferman, Sanjoy Mitter, and Hariharan Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 29(4):983–1049, 2016. doi: 10.1090/jams/852. Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319–2323, 2000. doi: 10.1126/science.290.5500.2319. Sam T. Roweis and Lawrence K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290 (5500):2323–2326, 2000. doi: 10.1126/science.290.5500.2323. Laurens van der Maaten and Geoffrey Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research, 9: 2579–2605, 2008. URL http://www.jmlr.org/papers/v9/vandermaaten08a.html. Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction, 2018a. URL http://arxiv.org/abs/1802.03426. cite arxiv:1802.03426Comment: Reference implementation available at http://github.com/lmcinnes/umap. Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger. Umap: Uniform manifold approximation and projection. Journal of Open Source Software, 3(29):861, 2018b. doi: 10.21105/joss.00861. URL https: //doi.org/10.21105/joss.00861. Yasuhiko Asao and Yuichi Ike. Curvature of point clouds through principal component analysis. arXiv preprint arXiv:2106.09972, 2022. doi: 10.48550/arXiv.2106.09972. URL https://arxiv.org/abs/2106.09972. Ming-Yen Cheng and Hau-Tieng Wu. Efficient Weingarten map and curvature estimation on manifolds. Machine Learning, 110(6):1467–1510, 2021. doi: 10.1007/s10994-021-05953-4. A. L. M. Levada. A mean curvature approach to boundary detection: Geometric insights for unsupervised learning. arXiv preprint arXiv:2605.04274, 2025. doi: 10.48550/arXiv.2605.04274. URL https://arxiv.org/abs/2605.04274. 30
P REPRINT - J UNE 5, 2026
Karish Grover, Geoffrey J. Gordon, and Christos Faloutsos. Curvgad: Leveraging curvature for enhanced graph anomaly detection. In Proceedings of the 42nd International Conference on Machine Learning, Proceedings of Machine Learning Research. PMLR, 2025. URL https://proceedings.mlr.press/v267/grover25a.html. Jia Guo, Shuai Lu, Weihang Zhang, Fang Chen, Huiqi Li, and Hongen Liao. Dinomaly: The less is more philosophy in multi-class unsupervised anomaly detection. In 2025 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 20405–20415, 2025. doi: 10.1109/CVPR52734.2025.01900. Yusuke Imoto, Tomonori Nakamura, Emerson G. Escolar, Michio Yoshiwaki, Yoji Kojima, Yukihiro Yabuta, Yoshitaka Katou, Takuya Yamamoto, Yasuaki Hiraoka, and Mitinori Saitou. Resolution of the curse of dimensionality in singlecell RNA sequencing data analysis. Life Science Alliance, 5(12):e202201591, 2022. doi: 10.26508/lsa.202201591. Sintija Stevanoska, Jurica Levatić, and Sašo Džeroski. Semi-supervised learning from tabular data with autoencoders: when does it work? Machine Learning, 114(11), October 2025. ISSN 0885-6125. doi: 10.1007/s10994-025-06898-8. Alaa Tharwat and Wolfram Schenck. A survey on active learning: State-of-the-art, practical challenges and research directions. Mathematics, 11(4), 2023. ISSN 2227-7390. doi: 10.3390/math11040820. URL https://www.mdpi. com/2227-7390/11/4/820. Michael Spivak. A Comprehensive Introduction to Differential Geometry. Publish or Perish, 1999. Manfredo P. do Carmo. Differential Geometry of Curves and Surfaces: Revised and Updated Second Edition. Dover Publications, New York, UNITED STATES, 2016. ISBN 978-0-486-81797-2. Barrett O’Neill. Elementary Differential Geometry. Academic Press, 2nd edition, 2006. Loring W Tu. Differential geometry: connections, curvature, and characteristic classes, volume 275. Springer, 2017. John Oprea. Differential Geometry and its Applications. The Mathematical Association of America, New York, 2 edition, 2007. Tristan Needham. Visual Differential Geometry and Forms: A Mathematical Drama in Five Acts. Princeton University Press, Princeton, NJ, 2021. ISBN 978-0-691-20369-0. A. Singer and H.-T. Wu. Vector diffusion maps and the connection laplacian. Communications on Pure and Applied Mathematics, 65(8):1067–1144, 2012. doi: https://doi.org/10.1002/cpa.21395. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. ISBN 0-89871-447-8 (paperback). Elizabeth S. Meckes. The Random Matrix Theory of the Classical Compact Groups. Cambridge Tracts in Mathematics. Cambridge University Press, 2019. Joaquin Vanschoren, Jan N. van Rijn, Bernd Bischl, and Luis Torgo. OpenML: Networked science in machine learning. SIGKDD Explorations Newsletter, 15(2):49–60, 2013. doi: 10.1145/2641190.2641198. Bernd Bischl, Giuseppe Casalicchio, Matthias Feurer, Pieter Gijsbers, Frank Hutter, Michel Lang, Rafael G. Mantovani, Jan N. van Rijn, and Joaquin Vanschoren. OpenML benchmarking suites. In Advances in Neural Information Processing Systems (NeurIPS) Track on Datasets and Benchmarks, 2021. URL https://arxiv.org/abs/1708. 03731. Sourav Chatterjee. A new coefficient of correlation. Journal of the American Statistical Association, 116(536): 2009–2022, 2021. doi: 10.1080/01621459.2020.1758115.
31