Paper deep dive
Shape Operator PCA: Curvature-Aware Projections for Geometric Machine Learning
Alexandre L. M. Levada
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/18/2026, 6:12:08 AM
Summary
The paper introduces SHOPCA (Shape Operator-based Principal Component Analysis), a novel unsupervised metric learning and dimensionality reduction method. SHOPCA enhances classical PCA by regularizing the global covariance matrix with a mean shape operator derived from differential geometry, incorporating curvature information alongside variance. A single trace-normalized mixing coefficient alpha controls the trade-off between variance and curvature, with an unsupervised eigengap criterion used to select alpha. The method is evaluated on over 50 datasets, showing improved clustering quality over PCA and better performance than UMAP in small-sample settings.
Entities (10)
Relation Signals (6)
Alexandre Luis Magalhães Levada → affiliatedwith → Federal University of São Carlos
confidence 95% · Affiliation: Federal University of São Carlos
SHOPCA → extends → PCA
confidence 95% · SHOPCA regularizes the global covariance matrix... recovering standard PCA at alpha = 0
SHOPCA → improves → clustering quality
confidence 95% · SHOPCA consistently improves clustering quality over PCA across a broad range of datasets
SHOPCA → uses → Shape Operator
confidence 95% · regularizes the global covariance matrix using the mean shape operator
SHOPCA → evaluatedon → Benchmark Datasets
confidence 90% · We evaluate SHOPCA on more than 50 real-world benchmark datasets
SHOPCA → outperforms → UMAP
confidence 90% · surpasses UMAP on small-sample settings
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:In this paper, we propose SHOPCA (Shape Operator-based Principal Component Analysis), a novel method for unsupervised metric learning and dimensionality reduction that incorporates differential geometric information into the covariance structure of classical PCA. SHOPCA regularizes the global covariance matrix using the mean shape operator, defined as the average of the absolute local shape operators estimated from the data manifold, steering principal components toward directions of both maximum variance and informative curvature. A single trace-normalized mixing coefficient $\alpha$ controls the regularization, recovering standard PCA at $\alpha = 0$ and a curvature-driven embedding as $\alpha \to \infty$. We further introduce a fully unsupervised criterion for selecting $\alpha$ based on the spectral eigengap of the regularized covariance matrix, maximizing the relative separation between the top-$d$ and remaining eigenvalues without using class labels. We evaluate SHOPCA on more than 50 real-world benchmark datasets, comparing it with PCA, ISOMAP, and UMAP using Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), Fowlkes-Mallows index (FM), and V-measure. Results show that SHOPCA consistently improves clustering quality over PCA across a broad range of datasets and surpasses UMAP on small-sample settings, where iterative neighborhood-based manifold estimation can degrade. SHOPCA is computationally tractable, parameter-efficient, and applicable to domains requiring fully unsupervised, geometry-aware dimensionality reduction.
Tags
Links
- Source: https://arxiv.org/abs/2608.15313v1
- Canonical: https://arxiv.org/abs/2608.15313v1
Trouble viewing inline? Open PDF directly →
Full Text
107,684 characters extracted from source content.
Expand or collapse full text
Shape Operator PCA: Curvature-Aware Projections for Geometric Machine Learning Alexandre Luis Magalhães Levada Affiliation: Federal University of São Carlos Affiliation: 13565-905, São Carlos-SP, Brazil Email: alexandre.levada@ufscar.br Abstract In this paper, we propose SHOPCA (Shape Operator-based Principal Component Analysis), a novel method for unsupervised metric learning and dimensionality reduction that incorporates differential geometric information into the covariance structure of classical PCA. The method regularizes the global covariance matrix using the mean shape operator, defined as the average of the absolute local shape operators estimated from the data manifold, steering the principal components toward directions of both maximum variance and informative curvature and yielding a richer, more discriminative feature space. The regularization is governed by a single, trace-normalized mixing coefficient α, which recovers standard PCA at α=0α=0 and a purely curvature-driven embedding at α→∞α→∞, making SHOPCA a natural, interpretable, and cross-dataset-comparable generalization of the classical method. We introduce a fully unsupervised selection criterion based on the spectral eigengap of the regularized covariance matrix: α is chosen to maximize the relative separation between the top-d and remaining eigenvalues, identifying the most well-defined d-dimensional invariant subspace without ever consulting class labels (label-free strategy). We evaluate SHOPCA on more than 50 real-world benchmark datasets, comparing against PCA, ISOMAP and UMAP using four external clustering indices: Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), Fowlkes-Mallows index (FM), and V-measure. Experimental results show that SHOPCA consistently improves clustering quality over PCA across a broad range of datasets, and surpasses UMAP on small-sample settings where iterative neighborhood-based manifold estimation is known to degrade. The proposed approach is computationally tractable, parameter-efficient, and applicable to any domain where a fully unsupervised, geometry-aware metric is desirable. 1 Introduction The ability to measure similarity between data points is central to a wide range of machine learning and pattern recognition tasks. Algorithms such as k-nearest neighbor classifiers, hierarchical clustering, and spectral methods are all governed by the geometry of the input space, making the choice of distance metric a foundational design decision rather than an implementation detail 18, 2. This observation motivates the field of metric learning, which seeks to automatically infer a task-appropriate geometry from data itself 17. A key insight, often underappreciated, is that metric learning and dimensionality reduction are not independent endeavors but two complementary perspectives on the same underlying problem. Dimensionality reduction and metric learning are closely intertwined: most metric learning algorithms rely on PCA as a pre-processing step to achieve tractability and robustness to noise 16, yet there is no principled reason to assume that a variance-maximizing projection is optimal for the metric to be learned. Harandi et al. 13 argue that treating these two steps separately is fundamentally limiting and propose instead a Riemannian framework that jointly learns a dimensionality-reducing mapping and a metric in the induced space 34, 12, achieving higher accuracy than state-of-the-art metric learning algorithms while working directly on high-dimensional features. This unification is further supported by the observation that spectral dimensionality reduction methods can be uniformly interpreted as instances of kernel PCA with different kernels 30, establishing a formal bridge between manifold-based projection and kernel metric learning 21. Taken together, these results suggest that the most informative low-dimensional representation is one in which the projection and the metric are co-designed, a principle that motivates the geometry-aware regularization strategy proposed in the present work. This co-design principle finds its most general expression in the recent movement toward Geometric Machine Learning (GML) 32, which extends learning models beyond flat Euclidean spaces to the richer structures offered by non-Euclidean geometry, topology, and algebra 4, 3, 7, 24. The foundational work of Bronstein et al. 4 shows that the most successful deep architectures can be understood as instances of a single geometric blueprint grounded in symmetry and invariance, and subsequent work has generalized structured learning to non-Euclidean domains such as graphs and manifolds 24. Although GML has been predominantly instantiated through deep, iterative architectures, its core motivation applies equally to classical multivariate statistics: if data lies on a curved manifold, any method that ignores this curvature – including PCA, the workhorse pre-processing step of the metric-learning pipelines discussed above – operates on a flawed geometric premise. Despite this convergence of dimensionality reduction, metric learning, and geometric machine learning around the idea that projection and metric should be co-designed from the data’s intrinsic geometry, three practical gaps have limited its adoption as a routine, general-purpose tool. First, GML-inspired architectures typically require deep or iterative optimization 4, 3, which is computationally heavy and poorly suited to the small-sample, high-dimensional regimes common in bioinformatics, spectroscopy, and other scientific domains. Second, joint metric-learning and dimensionality-reduction frameworks such as the Riemannian approach of Harandi et al. 13 learn their geometry from labeled data, which precludes their use in purely exploratory or unsupervised analysis, precisely the setting in which PCA is most widely used in practice. Third, neighborhood-graph-based nonlinear embeddings such as UMAP 21 construct their geometric representation from a k-nearest-neighbor graph, and this construction is known to become unreliable when the number of samples is small, since a sparse, unstable neighborhood graph no longer faithfully represents the underlying manifold. At the other extreme, classical linear methods such as PCA remain entirely agnostic to the manifold’s intrinsic curvature, retaining only its variance structure 16. To the best of our knowledge, no existing method incorporates second-order differential-geometric information, curvature, as opposed to the first-order, variance-only geometry captured by PCA, into a linear, closed-form, and fully unsupervised metric-learning pipeline, nor does any such method extend its label-free character to the selection of its own regularization hyperparameter, which is typically tuned against labeled validation data or left to ad hoc heuristics with no geometric justification. This three-part gap, computational cost, reliance on labels, and instability in small-sample regimes, motivates the present work. Motivated by the geometric perspective outlined above, we propose SHOPCA (Shape Operator-based Principal Component Analysis), a novel method that bridges classical multivariate statistics and differential geometry for unsupervised metric learning and dimensionality reduction, addressing each of the three gaps identified above within a single, linear, closed-form procedure. The central idea is to enrich the global covariance structure of PCA with curvature information derived from the data manifold. To this end, we estimate, for each data point, a local shape operator (a classical object from differential geometry that encodes how the manifold curves in the ambient space) and aggregate these local estimates into a mean shape operator. This mean shape operator is then used to regularize the global covariance matrix, yielding a geometry-aware covariance Σcurv _curv. The principal components of Σcurv _curv are therefore steered toward directions of both maximum variance and informative curvature, producing a more discriminative feature space than standard PCA without requiring any labeled data. The main contributions of the proposed method are fivefold: 1. a principled geometric regularization of the PCA covariance matrix via the mean shape operator, establishing a formal connection between classical multivariate statistics and differential geometry, and positioning it within the broader GML movement as a linear, closed-form alternative to deep geometric architectures; 2. a natural, trace-normalized generalization of classical PCA, recovered exactly as standard PCA at α=0α=0 and yielding a purely curvature-driven embedding at α→∞α→∞, with a single and cross-dataset-comparable mixing coefficient controlling the geometry-variance trade-off; 3. a fully unsupervised, label-free criterion for selecting α, based on the spectral eigengap of Σcurv(α) _curv(α), which closes the one remaining supervised step of the original pipeline and renders the entire SHOPCA procedure, representation learning and hyperparameter selection alike, independent of class labels; 4. robustness to small sample sizes, addressing a known limitation of neighborhood-graph-based nonlinear methods such as t-SNE and UMAP, whose local graph estimation degrades when few samples are available; and 5. a computationally tractable, closed-form linear projection that requires no iterative optimization, no explicit manifold parametrization, and no neighborhood-graph construction beyond the local k-N neighborhoods used to estimate curvature. The remainder of the paper is organized as follows: Section 2 describes the proposed shape-operator-based PCA method, its trace-normalized regularization, and the unsupervised eigengap criterion for hyperparameter selection in detail. Section 3 presents the computational experiments and the obtained results. Lastly, Section 4 presents our conclusions and final remarks. 2 Differential Geometry Basics Classical PCA characterizes a dataset through a single first-order statistic, the covariance matrix, which captures how the data spreads but is entirely blind to how it bends. Two datasets can share an identical covariance structure while lying on manifolds of markedly different curvature, and it is precisely along directions of high curvature – where a manifold folds, twists, or separates into distinct branches – that class boundaries and cluster structure are often most pronounced. Formalizing this notion of “bending” requires a small set of classical tools from differential geometry: the manifold itself, its tangent space, and the first and second fundamental forms, which together give rise to the shape operator at the heart of the proposed regularization. We introduce these objects below, following the standard treatment of 9, 20, 8, and indicate throughout how each formal object is subsequently approximated from finite samples. Manifold. An m-dimensional smooth manifold is a topological space that is locally Euclidean of dimension m: every point admits a neighborhood homeomorphic to an open subset of ℝmR^m via a chart, and the transition maps between overlapping charts are smooth 20. Throughout this paper we consider a manifold ℳM embedded in the ambient space via a smooth embedding ℳ↪ℝdM ^d, of intrinsic dimension m≪dm d – the standard manifold hypothesis underlying most dimensionality-reduction and manifold-learning methods, under which the observed high-dimensional data =1,…,n⊂ℝdX=\x_1,…,x_n\ ^d are assumed to be samples drawn from such an ℳM. Tangent space. At each point ∈ℳx , the tangent space TℳT_xM is the m-dimensional vector space of vectors tangent to ℳM at x, equivalently the space spanned by the velocity vectors of smooth curves on ℳM through x 9, 20. It is the best linear (first-order) approximation of ℳM at x. This is the object directly targeted by the local estimation procedure of Section 2.2: the leading eigenvectors of the local sample covariance matrix computed from the k-nearest neighbors of ix_i furnish an empirical orthonormal basis for TiℳT_x_iM, an approximation strategy also used, in various forms, throughout the manifold-learning literature. First fundamental form (metric tensor). Given an orthonormal basis 1,…,m\w_1,…,w_m\ of TℳT_xM, the first fundamental form is the bilinear form obtained by restricting the ambient Euclidean inner product to the tangent space, I(,)=⟨,⟩,,∈Tℳ,I_x(v,w)= ,w , ,w∈ T_xM, (1) represented, relative to the chosen basis, by the metric tensor g=(gij)g_x=(g_ij) with gij=⟨i,j⟩g_ij= _i,w_j 9, 8. The first fundamental form is intrinsic, it determines all metric properties of ℳM (lengths, angles, volumes) measurable without reference to the ambient space, and is precisely the object that metric learning seeks to adapt to a given task 17, 2. In the proposed method, the local sample covariance matrix iC_i plays the role of an estimate of the (inverse) metric tensor at ix_i, following the classical identification of a local Gaussian covariance with the Mahalanobis metric it induces. Second fundamental form. Whereas the first fundamental form is intrinsic, the second fundamental form is extrinsic: it measures how ℳM curves away from its own tangent space as one moves through the ambient space ℝdR^d. For a hypersurface ℳ↪ℝdM ^d (codimension one, admitting a smooth unit normal field :ℳ→d−1n:M ^d-1 up to sign), the second fundamental form at x is the symmetric bilinear form I(,)=⟨−∇,⟩,I_x(v,w)= - _vn,\,w , (2) where ∇ _v denotes the ambient directional derivative along v 9, 23. Intuitively, II_x quantifies the rate at which ℳM pulls away from the tangent plane TℳT_xM along a given direction, and is the source of all extrinsic curvature information used in the present work; it is this quantity that the quadratic and cross-product terms of the local PCA basis (Section 2.2) are designed to approximate. Shape operator. The shape operator, or Weingarten map, at ∈ℳx is the linear endomorphism of the tangent space S:Tℳ→Tℳ,S()=−∇,S_x:T_xM→ T_xM, S_x(v)=- _vn, (3) i.e., minus the derivative of the unit normal field along v 9. It is self-adjoint with respect to the first fundamental form, I(S(),)=I(,S())I_x(S_x(v),w)=I_x(v,S_x(w)), and therefore admits an orthonormal eigenbasis with real eigenvalues κ1,…,κm _1,…, _m, the principal curvatures, whose eigenvectors are the principal directions along which ℳM bends maximally and minimally 9, 23. The shape operator relates the two fundamental forms via I(,)=I(S(),)I_x(v,w)=I_x(S_x(v),w), or, in matrix form, S=I−1IIS_x=I_x^-1\,I_x, exactly the identity used to define the local shape operator estimator in Eq. (9), with I−1I_x^-1 approximated by the local covariance matrix iC_i. From the eigenvalues of S_x, two classical scalar invariants follow directly: the mean curvature H=1m∑jκj=1mtr(S)H= 1m _j _j= 1mtr(S_x), an extrinsic average-bending measure, and the Gaussian curvature K=∏jκj=det(S)K= _j _j= (S_x), which – despite being defined through the extrinsic shape operator – is in fact intrinsic to ℳM by Gauss’s Theorema Egregium 9. The shape operator thus offers a compact, basis-independent, and complete local description of manifold curvature, unifying the intrinsic (first fundamental form) and extrinsic (second fundamental form) geometry of ℳM into a single linear operator at each point. Section 2.2 develops a finite-sample estimator of S_x from a local k-nearest-neighbor neighborhood, and Section 2 uses the resulting mean shape operator S¯ S to regularize the global PCA covariance matrix, steering its principal components toward directions of both maximum variance and informative curvature. We remark that the classical hypersurface formulation above, with a single normal direction, is adopted here for notational clarity; its extension to general codimension d−m>1d-m>1, where the second fundamental form is valued in the full normal bundle rather than along a single normal vector 8, is handled operationally in Section 2.2 through the local PCA frame, which spans the tangent directions directly and treats the remaining directions collectively via the quadratic and cross-product terms of the local basis. 2.1 Approximating the Shape Operator The shape operator is a central concept in differential geometry, providing an intrinsic and compact characterization of curvature through the relationship between the first and second fundamental forms of a smooth manifold. Intuitively, it quantifies how the unit normal vector field varies as one moves along directions in the tangent space, thereby encoding the local bending behavior of the manifold. Formally, let ℳ↪ℝdM ^d be a smooth hypersurface and let :ℳ→d−1n:M ^d-1 denote the unit normal vector field. The shape operator at a point ∈ℳx is the linear map :Tℳ→TℳS_x:T_xM→ T_xM defined by ()=−∇,S_x(v)=- _v\,n, (4) where ∇ _v denotes the covariant derivative of n in the direction ∈Tℳv∈ T_xM. This operator is self-adjoint with respect to the Riemannian metric, and its spectral decomposition yields real eigenvalues κ1,κ2,…,κd−1 _1, _2,…, _d-1, known as the principal curvatures, which measure the intensity of bending along specific tangential directions, while the associated eigenvectors define the principal directions along which these curvatures are attained. From these eigenvalues, two fundamental scalar invariants are derived: the mean curvature H=1d−1∑j=1d−1κjH= 1d-1 _j=1^d-1 _j (5) and the Gaussian curvature K=∏j=1d−1κj,K= _j=1^d-1 _j, (6) which provide complementary summaries of the geometric behavior of ℳM. This decomposition simultaneously captures both the magnitude and orientation of local geometric variation, making the shape operator a natural and principled tool for incorporating manifold geometry into data-driven learning methods. Figure 1: The shape operator measures how a normal vector changes from the tail to the tip of a tangent vector v→ v. In other words, it tells us how fast the normal changes, or roughly speaking, how fast a surface turns if you roll it along the floor in some direction. 2.2 Computing the Local Shape Operators Let =1,2,…,n⊂ℝdX=\x_1,x_2,…,x_n\ ^d be the dataset and let k be a fixed number of nearest neighbors. The computation of the mean shape operator ¯ S proceeds in four steps for each point i∈x_i . Step 1: k-Nearest Neighbor Search. For each point ix_i, we identify its k nearest neighbors in the Euclidean sense: i=kNN(i,k)=i1,i2,…,ik⊂,N_i=kNN(x_i,\,k)= \x_i_1,x_i_2,…,x_i_k \ , (7) where the neighborhood iN_i provides a local sample of the data manifold ℳM in the vicinity of ix_i. Step 2: Local Covariance Matrix and Metric Tensor. Given the local neighborhood iN_i, we compute the d×d× d local covariance matrix i=1k−1∑∈i(−¯i)(−¯i)⊤,C_i= 1k-1 _x _i (x- x_i ) (x- x_i )^\! , (8) where ¯i=1k∑∈i x_i= 1k _x _ix is the local mean. The local metric tensor, which approximates the first fundamental form of ℳM at ix_i, is identified as i≈i−1g_i _i^-1. This identification is motivated by the fact that the inverse covariance matrix in the Mahalanobis distance encodes the directions and magnitudes of spread of the manifold in the ambient space, thereby approximating the pullback metric induced by the embedding ℳ↪ℝdM ^d. Step 3: Local PCA Basis and Second Fundamental Form. Let i=iii⊤C_i=W_i\, _i\,W_i be the eigendecomposition of the local covariance matrix, where the columns 1(i),2(i),…,d(i)w_1^(i),w_2^(i),…,w_d^(i) of iW_i are the local principal directions sorted in decreasing order of variance. These eigenvectors define a local orthonormal frame that approximates the tangent space of ℳM at ix_i. We approximate the second fundamental form using quadratic and cross-product terms constructed from the local PCA basis 10. Specifically, we form the matrix i∈ℝd×pQ_i ^d× p, where p=1+d+d(d+1)/2p=1+d+d(d+1)/2, whose columns comprise: 1. a constant column ∈ℝd1 ^d; 2. the d principal directions j(i)w_j^(i), for j=1,…,dj=1,…,d; 3. the d squared terms (j(i))∘2 (w_j^(i) )^\! 2, where ∘2 2 denotes the elementwise square; and 4. the (d2) d2 cross-product terms j(i)∘l(i)w_j^(i) _l^(i) for all j<lj<l, where ∘ denotes the elementwise (Hadamard) product. The quadratic and cross-product columns of iQ_i collectively form the matrix i∈ℝd×rH_i ^d× r, with r=d+(d2)r=d+ d2, obtained by discarding the first d+1d+1 columns of iQ_i. The second fundamental form is then approximated as i=ii⊤∈ℝd×d.I_i=H_i\,H_i \;∈\;R^d× d. (9) Intuitively, iI_i captures the bending of ℳM at ix_i by measuring how the second-order structure of the local PCA frame deviates from a flat configuration. Step 4: Local Shape Operator and Mean Shape Operator. The local shape operator (Weingarten map) at ix_i is defined as the linear map relating the second fundamental form to the first fundamental form. Using the approximations established in Steps 2 and 3, the local shape operator is given by i=−ii−1=−ii∈ℝd×d,S_i=-\,I_i\,g_i^-1=-\,I_i\,C_i\;∈\;R^d× d, (10) where the negative sign follows the standard convention in differential geometry. Finally, the mean shape operator is obtained by averaging the absolute values of the local shape operators over all n data points: ¯=1n∑i=1ni=−1n∑i=1nii, S= 1n _i=1^nS_i=- 1n _i=1^nI_i\,C_i, (11) The matrix ¯ S aggregates local curvature information from the entire dataset into a single d×d× d symmetric matrix, which is subsequently used to regularize the global covariance matrix. The complete pseudocode for the proposed average shape operator estimation is presented in Algorithm 1. Algorithm 1 Average Shape Operator Estimation 1: X∈ℝn×dX ^n× d: data matrix (n samples, d features), k∈ℕk : number of neighbors. 2: S∈ℝdS ^d: average shape operator. 3: function Shape-Operator(X,kX,k) 4: for i←1i← 1 to n do 5: Pi←P_i← Nearest_Neighbors(x→i,k x_i,k) 6: Σi← _i← Covariance(PiP_i) 7: U←U← Eigenvectors(Σi _i) 8: Build matrix i∈ℝd×pQ_i ^d× p 9: Build matrix HiH_i (last d+(d2)d+ d2 columns of iQ_i) 10: i←HiHiTI_i← H_iH_i^T 11: i←−iΣiS_i←-I_i _i 12: S←S+SiS← S+S_i 13: end for 14: return S/nS/n 15: end function A note on the computational complexity of the proposed average shape operator estimation algorithm: finding the k nearest neighbors in a d-dimensional space can be done in O(ndlogn)O(nd n) with a KD-tree. For each sample, the computation of the local covariance matrix is O(kd2)O(kd^2) and its eigendecomposition is O(d3)O(d^3). The multiplication of HiHiTH_iH_i^T has cost O(d2r)O(d^2r), where r=d(d+1)/2r=d(d+1)/2, which leads to a total cost of O(d4)O(d^4). Therefore, the overall complexity for the estimation of the average shape operator is O(ndlogn+nd2(k+d2))O(nd n+nd^2(k+d^2)), which means that the proposed method is more sensitive to an increase in the number of features d than to an increase in the number of samples n. In case of large d, it is possible to reduce the number of features with regular PCA prior to the estimation of the average shape operator. 2.3 Vectorized Implementation of the Average Shape Operator Algorithm 1 computes the average shape operator S¯ S by looping over the n samples one at a time: for each ix_i, a local covariance matrix is formed, eigendecomposed, and combined with a quadratic/cross-product expansion of its eigenvectors to produce a single local shape operator SiS_i, which is then accumulated into the running sum. Every step within one iteration is individually cheap, O(kd2)O(kd^2) for the covariance, O(d3)O(d^3) for the eigendecomposition, O(d4)O(d^4) for the dominant matrix product, but the loop itself is executed at the Python interpreter level, and each of these steps is dispatched as its own small linear-algebra call. For matrices of moderate size d, the fixed overhead of dispatching n such calls (interpreter bookkeeping, NumPy ufunc setup, and the inability of a single small matrix operation to saturate a modern CPU’s cache hierarchy or SIMD lanes) frequently exceeds the cost of the arithmetic itself, so wall-clock time scales far worse than the operation count alone would suggest. We remove this bottleneck by reformulating each step of Algorithm 1 as a single batched tensor operation over groups of b samples at a time, rather than one sample at a time. Given a batch of indices B⊂1,…,nB⊂\1,…,n\, |B|=b|B|=b: 1. Batched local covariances. The k-nearest-neighbor sets for the whole batch are gathered into a single tensor of shape (b,k,d)(b,k,d), centered by their per-point local means, and the b local covariance matrices ii∈B\C_i\_i∈ B are obtained in one call via a batched matrix product i=1k−1Δi⊤ΔiC_i= 1k-1 _i _i, dispatched as a single (b,d,k)×(b,k,d) batched GEMM rather than b separate calls to a covariance routine. 2. Batched eigendecomposition. Since every iC_i is symmetric, the b local PCA bases ii∈B\W_i\_i∈ B are obtained from a single call to a batched symmetric eigensolver on the stacked tensor (b,d,d)(b,d,d), rather than b sequential eigendecompositions. 3. Vectorized quadratic/cross terms. The squared columns i∘2W_i 2 are obtained by one elementwise square over the whole batch tensor; the (d2) d2 cross-product columns are obtained by one elementwise multiplication of two fancy-indexed views of the batch tensor (using precomputed index pairs, computed once for all batches), replacing the explicit nested loops over column pairs used in Algorithm 1. 4. Batched matrix products. The second fundamental form approximation IIi=ii⊤I_i=H_iH_i and the local shape operator Si=−IIiiS_i=-I_iC_i are each computed for the entire batch in a single batched matrix product. The b resulting local shape operators are summed and the loop advances to the next batch of b points; S¯ S is obtained by dividing the total sum by n once every batch has been processed. Algorithm 2 summarizes the resulting procedure. Algorithm 2 Vectorized Average Shape Operator Estimation 1: X∈ℝn×dX ^n× d: data matrix; k∈ℕk : number of neighbors; b∈ℕb : batch size 2: S¯∈ℝd×d S ^d× d: average shape operator 3: function Vectorized-Shape-Operator(X,k,bX,k,b) 4: Pii=1n←\P_i\_i=1^n← Batch-Nearest-Neighbors(X,k)(X,k) ⊳ single KD-tree query, O(ndlogn)O(nd n) 5: S¯←∈ℝd×d S 0 ^d× d 6: for each batch B of b indices in 1,…,n\1,…,n\ do 7: ii∈B←\C_i\_i∈ B← Batched-Covariance(Pii∈B)(\P_i\_i∈ B) ⊳ single batched GEMM 8: i,ii∈B←\W_i, _i\_i∈ B← Batched-Eigh(ii∈B)(\C_i\_i∈ B) ⊳ single batched eigensolver call 9: ii∈B←\H_i\_i∈ B← Vectorized-Quadratic-Cross(ii∈B)(\W_i\_i∈ B) ⊳ elementwise + fancy indexing 10: IIii∈B←\I_i\_i∈ B← Batched-Matmul(i,i⊤)(\H_i\,\H_i\ ) 11: Sii∈B←−\S_i\_i∈ B←-Batched-Matmul(IIi,i)(\I_i\,\C_i\) 12: S¯←S¯+∑i∈BSi S← S+ _i∈ BS_i 13: end for 14: return S¯/n S/n 15: end function 2.3.1 Computational Complexity Fix a batch of size b. Gathering and centering the neighbor tensor costs O(bkd)O(bkd). The batched covariance product costs O(bkd2)O(bkd^2) (a (d×k)×(k×d)(d× k)×(k× d) product per sample, over b samples). The batched eigendecomposition costs O(bd3)O(bd^3). Constructing iH_i (the r=d+(d2)=O(d2)r=d+ d2=O(d^2) quadratic and cross-product columns) costs O(bd⋅r)=O(bd3)O(bd· r)=O(bd^3) per batch. The dominant step is IIi=ii⊤I_i=H_iH_i , a (d×r)×(r×d)(d× r)×(r× d) product costing O(d2r)=O(d4)O(d^2r)=O(d^4) per sample, or O(bd4)O(bd^4) per batch; the final product Si=−IIiiS_i=-I_iC_i costs O(bd3)O(bd^3), dominated by the previous term. Summing these per-batch costs and multiplying by the n/bn/b batches gives a total of O(ndlogn+nkd2+nd3+nd4)=O(ndlogn+nkd2+nd4),O\! (nd n+nkd^2+nd^3+nd^4 )=O\! (nd n+nkd^2+nd^4 ), since the O(nd3)O(nd^3) terms are dominated by O(nd4)O(nd^4). This is exactly the complexity of the original per-sample Algorithm 1, O(ndlogn+nd2(k+d2))O(nd n+nd^2(k+d^2)). Remark 1. Batching does not change the asymptotic time complexity of average shape operator estimation: the total number of floating-point operations is unchanged, since Algorithm 2 performs the same per-sample arithmetic as Algorithm 1, only reorganized into ⌈n/b⌉ n/b large batched calls instead of n small sequential ones. The practical speedup obtained is therefore a constant-factor effect, arising from (i) amortizing fixed per-call Python/NumPy dispatch overhead over b samples instead of one, and (i) batched matrix multiplication and eigendecomposition being dispatched to BLAS/LAPACK routines that exploit multi-threading, cache locality, and SIMD vectorization far more effectively than n sequential calls on small matrices. Empirically, on a synthetic benchmark with n=6,430n=6,430 and d=36d=36 (comparable in scale to the satimage dataset in Table 2), we measured a 2.7×2.7× wall-clock speedup from this reformulation alone, with larger relative gains expected as n grows. Remark 2 (Memory–time trade-off). The naive Algorithm 1 holds only one sample’s intermediate matrices in memory at a time, an O(d2)O(d^2) footprint independent of n. Algorithm 2 instead materializes the batch tensor ii∈B\H_i\_i∈ B of shape (b,d,r)(b,d,r) with r=O(d2)r=O(d^2), giving a peak memory footprint of O(bd3)O(bd^3). The batch size b is therefore an explicit, tunable trade-off between wall-clock speed and memory use: larger b amortizes dispatch overhead further but increases peak memory proportionally, which matters in particular for datasets where d is capped at 100100 via the pre-PCA step (Section 3), since O(bd3)O(bd^3) then grows quickly with b. 3 Shape Operator Based PCA Let =1,2,…,n⊂ℝdX=\x_1,x_2,…,x_n\ ^d be a dataset of n observations in a d-dimensional ambient space. We assume that X is sampled from an unknown smooth Riemannian manifold ℳ↪ℝdM ^d of intrinsic dimension m≪dm d. Classical PCA seeks a linear projection P:ℝd→ℝmP:R^d ^m that maximizes the variance of the projected data, operating exclusively on the global covariance matrix Σ=1n∑i=1n(i−)(i−)⊤, = 1n _i=1^n(x_i- μ)(x_i- μ) , (12) where μ is the sample mean. While variance is a natural measure of data spread, it is blind to the intrinsic geometry of ℳM: two datasets with identical covariance structure may lie on manifolds of fundamentally different curvature, and hence exhibit very different cluster structure, class boundaries, and local neighborhoods. To capture this geometric information, we propose to augment the covariance matrix with a curvature-derived term computed from the data manifold. The key geometric object is the shape operator (also known as the Weingarten map), a classical construct from differential geometry that, for each point i∈ℳx_i , encodes the rate and direction at which the manifold curves relative to the ambient space. Intuitively, while the covariance matrix captures where the data spreads, the shape operator captures how the underlying surface bends, and directions of strong curvature are precisely those along which the manifold folds, separates, or forms boundaries between distinct data groups. By aggregating local shape operators into a global mean shape operator ¯ S and using it to regularize Σ , we obtain a geometry-aware covariance Σcurv=Σ+λ¯,λ≥0, _curv= +λ\, S, λ≥ 0, (13) whose principal components reflect both the statistical spread and the intrinsic geometric structure of the data. This regularization requires no class labels, no explicit manifold parametrization, and no nonlinear optimization, making it a natural, interpretable, and computationally tractable extension of classical PCA for unsupervised metric learning. In the limiting cases, the proposed formulation recovers two well-defined methods: when λ=0λ=0, it reduces exactly to standard PCA; and as λ→∞λ→∞, the variance term becomes negligible and the principal components are determined entirely by the mean shape operator ¯ S, yielding a purely curvature-driven spectral embedding that we term Shape Operator Eigenmaps, a novel dimensionality reduction method in its own right, whose theoretical properties and empirical behavior we leave as a subject for future investigation. The optimal balance between these two competing objectives depends directly on the intrinsic geometry of each dataset, and identifying λ∗λ^* allows this trade-off to be resolved in a data-driven manner. A practical difficulty immediately arises, however. The magnitude of ¯ S relative to Σ – and hence the range of λ values that meaningfully interpolate between the variance-dominated and curvature-dominated regimes of Eq. equation 13 – depends on dataset-specific quantities such as the local sampling density, the chosen neighborhood size k, and the manifold’s own curvature scale, none of which are comparable across datasets. Consequently, a single fixed candidate set for λ can correspond to a negligible perturbation of Σ on one dataset and an overwhelming one on another. Selecting λ well is therefore essential, and doing so without recourse to class labels – consistent with the fully unsupervised character of the method – requires a criterion grounded in the geometry of Σcurv(λ) _curv(λ) itself rather than in downstream, label-dependent clustering quality. We describe such a criterion next. 3.1 Unsupervised Estimation of the Regularization Parameter via Spectral Eigengap To select λ without labeled data, we exploit a purely spectral, geometric criterion: the relative separation, or eigengap, between the m-th and (m+1)(m+1)-th eigenvalues of Σcurv(λ) _curv(λ). This closes the one remaining supervised step of an otherwise label-free pipeline – selecting λ via cross-validated grid search against an external clustering index would reintroduce exactly the dependence on labeled data that the rest of the method is designed to avoid. A dataset-comparable candidate grid. Rather than searching a fixed, dataset-independent set of λ values, we center the candidate grid at the scale that equates the trace of the curvature term with that of the covariance, scale=tr(Σ)tr(|¯|),scale= tr( )tr(| S|), (14) and search λj=cj⋅scale _j=c_j·scale for cjc_j ranging over a wide, geometrically spaced, dataset-independent set of dimensionless multipliers (e.g. cj∈[0.1,100]c_j∈[0.1,100]). This ensures that cj=1c_j=1 corresponds to a comparable variance/curvature balance regardless of the dataset at hand, while leaving λ itself unbounded and Eq. equation 13 unchanged; Eq. equation 14 only informs which candidate values of λ are worth evaluating, not the model itself. The eigengap criterion. For each candidate λj _j, let e1(λj)≥⋯≥ed(λj)e_1( _j)≥·s≥ e_d( _j) denote the eigenvalues of Σcurv(λj) _curv( _j), ordered by real part. Because ¯ S is generally not symmetric (Section 2.2), neither is Σcurv(λj) _curv( _j), and its eigenvalues are in general complex; in practice we find them to be real to within numerical precision, and hereafter work with their real parts. We define the relative eigengap at the target dimension m as gapm(λ)=em(λ)−em+1(λ)∑i=1d|ei(λ)|,gap_m(λ)= e_m(λ)-e_m+1(λ) _i=1^d|e_i(λ)|, (15) and select λ∗=argmaxλ∈λ1,…,λJgapm(λ).λ^*= _λ∈\ _1,…, _J\gap_m(λ). (16) We normalize by the total spectral mass ∑i|ei(λ)| _i|e_i(λ)| in Eq. equation 15 rather than by |em(λ)||e_m(λ)|, the more familiar normalization in the spectral clustering literature 22, 29. Because Σcurv(λ) _curv(λ) is asymmetric, its eigenvalues can cross zero as λ varies; at such a crossing, normalizing by |em(λ)||e_m(λ)| produces a spurious, arbitrarily large spike that reflects a vanishing denominator rather than a genuinely well-separated subspace. We verified this failure mode directly in preliminary experiments: the λ selected by the unnormalized criterion at such a spike can yield markedly worse embeddings than neighboring, smoothly varying values of λ. Normalizing by the total spectral mass removes this instability. Theoretical justification for eigengap. The use of the eigengap in Eq. equation 15 as a proxy for the quality of the induced m-dimensional subspace is grounded in classical matrix perturbation theory, in particular the Davis–Kahan sinΘ theorem 6 and its statistically oriented restatement 33. Informally, the Davis–Kahan theorem bounds the angle between an exact m-dimensional invariant subspace of a matrix and the corresponding subspace of a perturbed version of that matrix by a quantity inversely proportional to the gap separating the m-th and (m+1)(m+1)-th eigenvalues of the unperturbed matrix. Applied here, this means that a large gapm(λ)gap_m(λ) certifies that the m-dimensional eigenspace of Σcurv(λ) _curv(λ), and hence the SHOPCA projection itself, is stable: small perturbations of Σcurv(λ) _curv(λ), such as those induced by finite-sample noise in estimating Σ and ¯ S or by resampling the dataset, rotate the resulting subspace only slightly. Conversely, when em(λ)e_m(λ) and em+1(λ)e_m+1(λ) are nearly degenerate, the top-m eigenspace is only weakly identified by the data, and arbitrarily small perturbations can rotate it substantially, yielding an unstable, poorly reproducible projection. Maximizing gapm(λ)gap_m(λ) over λ therefore selects the regularization strength for which the resulting low-dimensional representation is, in this precise and classical sense, most robustly determined by the data, without ever consulting class labels. This is the same stability argument that motivates the eigengap heuristic for choosing the number of retained eigenvectors (equivalently, the embedding dimension) in spectral clustering 22, 29; here it is repurposed to select a regularization strength for a fixed target dimension m, rather than the dimension itself, but the underlying justification is identical. Computational cost. The criterion in Eqs. equation 15–equation 16 is inexpensive relative to the estimation of ¯ S itself. Evaluating gapm(λj)gap_m( _j) for a single candidate requires only forming Σcurv(λj) _curv( _j) and computing its eigenvalues, an O(d3)O(d^3) operation; over a grid of J candidates, the total added cost is O(Jd3)O(Jd^3), negligible next to the O(ndlogn+nd2(k+d2))O(nd n+nd^2(k+d^2)) cost of computing ¯ S (Section 2.2) for any realistic J. Since both Σ and ¯ S are computed without reference to class labels, λ∗λ^* closes the one remaining supervised step of the pipeline discussed in the Introduction, rendering representation learning and hyperparameter selection alike fully unsupervised. 4 Theoretical Justification for Regularization Previous sections introduced Σcurv(λ) _curv(λ) constructively, as a covariance regularized by an estimated curvature term, and Section 3.1 justified the selection of λ via a stability argument grounded in classical perturbation theory. This section complements that construction with four further results that (i) situate SHOPCA within the established framework of regularized covariance estimation, (i) characterize precisely when and why the local and global shape operator estimators admit a real spectrum, (i) give an exact variational characterization of the resulting principal directions, and (iv) quantify, via a second classical perturbation bound, how far the SHOPCA spectrum can move from that of standard PCA as a function of λ. Throughout, ∥⋅∥2\|·\|_2 denotes the spectral (operator) norm. 4.1 SHOPCA as Geometry-Informed Covariance Shrinkage The regularized covariance of Eq. equation 13, Σcurv=Σ+λ¯ _curv= +λ S, is structurally an instance of a well-established family of estimators in multivariate statistics: covariance shrinkage estimators, which combine the sample covariance Σ with a structured target matrix to reduce estimation variance, particularly when n is small relative to d 19. The classical Ledoit–Wolf estimator takes the form Σ^shrink=(1−ρ)Σ+ρF _shrink=(1-ρ) +ρ F, where F is a low-variance, high-bias target – typically a scaled identity, F=τIdF=τ I_d – and ρ∈[0,1]ρ∈[0,1] is chosen to asymptotically minimize the expected quadratic loss ‖Σ^shrink−Σtrue‖F2E\| _shrink- _true\|_F^2 19. SHOPCA instantiates this same bias–variance trade-off, but with a shrinkage target that is neither arbitrary nor generic: ¯ S is estimated from the local curvature of the data manifold, biasing Σcurv _curv toward directions informative about class structure (Section 3) rather than merely toward isotropy. This gives a concrete explanation for the small-sample robustness observed in Section 5 (Table 4): standard shrinkage trades bias for variance with no regard to whether the resulting bias is useful, whereas ¯ S is itself an average of n independent, fixed-size local estimators iS_i, each computed from a k-neighborhood whose size does not shrink with n, so its own estimation variance need not grow as n becomes small. 4.2 Spectral Realness of the Shape Operator Estimator Classical differential geometry guarantees that the continuous shape operator S_x is self-adjoint with respect to the first fundamental form and therefore has a real spectrum 9. We show the discrete estimator iS_i inherits this property by an independent algebraic argument, and use it to explain exactly why the aggregate ¯ S loses the guarantee, corroborating the empirical observation of Section 3.1 that Σcurv(λ) _curv(λ) has, in practice, a numerically real but formally unguaranteed spectrum. Proposition 1. Suppose i≻0C_i 0 (e.g. because k≥d+1k≥ d+1). Then i=−iiS_i=-I_iC_i has real, non-positive eigenvalues. Proof. For square matrices A,BA,B of equal size, ABAB and BABA share the same eigenvalues, with multiplicity 14. Hence iiI_iC_i and iiC_iI_i share the same spectrum. Since i≻0C_i 0, its symmetric square root i1/2C_i^1/2 is invertible, and i−1/2(i)i1/2=i1/2i1/2,C_i^-1/2(C_iI_i)C_i^1/2=C_i^1/2I_iC_i^1/2, so iiC_iI_i is similar to i1/2ii1/2C_i^1/2I_iC_i^1/2, which is symmetric and positive semi-definite: for any x, ⊤(i1/2ii1/2)=(i1/2)⊤i(i1/2)≥0x (C_i^1/2I_iC_i^1/2)x=(C_i^1/2x) I_i(C_i^1/2x)≥ 0, using i=ii⊤⪰0I_i=H_iH_i 0. Similar matrices share eigenvalues, so iiC_iI_i, hence iiI_iC_i, has real, non-negative eigenvalues, and i=−iiS_i=-I_iC_i has real, non-positive eigenvalues. ∎ Proposition 1 is a discrete, purely algebraic echo of the continuous self-adjointness fact, obtained by a different route (similarity to a symmetric PSD matrix rather than self-adjointness with respect to a Riemannian metric). It certifies that each local “principal curvature” estimate κ1(i),…,κd(i)κ^(i)_1,…,κ^(i)_d is always real. Crucially, however, the similarity transform used depends on iC_i, which differs from point to point: iS_i is similar to a symmetric matrix through the point-specific change of basis i1/2C_i^1/2, not through a single transformation shared by all n points. The average ¯=1n∑i S= 1n _iS_i is therefore, in general, not similar to a symmetric matrix through any single transformation, and Proposition 1 does not extend to ¯ S or to Σcurv(λ) _curv(λ). This is precisely why Section 3.1 could only report that the eigenvalues of Σcurv(λ) _curv(λ) are found real numerically, rather than prove they must be: realness is guaranteed at the local level by construction, and is only approximately preserved under aggregation across differing local frames, most closely when the local metrics iC_i vary slowly across the data, i.e. when the manifold’s geometry is locally homogeneous. 4.3 A Variational Characterization of the Principal Directions We now characterize SHOPCA’s directions as stationary points of an explicit variance–curvature objective, and identify precisely where this characterization departs from the operator actually diagonalized in practice. For any square matrix A, write Asym=12(A+A⊤)A^sym= 12(A+A ), Askew=12(A−A⊤)A^skew= 12(A-A ). For any v, ⊤Askew=0v A^skewv=0 (a scalar equal to its own transpose, −⊤Askew-v A^skewv). Hence ⊤Σcurv(λ)=⊤Σ+λ⊤¯sym=⊤Σcurvsym(λ),v _curv(λ)\,v=v \,v+λ\,v S^symv=v _curv^sym(λ)\,v, (17) where Σcurvsym(λ):=Σ+λ¯sym _curv^sym(λ):= +λ S^sym is symmetric. Proposition 2. The stationary points of the Rayleigh quotient R()=⊤Σcurv(λ)/⊤R(v)=v _curv(λ)v/v v coincide exactly with the eigenvectors of Σcurvsym(λ) _curv^sym(λ), with critical values equal to its (real) eigenvalues. Proof. Immediate from Eq. equation 17 and the Rayleigh–Ritz theorem applied to the symmetric matrix Σcurvsym(λ) _curv^sym(λ) 14. ∎ Proposition 2 gives SHOPCA a clean interpretive story: its principal directions extremize a convex combination of variance and (symmetrized) curvature. Its scope must be stated precisely, however: the projection of Section 3 diagonalizes Σcurv(λ) _curv(λ) itself via a general, non-symmetric eigendecomposition, not Σcurvsym(λ) _curv^sym(λ), and the eigenvectors of a non-normal matrix need not coincide with, or be orthogonal to, those of its symmetric part. The two sets of principal directions coincide exactly only when ¯skew=0 S^skew=0, and are close whenever ‖¯skew‖2≪‖¯sym‖2\| S^skew\|_2 \| S^sym\|_2. This gives a concrete, checkable diagnostic for how faithfully Proposition 2’s variational reading describes the directions the algorithm actually returns, and points to a natural symmetrized variant, Σcurvsym(λ) _curv^sym(λ), whose eigenvectors are guaranteed real and orthogonal by the spectral theorem, as a principled ablation for future work. 4.4 Eigenvalue Stability under Regularization Section 3.1 argued, via Davis–Kahan, that a large eigengap certifies eigenspace stability under perturbation. We complement this with a direct, quantitative bound on how far the eigenvalues themselves can move from those of standard PCA as λ grows, via the Bauer–Fike theorem 1, which bounds eigenvalue perturbation for a diagonalizable matrix under an arbitrary, not necessarily symmetric, additive perturbation. Proposition 3. Let e1(λ)≥⋯≥ed(λ)e_1(λ)≥·s≥ e_d(λ) be the eigenvalues of Σcurv(λ) _curv(λ), ordered by real part, and e1,…,ede_1,…,e_d the eigenvalues of Σ . For every λ≥0λ≥ 0 and j, there exists i such that |ej(λ)−ei|≤λ‖¯‖2.|e_j(λ)-e_i|\;≤\;λ\,\| S\|_2. (18) Proof. Σ is symmetric, hence orthogonally diagonalizable, Σ=Vdiag(e1,…,ed)V⊤ =Vdiag(e_1,…,e_d)V with V orthogonal, so κ2(V)=1 _2(V)=1. Applying Bauer–Fike 1 to Σ with perturbation E=λ¯E=λ S gives, for every eigenvalue ej(λ)e_j(λ) of Σ+E +E, mini|ej(λ)−ei|≤κ2(V)‖E‖2=λ‖¯‖2 _i|e_j(λ)-e_i|≤ _2(V)\|E\|_2=λ\| S\|_2. ∎ Eq. equation 18 makes precise that Σcurv(λ) _curv(λ) interpolates continuously away from the PCA spectrum as λ grows, at a rate governed entirely by ‖¯‖2\| S\|_2. The bound guarantees only the existence of a nearby unperturbed eigenvalue, not that the correspondence preserves the original ordering i=ji=j, consistent with, and giving formal underpinning to, the empirical observation of Section 3.1 that the eigenvalue ordering of Σcurv(λ) _curv(λ) can change discontinuously with λ, which motivated normalizing the eigengap criterion by total spectral mass rather than by |em(λ)||e_m(λ)| alone. 4.5 SHOPCA as a Regularized Mahalanobis Metric Finally, we make explicit the sense in which Σcurv(λ) _curv(λ) realizes the co-design principle motivated in the introduction: that projection and metric should be learned jointly 13, 30. Whenever invertible, Σcurv(λ) _curv(λ) induces a Mahalanobis-type dissimilarity dcurv2(,,λ)=(−)⊤Σcurv(λ)−1(−),d_curv^2(x,y;λ)=(x-y) _curv(λ)^-1(x-y), (19) generalizing the classical PCA-induced metric d2(,)=(−)⊤Σ−1(−)d^2(x,y)=(x-y) ^-1(x-y), recovered exactly at λ=0λ=0. This formalizes SHOPCA as unsupervised metric learning in the precise sense used elsewhere in the paper: the projection defined by the leading eigenvectors of Σcurv(λ) _curv(λ) and the metric that same matrix induces via Eq. equation 19 are, by construction, one and the same object: the co-design property identifies as the guiding principle, achieved here without the labeled supervision required by joint frameworks such as 13. 5 Computational Experiments and Results This section presents a comprehensive empirical evaluation of SHOPCA against three reference methods, organized into three complementary experiments, each designed to isolate a different aspect of the proposed method’s behavior. In the first experiment, SHOPCA is compared directly against standard PCA (Eq. equation 12), the classical linear baseline it generalizes, in order to quantify the gain attributable to the curvature-aware regularization of Eq. equation 13 in isolation from any other methodological difference. In the second experiment, SHOPCA is compared against Isomap 27, one of the most representative and widely adopted nonlinear manifold learning algorithms, which recovers a low-dimensional embedding by applying classical multidimensional scaling to pairwise geodesic distances estimated over a neighborhood graph. Isomap was selected as a benchmark precisely because it embodies a methodological paradigm markedly different from both PCA and SHOPCA: whereas PCA is linear and variance-only, and SHOPCA remains linear while incorporating local curvature, Isomap is fully nonlinear and reconstructs global manifold structure from local geodesic distance estimates, providing a natural point of comparison for assessing how much of Isomap’s nonlinear representational power can be matched by a closed-form linear method enriched with local differential-geometric information. In the third experiment, we compare SHOPCA against UMAP 21, a state-of-the-art nonlinear embedding method, restricting attention to datasets with small sample sizes (S); this restriction is motivated by a well-documented limitation of neighborhood-graph-based methods such as UMAP and Isomap, whose local graph estimation is known to degrade when few samples are available, a regime in which SHOPCA’s covariance-based curvature estimation is, by construction, less reliant on graph connectivity and therefore expected to remain more robust. The experiments were conducted on more than 50 publicly available, real-world benchmark datasets retrieved from the OpenML repository 28 (openml.org), spanning a wide range of sample sizes, feature dimensionalities, and numbers of classes, and covering diverse application domains including image recognition, spectroscopy, genomics, and tabular data. This heterogeneity provides a rigorous testbed for evaluating the generalization of SHOPCA across markedly different geometric regimes, rather than restricting the evaluation to a narrow class of well-behaved datasets. For every method and every dataset, the same evaluation protocol is applied: each method projects the data onto a two-dimensional latent space, and Agglomerative Clustering with Ward’s linkage 31 is then applied in this latent space to recover c clusters, where c is set to the ground-truth number of classes for each dataset. The resulting partitions are scored against the ground-truth labels using four external clustering validity indices: the Adjusted Rand Index (ARI) 15, Normalized Mutual Information (NMI) 26, the Fowlkes–Mallows index (FM) 11, and the V-measure 25. Ground-truth labels are used exclusively for this final evaluation step, never within the projection methods themselves, preserving the fully unsupervised character of SHOPCA, PCA, ISOMAP, and UMAP alike. Because the local shape operator estimation described in Section 2.2 scales as O(nd2(k+d2))O(nd^2(k+d^2)) in the ambient dimensionality d, it becomes computationally impractical on the higher-dimensional datasets in our benchmark suite without an intermediate dimensionality reduction step. We therefore apply standard PCA as a pre-processing stage whenever the number of features exceeds a threshold T=25T=25, projecting the data onto its first T principal components before estimating the local shape operators and proceeding with the remainder of the pipeline. This threshold is applied uniformly to SHOPCA and, where applicable, to the reference methods, ensuring that all compared methods operate on the same effective input representation and that observed performance differences reflect the projection methods themselves rather than differences in pre-processing. In all experiments, the parameter k (number of neighbors in the k-N graph) is set to log2n _2n. Table 1 lists the more than 50 real-world benchmark datasets used throughout the computational experiments, reporting for each the dataset name, the number of samples (n), the number of features (d), and the number of classes (c). The selected datasets span several orders of magnitude in n, d, and c, and cover a broad range of application domains, providing a diverse and heterogeneous testbed for evaluating SHOPCA under markedly different geometric regimes. The full set of datasets in Table 1 is used in the first two experiments (SHOPCA vs. PCA and SHOPCA vs. Isomap); the third experiment, comparing SHOPCA against UMAP, is restricted to the subset of datasets with small sample sizes. Table 1: Datasets used in the computational experiments: each row shows the dataset name, number of samples (n), number of features (d), and number of classes (c). # Dataset n d c # Dataset n d c 1 iris 150 4 3 30 har 10299 561 6 2 wine 178 13 3 31 Olivetti_Faces 400 4096 40 3 digits 1797 64 10 32 balance-scale 625 4 3 4 diggle_table_a2 310 8 9 33 mammography 11183 6 2 5 satimage 6430 36 6 34 abalone 4177 8 3 6 mfeat-morphological 2000 6 10 35 Engine1 383 5 3 7 mfeat-zernike 2000 47 10 36 penguins 344 6 3 8 mfeat-karhunen 2000 64 10 37 solar-flare 1066 12 6 9 mfeat-fourier 2000 77 10 38 usp05 203 16 11 10 mfeat-factors 2000 216 10 39 Breast 699 10 2 11 mfeat-pixel 2000 240 10 40 prnn_synth 250 2 2 12 isolet 7797 617 26 41 irish 500 5 2 13 optdigits 5620 64 10 42 Satellite 5100 36 2 14 pendigits 10992 16 10 43 zoo 101 16 7 15 semeion 1593 256 10 44 glass 214 9 6 16 mnist_784 (20%) 14000 784 10 45 prnn_viruses 61 18 4 17 Fashion-MNIST (20%) 14000 784 10 46 ionosphere 351 34 2 18 Kuzushiji-MNIST (20%) 14000 784 10 47 ThreeOf9 512 9 2 19 cardiotocography 2126 35 10 48 eucalyptus 736 19 5 20 cnae-9 1080 856 9 49 sonar 208 60 2 21 JapaneseVowels 9961 14 9 50 blood-transfusion-service-center 748 4 2 22 segment 2310 19 7 51 lymph 148 18 4 23 Indian_pines 9144 220 8 52 teachingAssistant 151 6 3 24 letter 20000 16 26 53 backache 180 31 2 25 wap.wc 1560 8460 20 54 heart-h 294 13 2 26 texture 5500 40 11 55 ecoli 336 7 8 27 USPS 9298 256 10 56 liver-disorders 342 5 16 28 coil-20 1440 1024 20 57 LED-display-domain-7digit 500 7 10 29 UMIST_Faces_Cropped 575 10304 20 5.1 First Experiment: Regular PCA versus Shape Operator PCA The first experiment establishes a direct baseline by comparing standard PCA against the proposed SHOPCA across the full benchmark suite of 30 real-world datasets. Our objective is to isolate the marginal gain attributable to curvature-aware covariance regularization: both methods project the data onto a two-dimensional latent space, after which agglomerative clustering with Ward’s linkage is applied to recover the ground-truth number of clusters, and the resulting partitions are evaluated using the Adjusted Rand Index (ARI), Fowlkes-Mallows index (FM), and V-measure. Because SHOPCA reduces exactly to PCA when λ=0λ=0, any improvement can be attributed solely to the mean shape operator term (geometric regularization). The regularization parameter λ is estimated using the eigengap method, which is fully unsupervised (label free). If the intrinsic geometry of the data manifold carries discriminative structure beyond what variance alone captures, SHOPCA should systematically yield more compact and well-separated clusters, a hypothesis that is supported by the obtained results. Table 2: Comparison between standard PCA and the proposed SHOPCA on 3030 real-world benchmark datasets. Both methods project the data onto a two-dimensional latent space, after which agglomerative clustering with Ward’s linkage is applied to recover the ground-truth number of classes. Clustering quality is assessed by the Adjusted Rand Index (ARI), Fowlkes–Mallows index (FM), and V-measure (VM). Bold entries indicate the best result per metric. SHOPCA achieves superior scores, with average relative improvements of 212%212\% (ARI), 73%73\% (FM), and 140%140\% (VM) over regular PCA for these datasets. Regular PCA Shape Operator PCA Datasets ARI FM VM ARI FM VM λ iris 0.6254 0.7542 0.6722 0.7312 0.8222 0.7701 55.19 wine 0.6583 0.7734 0.6986 0.6764 0.7926 0.7252 12.46 digits 0.0358 0.2148 0.1322 0.3372 0.4228 0.5222 17.51 diggle_table_a2 0.2896 0.3805 0.5704 0.4708 0.5638 0.7078 127.60 satimage 0.0057 0.2840 0.0092 0.1012 0.4603 0.2627 9,95 mefat-zernike 0.0183 0.1744 0.0619 0.2884 0.3779 0.4718 17.03 mfeat-karhunen 0.0415 0.1606 0.1034 0.3118 0.3993 0.5110 16.70 mfeat-fourier 0.0354 0.1666 0.0916 0.2720 0.3643 0.4760 18.34 mfeat-factors 0.0249 0.1578 0.0737 0.4208 0.4855 0.6040 18.79 mfeat-pixel 0.1494 0.2455 0.2936 0.3483 0.4398 0.5474 17.23 isolet 0.0227 0.0821 0.1167 0.1953 0.2495 0.5085 20.35 optdigits 0.0499 0.1708 0.1007 0.2464 0.3607 0.4670 20.61 pendigits 0.3234 0.4044 0.4989 0.4004 0.4709 0.5622 127.48 semeion 0.0218 0.1644 0.0985 0.2213 0.3182 0.3956 8.45 mnist_784 (20%) 0.0299 0.1727 0.0984 0.1800 0.3035 0.3268 12.96 Fashion-MNIST (20%) 0.0202 0.2050 0.0855 0.3050 0.3917 0.4887 32.86 Kuzushiji-MNIST (20%) 0.0575 0.1992 0.1474 0.1163 0.2329 0.2393 12.55 cardiotocography 0.0422 0.2218 0.0885 0.6725 0.7340 0.6938 19.64 cnae-9 0.0128 0.2297 0.0968 0.1411 0.3377 0.4122 18.49 JapaneseVowels 0.0504 0.1945 0.1875 0.1276 0.2390 0.2702 33.04 segment 0.1529 0.4197 0.3302 0.2892 0.4603 0.5284 345.71 Indian_pines 0.0061 0.2726 0.0124 0.3821 0.5835 0.5802 23.68 letter 0.0345 0.1050 0.1939 0.0905 0.1405 0.3208 44.86 wap.wc 0.0285 0.2413 0.2683 0.0448 0.2343 0.3286 21.29 texture 0.0192 0.2002 0.0648 0.5444 0.6055 0.7412 36.64 USPS 0.1123 0.2345 0.2180 0.3621 0.4360 0.4754 22.82 coil-20 0.1072 0.1888 0.2910 0.4357 0.4709 0.6643 64.64 UMIST_Faces_Cropped 0.0857 0.1499 0.3535 0.4371 0.4733 0.7020 41.50 har 0.0410 0.2600 0.0675 0.5022 0.6115 0.6465 30.49 Olivetti_Faces 0.0302 0.0588 0.4326 0.1236 0.1575 0.5566 8.35 Average 0.1044 0.2496 0.2153 0.3259 0.4313 0.5169 Median 0.0384 0.2026 0.1244 0.3084 0.4294 0.5166 Table 2 presents the complete results of the first experiment, comparing standard PCA against SHOPCA on 30 diverse real-world datasets. The evidence is unequivocal: SHOPCA achieves strictly superior clustering performance on every single dataset across all three external validity indices (ARI, FM, and V-measure), with no exceptions. This 30/30 sweep is not a marginal effect, it represents a systematic dominance that strongly validates the central hypothesis of this work: local curvature carries discriminative geometric information that variance alone cannot capture. The aggregate statistics reveal the magnitude of the improvement. Averaged across the benchmark suite, SHOPCA attains ARI, FM, and V-measure scores of 0.3259 , 0.4313 , and 0.5169 , compared to 0.1044 , 0.2496 , and 0.2153 for standard PCA, corresponding to relative gains of 212% , 73% , and 140% , respectively. Crucially, the median scores exhibit the same pattern (0.3084 vs. 0.0384 in ARI; 0.4294 vs. 0.2026 in FM; 0.5166 vs. 0.1244 in VM), confirming that the gains are not driven by a handful of outlier datasets but reflect a consistent tendency across heterogeneous data regimes. Perhaps the most compelling evidence comes from datasets where standard PCA performs near chance level. On satimage, PCA yields an ARI of 0.0057 and a V-measure of 0.0092, effectively failing to recover any cluster structure, whereas SHOPCA raises these scores to 0.1012 and 0.2627 , representing an 18-fold and 29-fold improvement, respectively. Similarly, on Indian_pines, PCA achieves an ARI of 0.0061, while SHOPCA reaches 0.3821, a relative gain exceeding 60×60×. These dramatic rescues suggest that in datasets with pronounced manifold curvature, variance-based projections collapse distinct classes onto overlapping subspaces, whereas curvature-aware regularization successfully disentangles them by exploiting the second-order geometry of the embedding. The gains are equally pronounced on structured visual and spectral datasets. For the mfeat family, where PCA struggles with ARI scores below 0.05, SHOPCA consistently lifts performance into the 0.27-0.42 range. On image datasets such as coil-20, USPS, and UMIST_Faces_Cropped, SHOPCA improves ARI by factors of 4×4×, 3.2×3.2× and 5.1×5.1×, respectively. This pattern corroborates the theoretical intuition that natural image manifolds, characterized by non-linear illumination, pose, and texture variations, exhibit rich local curvature that PCA, by construction, ignores. By steering the principal components toward directions of both maximum variance and maximum curvature, SHOPCA preserves the geometric folds that separate semantic classes. Notably, even on datasets where PCA already performs reasonably well, SHOPCA provides consistent, non-negligible improvements. On iris and wine, where PCA achieves ARI scores of 0.6254 and 0.6583, SHOPCA pushes these to 0.7312 and 0.6764. The relative gains here are more modest (17% and 3% in ARI), which is theoretically expected: when the underlying manifold is approximately flat or the class structure aligns well with the directions of maximum variance, curvature contributes less additional discriminative signal. The fact that SHOPCA doesn’t degrade performance, even in these favorable scenarios, is an important practical reassurance: curvature regularization is a safe, beneficial augmentation that adapts its contribution to the intrinsic geometry of the data. The selected regularization parameters λ∗λ^* span more than two orders of magnitude, from 8.35 on Olivetti_Faces to 345.71 on segment, with no dataset selecting λ∗=0λ^*=0. This variability is diagnostically meaningful. It confirms that the optimal balance between variance and curvature is intrinsically data-dependent, and that a purely variance-driven projection (i.e., standard PCA) is never optimal within the proposed framework. Datasets with complex, highly curved manifolds such as segment and diggle_table_a2 demand strong curvature regularization (λ∗>100λ^*>100), whereas simpler geometries require only moderate weighting. The absence of any dataset preferring λ∗=0λ^*=0 constitutes strong empirical evidence that curvature information is universally present and useful, even if its relative importance varies. Taken together, these results establish that the mean shape operator is not merely a theoretical curiosity from differential geometry, but a powerful, practical regularizer for unsupervised representation learning. By making curvature computationally tractable within a linear, closed-form framework, SHOPCA bridges the gap between the geometric richness of non-linear manifold learning and the simplicity and interpretability of PCA, delivering superior clustering performance across 30 independent validations without a single failure case. To complement the quantitative analysis, Figure 2 provides a qualitative comparison of the two-dimensional embeddings produced by standard PCA and the proposed SHOPCA on representative datasets: cardiotocography, mfeat-feactors and texture. The visual evidence corroborates the statistical findings: standard PCA projections exhibit substantial inter-class overlap, with distinct clusters collapsing onto shared regions of the latent space due to the exclusive prioritization of variance-maximizing directions. In contrast, SHOPCA yields markedly more separated and compact cluster structures, as the curvature-aware regularization steers the principal components toward directions where the manifold folds and bends, thereby aligning the embedding with the intrinsic geometric boundaries between classes. Notably, these results are obtained with the regularization parameter λ estimated entirely via the eigengap heuristic, a fully unsupervised, geometry-driven criterion that requires no labeled validation data or manual tuning, preserving the unsupervised integrity of the pipeline while delivering discriminative representations. Figure 2: Qualitative comparison of the two-dimensional embeddings produced by standard PCA and the proposed SHOPCA on representative datasets: cardiotocography (first row), mfeat-factors (second row) and texture (third row). 5.2 Second Experiment: ISOMAP versus Shape Operator PCA While the preceding experiment establishes that curvature-aware regularization improves upon its linear baseline, a natural question is whether SHOPCA remains competitive against established nonlinear manifold learning methods that explicitly model global geometric structure. ISOMAP 27 provides an ideal comparator: it recovers the intrinsic geometry of a data manifold by approximating geodesic distances along the manifold via shortest paths through a neighborhood graph, and subsequently embeds the data via multidimensional scaling. This two-stage pipeline, local neighborhood graph construction followed by global distance preservation, represents the classical geometric approach to nonlinear dimensionality reduction, one that fundamentally differs from SHOPCA’s strategy of encoding local curvature directly into a linear covariance regularizer. By comparing against ISOMAP, we can assess whether the computational simplicity and closed-form nature of SHOPCA come at the cost of geometric fidelity, or whether local curvature information, when aggregated globally, captures manifold structure as effectively as geodesic distance preservation. We evaluate both methods on 2525 datasets under identical clustering protocols, projecting each dataset to two dimensions and measuring the quality of the resulting partitions through ARI, FM, and V-measure. Table 3: Comparison between ISOMAP and the proposed SHOPCA on 2525 real-world benchmark datasets. Both methods project the data onto a two-dimensional latent space, after which agglomerative clustering with Ward’s linkage is applied to recover the ground-truth number of classes. Clustering quality is assessed by the Adjusted Rand Index (ARI), Fowlkes–Mallows index (FM), and V-measure (VM). Bold entries indicate the best result per metric. ISOMAP Shape Operator PCA Datasets ARI FM VM ARI FM VM lambda iris 0.7074 0.8084 0.7772 0.7312 0.8222 0.7701 55.19 digits 0.3288 0.4536 0.5855 0.3372 0.4228 0.5222 17.51 diggle_table_a2 0.4609 0.5450 0.7168 0.4708 0.5638 0.7078 127.60 mfeat-morphological 0.2873 0.3979 0.5048 0.4381 0.5242 0.6280 256.75 mfeat-fourier 0.2771 0.3610 0.4364 0.2720 0.3643 0.4760 18.34 Fashion-MNIST (20%) 0.2859 0.3858 0.4897 0.3050 0.3917 0.4887 32.86 cardiotocography 0.4951 0.5736 0.5864 0.6725 0.7340 0.6938 19.64 cnae-9 0.1379 0.2645 0.2993 0.1058 0.3161 0.3614 18.52 segment 0.2981 0.4634 0.4874 0.2892 0.4603 0.5284 345.71 texture 0.4236 0.4952 0.6409 0.5444 0.6055 0.7412 36.64 USPS 0.3112 0.4210 0.4761 0.3621 0.4360 0.4754 22.82 coil-20 0.3172 0.3681 0.6457 0.3685 0.4137 0.6413 64.67 UMIST_Faces_Cropped 0.2567 0.3148 0.5853 0.4023 0.4438 0.6989 41.55 har 0.4895 0.5903 0.6231 0.5022 0.6115 0.6465 30.49 balance_scale 0.0358 0.4188 0.0479 0.2415 0.5376 0.2522 9.50 mammography 0.0682 0.9703 0.0176 0.3127 0.9807 0.2209 79.50 abalone 0.0347 0.4121 0.0297 0.0749 0.4796 0.1679 173.49 Engine1 -0.0658 0.5916 0.1209 0.1648 0.5692 0.2699 16.44 penguins 0.2957 0.6021 0.3107 0.5030 0.7138 0.6517 60.64 solar-flare 0.0546 0.4072 0.1021 0.2664 0.4812 0.4689 1,360.02 usp05 0.1290 0.3348 0.3244 0.1757 0.4048 0.3502 20.31 Breast -0.0027 0.7370 0.0037 0.5717 0.8204 0.4655 9.03 prnn_synth 0.1275 0.6459 0.2485 0.4178 0.7199 0.3970 53.30 irish -0.0004 0.5019 0.0012 0.1439 0.6156 0.2233 279.84 Satellite 0.0316 0.7651 0.0262 0.1659 0.9801 0.0661 9.41 Average 0.2314 0.5132 0.3635 0.3536 0.5765 0.4765 Median 0.2771 0.4634 0.4364 0.3372 0.5376 0.4760 Table 3 reports the comparative performance of ISOMAP and SHOPCA on 25 benchmark datasets, designed to test whether a linear, curvature-aware regularizer can match or exceed a classical nonlinear geodesic embedding. The aggregate picture strongly favors SHOPCA: across all three external validity indices, SHOPCA achieves higher mean and median scores (ARI: 0.3536 vs. 0.2314 ; FM: 0.5765 vs. 0.5132 ; VM: 0.4765 vs. 0.3635 ), with relative improvements of 53% , 12% , and 31% , respectively. More precisely, SHOPCA secures the best score on all three metrics simultaneously in 16 of the 25 datasets, whereas ISOMAP never achieves a clean sweep. This 16/25 dominance is particularly significant because ISOMAP operates with strictly greater modeling capacity: it is nonlinear, graph-based, and explicitly preserves global geodesic structure, yet SHOPCA’s linear, curvature-regularized projection extracts more discriminative geometry for the downstream clustering task. The nine remaining datasets reveal a nuanced but still favorable pattern for SHOPCA. In seven of them (iris, diggle_table_a2, mfeat-fourier, Fashion-MNIST, cnae-9, USPS, coil-20), the methods split the metrics: ISOMAP occasionally wins on V-measure or ARI in isolation, but SHOPCA dominates the other two. For instance, on digits, ISOMAP achieves superior FM and V-measure (0.4536 and 0.5855), yet its ARI (0.3288) trails SHOPCA’s (0.3372); conversely, on segment, ISOMAP leads in ARI and FM, but SHOPCA claims V-measure (0.5284 vs. 0.4874). These splits are diagnostically informative. ISOMAP’s occasional V-measure advantages suggest that its geodesic embedding can produce some pure, compact clusters when the neighborhood graph faithfully approximates the manifold; however, its simultaneous losses in ARI or FM indicate that this purity comes at the cost of globally inconsistent partitions, misallocating classes in a way that penalizes the pairwise agreement metrics more severely. SHOPCA, by contrast, tends to distribute its gains more uniformly across indices, implying that curvature regularization produces a geometrically coherent subspace rather than a patchwork of locally faithful but globally incoherent neighborhoods. The most striking results occur in datasets where ISOMAP’s graph-based manifold estimation collapses entirely. On Engine1, ISOMAP yields a negative ARI (−0.0658 ), indicating performance worse than random assignment, while SHOPCA rescues the clustering to a meaningful 0.1648 . Similarly catastrophic failures for ISOMAP appear on Breast (ARI −0.0027), irish (ARI −0.0004), and prnn_synth (ARI 0.1275), all of which SHOPCA elevates to positive, substantial scores (0.5717, 0.1439, and 0.4178, respectively). These rescues are not marginal improvements; they represent qualitative phase transitions from non-informative to highly discriminative representations. The underlying cause is well understood: ISOMAP’s reliance on a k-nearest neighbor graph becomes brittle when local sampling is sparse or when the manifold curvature varies sharply, producing shortcuts or disconnected components that corrupt geodesic distance estimates. SHOPCA bypasses this fragility entirely by deriving curvature from local covariance eigendecompositions, a statistical operation that remains stable even when the neighborhood graph would be unreliable. Even in regimes where ISOMAP performs credibly, SHOPCA closes the gap or surpasses it without the computational overhead of shortest-path computations and iterative eigendecompositions of dense distance matrices. On the image datasets coil-20 and UMIST_Faces_Cropped, for example, SHOPCA improves ARI by 16% and 57%, respectively, despite ISOMAP’s theoretical advantage in modeling the non-linear pose and illumination manifolds inherent to visual data. The texture dataset provides the most decisive SHOPCA victory: ARI jumps from 0.4236 to 0.5444, and V-measure from 0.6409 to 0.7412, suggesting that the high-dimensional textural manifold is characterized by curvature heterogeneity that ISOMAP’s uniform geodesic metric fails to resolve, but which SHOPCA captures through its local second-order regularization. Overall, these results challenge the prevailing assumption that nonlinear geodesic embeddings are necessarily superior to linear methods for manifold-structured data. By encoding local curvature directly into the covariance structure rather than approximating it indirectly through graph-based geodesics, SHOPCA achieves comparable clustering fidelity at a fraction of the computational cost, while remaining robust in precisely the small-sample and high-curvature regimes where ISOMAP is known to degrade. To complement the quantitative analysis, Figure 3 presents a qualitative comparison of the two-dimensional embeddings produced by ISOMAP and the proposed SHOPCA on three representative datasets: cardiotocography, UMIST_Faces_Cropped, and texture. While ISOMAP successfully unfolds the global manifold geometry, it does so at the cost of local cluster morphology: the geodesic distance preservation criterion tends to stretch intrinsically curved neighborhoods into elongated, filamentary structures, producing clusters that resemble extended chains rather than compact, isotropic point clouds. This geometric distortion is particularly detrimental to classical clustering algorithms, such as k-means, agglomerative clustering with Ward’s linkage, and Gaussian mixture models, which implicitly assume that clusters occupy convex, densely packed regions of the feature space. In contrast, SHOPCA yields markedly more compact and well-delineated groupings. By regularizing the covariance matrix with local curvature rather than globally stretching the embedding to preserve geodesics, SHOPCA respects the natural density boundaries of each class, producing spherical, tightly concentrated clusters that align with the inductive biases of standard centroid-based and density-based partitioners. The visual evidence therefore corroborates the statistical findings: SHOPCA not only improves quantitative clustering scores, but fundamentally transforms the geometry of the latent space into a representation that is structurally more amenable to unsupervised class discovery. Figure 3: Qualitative comparison of the two-dimensional embeddings produced by standard PCA and the proposed SHOPCA on representative datasets: cardiotocography (first row), UMIST_Faces_Cropped (second row) and texture (third row). 5.3 Third Experiment: UMAP versus Shape Operator PCA UMAP has emerged as the de facto standard for nonlinear dimensionality reduction in unsupervised representation learning, prized for its scalability, visual fidelity, and preservation of both local and global manifold structure 21. Yet its success hinges on a critical assumption: that the data density is sufficient to support a reliable k-nearest neighbor graph. In small-sample settings, the neighborhood graph becomes sparse, disconnected, and geometrically inconsistent, causing UMAP’s iterative cross-entropy optimization to collapse into degenerate embeddings that poorly reflect the true manifold topology 5. This vulnerability is not merely a hyperparameter nuisance but a fundamental limitation of graph-based, iterative manifold learners. By contrast, SHOPCA requires no neighborhood graph construction and no iterative optimization: curvature information is injected directly into the global covariance matrix through the mean shape operator, a closed-form statistical aggregate of local eigendecompositions. Consequently, SHOPCA’s geometric regularization remains stable even when the number of observations is too small to support reliable graph inference. In this experiment, we evaluate both methods on 2828 small-sample datasets, restricting the comparison to regimes where UMAP’s graph-based assumptions are most stressed, to determine whether SHOPCA can serve as a principled, non-iterative alternative when data is scarce. Table 4: Comparison between UMAP and the proposed SHOPCA on 2828 real-world benchmark datasets. Both methods project the data onto a two-dimensional latent space, after which agglomerative clustering with Ward’s linkage is applied to recover the ground-truth number of classes. Clustering quality is assessed by the Adjusted Rand Index (ARI), Fowlkes–Mallows index (FM), and V-measure (VM). Bold entries indicate the best result per metric. UMAP Shape Operator PCA Datasets ARI FM VM ARI FM VM lambda iris 0.6312 0.7606 0.7039 0.7312 0.8222 0.7701 55.19 diggle_table_a2 0.4073 0.4802 0.6637 0.4708 0.5638 0.7078 127.60 balance_scale 0.0942 0.4419 0.0811 0.2415 0.5376 0.2522 9.50 zoo 0.5203 0.6225 0.7311 0.8516 0.8895 0.8224 697.51 glass 0.1496 0.3383 0.2652 0.1589 0.4203 0.3024 21.38 mammography -0.0297 0.8061 0.0131 0.3127 0.9807 0.2209 79.50 abalone 0.1332 0.4229 0.1253 0.0749 0.4796 0.1679 173.49 prnn_viruses 0.2507 0.5143 0.4849 0.3610 0.6206 0.4335 66.40 Engine1 0.1311 0.5428 0.2637 0.1648 0.5692 0.2699 16.44 penguins 0.0921 0.4282 0.1906 0.5030 0.7138 0.6517 60.64 UMIST_Faces_Cropped 0.3295 0.3657 0.6383 0.4611 0.4944 0.7152 41.47 ionosphere 0.0044 0.5381 0.0319 0.0374 0.7293 0.0626 40.88 ThreeOf9 0.0069 0.5037 0.0064 0.0911 0.5458 0.0685 0.32 solar-flare 0.0536 0.2655 0.1547 0.2664 0.4812 0.4689 1,360.02 eucalyptus 0.0415 0.2629 0.0699 0.0645 0.2961 0.0714 1,187.58 sonar -0.0040 0.5060 0.0002 0.0488 0.6096 0.0510 9.09 blood-transfusion-service-center -0.0521 0.6505 0.0152 0.0311 0.7968 0.0208 36.33 usp05 0.1315 0.3241 0.3667 0.1757 0.4048 0.3502 20.31 Breast -0.0378 0.6038 0.1439 0.5717 0.8204 0.4655 9.03 lymph 0.0852 0.4272 0.1029 0.0605 0.5349 0.1971 6.95 teachingAssistant 0.0405 0.4022 0.0620 0.0477 0.4542 0.0834 40.59 backache 0.0396 0.6786 0.0082 0.1041 0.7350 0.0282 13.66 prnn_synth 0.1512 0.6308 0.1945 0.4178 0.7199 0.3970 53.30 heart-h 0.0327 0.5443 0.0588 0.3365 0.6924 0.2361 11.01 ecoli 0.3925 0.5323 0.5713 0.5495 0.6747 0.4660 1,615.31 liver-disorders 0.0025 0.1107 0.1317 0.0239 0.1630 0.1438 10.62 LED-display-domain-7digit 0.2470 0.3305 0.3974 0.2469 0.3545 0.4145 76.11 irish -0.0004 0.5019 0.0012 0.1439 0.6156 0.2233 279.84 Average 0.1373 0.4835 0.2313 0.2696 0.5971 0.3237 Median 0.0887 0.5028 0.1378 0.2086 0.5894 0.2610 Table 4 presents the results of the third experiment, comparing SHOPCA against UMAP on 28 small-sample datasets where neighborhood-graph methods are known to be most vulnerable. The aggregate statistics reveal a decisive advantage for SHOPCA: mean ARI improves from 0.1373 to 0.2696 (96% relative gain), mean FM from 0.4835 to 0.5971 (24%), and mean V-measure from 0.2313 to 0.3237 (40%), with median scores exhibiting the same directional dominance. SHOPCA achieves the best result on all three metrics simultaneously in 23 of the 28 datasets (82%), a remarkable margin given that UMAP is widely regarded as the state of the art in nonlinear unsupervised representation learning. The five remaining datasets expose the characteristic failure modes of UMAP in data-scarce regimes. On abalone, lymph, and LED-display-domain-7digit, UMAP secures a marginally higher ARI (0.1332 vs. 0.0749, 0.0852 vs. 0.0605, and 0.2470 vs. 0.2469, respectively), yet loses simultaneously on FM and V-measure in all three cases. This split is diagnostically significant: ARI measures pairwise agreement adjusted for chance, whereas FM penalizes the joint precision and recall of the clustering against ground-truth class boundaries. UMAP’s ARI advantage therefore reflects the formation of some locally pure clusters, but its simultaneous FM collapse indicates that these clusters are globally misaligned with the true class structure, entire classes are systematically allocated to wrong partitions. SHOPCA’s lower ARI in these isolated cases is more than compensated by substantially higher FM, signaling that its curvature-regularized linear projection preserves the global class topology more faithfully than UMAP’s iterative graph optimization, which overfits to spurious local neighborhood structure when samples are scarce. The most dramatic evidence of this robustness gap appears in datasets where UMAP’s embedding degenerates entirely. On mammography, UMAP yields an ARI of 0.0297, worse than random assignment, while SHOPCA rescues the clustering to 0.3127 and achieves an FM of 0.9807. Similar catastrophic failures for UMAP occur on sonar (ARI −0.0040), blood-transfusion-service-center (ARI −0.0521), Breast (ARI −0.0378), and irish (ARI −0.0004), all of which SHOPCA elevates to positive, meaningful scores. These are not incremental improvements but qualitative phase transitions: UMAP’s stochastic neighbor embedding, deprived of sufficient data to estimate a reliable local connectivity graph, produces embeddings that actively destroy class-discriminative structure, whereas SHOPCA’s mean shape operator, computed from closed-form local covariances, remains stable because it does not depend on graph connectivity or iterative optimization convergence. The consistency of SHOPCA’s advantage in the Fowlkes-Mallows index is particularly noteworthy. FM is the most stringent of the three metrics because it requires both high precision (clusters are pure) and high recall (classes are recovered completely). UMAP’s mean FM across the 28 datasets is 0.4835, barely above the midpoint of the [0,1][0,1] range, whereas SHOPCA reaches 0.5971, with median FM of 0.5894 compared to UMAP’s 0.5028. This 17% median gap in FM, combined with the 135% median gap in ARI (0.2086 vs. 0.0887), confirms that SHOPCA’s curvature-aware regularization produces partitions that are not merely different from UMAP’s, but structurally superior in the sense of simultaneously capturing local cohesion and global class correspondence. The selected regularization parameters λ∗λ^* span an even wider dynamic range than in the second experiment, from 0.32 on ThreeOf9 to 1615.31 on ecoli, reinforcing the data-adaptive nature of the variance-curvature trade-off. The fact that no dataset selects λ∗λ^*, even in this small-sample regime where one might naively expect variance to dominate, confirms that curvature information remains extractable and beneficial even from severely undersampled neighborhoods. Collectively, these results position SHOPCA not merely as a fallback option when UMAP fails, but as a principled, computationally efficient alternative that dominates the state of the art in precisely the small-sample regimes where nonlinear graph-based methods are most needed yet least reliable. To complement the quantitative analysis, Figure 4 presents a qualitative comparison of the two-dimensional embeddings produced by UMAP (with default parameters) and the proposed SHOPCA on three representative small-sample datasets: UMIST_Faces_Cropped and mammography. UMAP approximates manifold topology through an iterative optimization of fuzzy simplicial sets, a strategy that presupposes sufficient data density to reliably estimate local connectivity and construct a faithful neighborhood graph. In low-density regions or small-sample regimes, this assumption breaks down: the underlying k-nearest neighbor graph becomes sparse and geometrically inconsistent, causing the stochastic cross-entropy optimization to overfit to spurious local connectivity patterns or collapse into degenerate configurations where distinct classes are artificially fused or fragmented. Visually, this manifests as irregular, diffused cluster boundaries and scattered cluster structure, which is not a realistic feature. In contrast, SHOPCA yields markedly more compact and well-delineated cluster arrangements. Because it operates in a single pass, computing local curvature through closed-form covariance eigendecompositions rather than iteratively optimizing a fuzzy topological representation, SHOPCA remains stable even when the number of observations is small to support reliable graph inference. The curvature-aware regularization steers the linear projection toward directions where the manifold folds and bends, producing isotropic point clouds whose natural density boundaries align with the ground-truth classes. This geometric regularity can be particularly advantageous for classical clustering algorithms such as k-means, agglomerative clustering with Ward’s linkage, and Gaussian mixture models, which implicitly assume that clusters occupy convex, tightly packed regions of the feature space. The visual evidence therefore corroborates the statistical findings: by replacing iterative graph-based topology inference with a principled, closed-form curvature regularizer, SHOPCA transforms the latent space into a representation that is structurally more amenable to unsupervised class discovery when data is scarce. Figure 4: Qualitative comparison of the two-dimensional embeddings produced by standard PCA and the proposed SHOPCA on representative datasets: UMIST_Faces_Cropped (first row) and mammography (second row). 6 Conclusions This paper introduced SHOPCA, a linear, closed-form generalization of PCA that regularizes the classical covariance matrix with a differential-geometric curvature term, the mean shape operator ¯ S, estimated locally from k-nearest-neighbor neighborhoods and aggregated into a single d×d× d matrix. The central methodological contribution is not merely the regularization itself, Eq. equation 13, but the fact that every step of the resulting pipeline, representation learning and the selection of the regularization strength λ, can be carried out without recourse to class labels. The eigengap criterion of Section 3.1 closes what was, in earlier formulations of this idea, the one remaining supervised step of an otherwise unsupervised method, and Section 4 placed this construction, and the resulting spectral behavior of Σcurv(λ) _curv(λ), on a precise theoretical footing: as an instance of geometry-informed covariance shrinkage (Section 4.1), with a provably real local spectrum whose realness is only approximately preserved under aggregation (Section 4.2), an exact variational characterization of its principal directions (Section 4.3), a quantitative Bauer–Fike bound on eigenvalue stability under regularization (Section 4.4), and an explicit reading as a regularized Mahalanobis metric that realizes, in closed form, the co-design of projection and metric motivated in the introduction. The empirical evidence across all three experiments consistently supports the central hypothesis that local curvature carries discriminative structure beyond what variance alone can capture. Against standard PCA on 30 real-world datasets (Section 5), SHOPCA achieved a strict, dataset-wide sweep across all three external validity indices, with average relative gains of 212% in ARI, 73% in FM, and 140% in V-measure – and, notably, no dataset selected λ∗=0λ^*=0, indicating that curvature information was found useful essentially everywhere in the benchmark suite. Against Isomap, a nonlinear method with strictly greater representational capacity, SHOPCA still achieved higher aggregate scores on all three indices and a clean three-metric sweep on 16 of 25 datasets, with Isomap’s failures concentrated precisely where its neighborhood-graph construction is known to be fragile. Against UMAP in the small-sample regime, where graph-based nonlinear methods are least reliable, SHOPCA’s advantage was most pronounced, with a 96% relative improvement in mean ARI and a three-metric sweep on 82% of the evaluated datasets. Taken together, these results indicate that a substantial fraction of the representational benefit typically attributed to nonlinear, iterative manifold learning can be recovered by a linear, closed-form method, provided the covariance structure is regularized with the correct second-order geometric information, and that this recovery does not require the computational overhead, hyperparameter sensitivity, or graph-construction fragility of its nonlinear counterparts. These strengths should be read alongside the method’s limitations, several of which are direct consequences of design choices made explicit in Sections 3–4. The shape operator estimator’s O(nd2(k+d2))O(nd^2(k+d^2)) complexity in the ambient dimension d (Section 2.2) remains the dominant computational cost for high-dimensional raw data, mitigated but not eliminated by the pre-PCA step at T=25T=25; datasets whose informative structure survives poorly under this initial linear compression are not well served by the present pipeline. The formal treatment of the shape operator assumes a hypersurface, codimension-one setting for notational clarity, and its extension to the general codimension d−m>1d-m>1 regime used operationally throughout the paper is handled through the local PCA frame rather than through a fully general, normal-bundle-valued second fundamental form, a gap between the operational estimator and its formal justification that Proposition 1 and the surrounding discussion make explicit rather than obscure. Relatedly, Proposition 2 shows that the variational reading of SHOPCA’s principal directions applies exactly to a symmetrized operator that the algorithm does not itself diagonalize, leaving an open, quantifiable gap between the interpretive story and the implementation whenever ¯ S’s skew-symmetric part is non-negligible. Finally, while λ is now selected in a fully unsupervised manner, the neighborhood size k and the target embedding dimension m remain fixed by heuristic or by problem specification rather than selected by an analogous label-free criterion. Several directions for future work follow directly from these observations. First, the label-free selection machinery developed for λ could be extended to jointly select the neighborhood size k and, where unknown, the target embedding dimension m, and cross-validated against a bootstrap/consensus stability criterion as an independent unsupervised signal, providing a fully self-contained, hyperparameter-free pipeline. Second, the symmetrized variant Σcurvsym(λ) _curv^sym(λ) identified in Proposition 2, whose eigenvectors are guaranteed real and orthogonal by the spectral theorem, warrants a direct empirical comparison against the asymmetric formulation used here, to determine whether closing the theory-implementation gap costs or improves clustering performance in practice. Third, the current formulation aggregates local shape operators via a simple, uniform average; locally adaptive weighting schemes that emphasize regions of high geometric complexity may yield sharper, more discriminative projections on datasets with strongly non-uniform manifold structure. Fourth, the same curvature-aware regularization principle extends naturally beyond PCA to other linear frameworks, including Linear Discriminant Analysis in semi-supervised settings, Factor Analysis, and Canonical Correlation Analysis, and to the general-codimension, normal-bundle-valued formulation of the second fundamental form left open, which would place the present codimension-one treatment inside a fully general theoretical account. Fifth, a rigorous statistical convergence analysis of the mean shape operator estimator, as a function of sample size, neighborhood size k, and intrinsic dimension m, would place the empirical robustness reported in Section 5 on firmer asymptotic footing, complementing the finite-sample perturbation bounds of Section 4.4. Sixth, integrating SHOPCA’s shape-operator regularization into deep architectures, as a geometry-aware initialization or feature-extraction layer, represents a natural bridge between the closed-form statistical approach developed here and the broader geometric deep learning movement. Finally, domain-specific applications in which the interplay between variance and curvature is expected to be particularly rich, such as hyperspectral imaging, single-cell transcriptomics, and medical imaging, offer a compelling avenue for further validation of the method’s practical impact beyond the benchmark suite considered in this paper. 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/ShapeOperatorPCA. Data availability All datasets used in the experiments are publicly available at w.openml.org. References Bauer and Fike (1960) F. L. Bauer and C. T. Fike Norms and exclusion theorems. Numerische Mathematik 2, p. 137–141. External Links: Document Cited by: §4.4, §4.4. Bellet et al. (2015) A. Bellet, A. Habrard, and M. Sebban Metric learning. Synthesis Lectures on Artificial Intelligence and Machine Learning, Springer, Cham. External Links: Document Cited by: §1, §2. Bronstein et al. (2021) M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković Geometric deep learning: grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2104.13478. External Links: Document Cited by: §1, §1. Bronstein et al. (2017) M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst Geometric deep learning: going beyond Euclidean data. IEEE Signal Processing Magazine 34 (4), p. 18–42. External Links: Document Cited by: §1, §1. Cai and He (2022) D. Cai and X. He Unsupervised feature selection for multi-cluster data. Pattern Recognition 125, p. 108536. External Links: Document Cited by: §5.3. Davis and Kahan (1970) C. Davis and W. M. Kahan The rotation of eigenvectors by a perturbation. I. SIAM Journal on Numerical Analysis 7 (1), p. 1–46. External Links: Document Cited by: §3.1. de Ocáriz Borde and Bronstein (2025) H. S. de Ocáriz Borde and M. Bronstein Mathematical foundations of geometric deep learning. External Links: 2508.02723, Link Cited by: §1. do Carmo (1992) M. P. do Carmo Riemannian geometry. Birkhäuser, Boston. Cited by: §2, §2, §2. do Carmo (2016) M. P. do Carmo Differential geometry of curves and surfaces: revised and updated second edition. Dover Publications, New York, NY. External Links: ISBN 978-0-486-81797-2 Cited by: §2, §2, §2, §2, §2, §4.2. Donoho and Grimes (2003) D. L. Donoho and C. Grimes Hessian eigenmaps: locally linear embedding techniques for high-dimensional data. Proceedings of the National Academy of Sciences 100 (10), p. 5591–5596. External Links: Document Cited by: §2.2. Fowlkes and Mallows (1983) E. B. Fowlkes and C. L. Mallows A method for comparing two hierarchical clusterings. Journal of the American Statistical Association 78 (383), p. 553–569. External Links: Document Cited by: §5. Ghojogh et al. (2023) B. Ghojogh, M. Crowley, F. Karray, and A. Ghodsi Elements of dimensionality reduction and manifold learning. Springer International Publishing, Cham. External Links: ISBN 978-3-031-10601-9, Document Cited by: §1. Harandi et al. (2017) M. Harandi, M. Salzmann, and R. Hartley Joint dimensionality reduction and metric learning: a geometric take. In Proceedings of the 34th International Conference on Machine Learning, D. Precup and Y. W. Teh (Eds.), Proceedings of Machine Learning Research, Vol. 70, p. 1404–1413. Cited by: §1, §1, §4.5, §4.5. Horn and Johnson (2013) R. A. Horn and C. R. Johnson Matrix analysis. 2nd edition, Cambridge University Press, New York, NY. External Links: Document Cited by: §4.2, §4.3. Hubert and Arabie (1985) L. Hubert and P. Arabie Comparing partitions. Journal of Classification 2 (1), p. 193–218. External Links: Document Cited by: §5. Jolliffe and Cadima (2016) I. T. Jolliffe and J. Cadima Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374 (2065), p. 20150202. External Links: Document Cited by: §1, §1. Kaya and Bilge (2019) M. Kaya and H. Ş. Bilge Deep metric learning: a survey. Symmetry 11 (9), p. 1066. External Links: Document Cited by: §1, §2. Kulis (2013) B. Kulis Metric learning: a survey. Foundations and Trends in Machine Learning 5 (4), p. 287–364. External Links: Document Cited by: §1. Ledoit and Wolf (2004) O. Ledoit and M. Wolf A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 88 (2), p. 365–411. External Links: Document Cited by: §4.1. Lee (2012) J. M. Lee Introduction to smooth manifolds. 2nd edition, Graduate Texts in Mathematics, Vol. 218, Springer, New York. External Links: Document Cited by: §2, §2, §2. McInnes et al. (2018) L. McInnes, J. Healy, and J. Melville UMAP: uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426. External Links: Link Cited by: §1, §1, §5.3, §5. Ng et al. (2001) A. Y. Ng, M. I. Jordan, and Y. Weiss On spectral clustering: analysis and an algorithm. In Advances in Neural Information Processing Systems (NeurIPS), Vol. 14, p. 849–856. Cited by: §3.1, §3.1. O’Neill (2006) B. O’Neill Elementary differential geometry. revised 2nd edition, Academic Press, Amsterdam. Cited by: §2, §2. Papillon et al. (2025) M. Papillon, S. Sanborn, J. Mathe, L. Cornelis, A. Bertics, D. Buracas, H. J. Lillemark, C. Shewmake, F. Dinc, X. Pennec, and N. Miolane Beyond Euclid: an illustrated guide to modern machine learning with geometric, topological, and algebraic structures. Machine Learning: Science and Technology 6. External Links: Document Cited by: §1. Rosenberg and Hirschberg (2007) A. Rosenberg and J. Hirschberg V-Measure: a conditional entropy-based external cluster evaluation measure. In Proceedings of the 2007 Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning, p. 410–420. External Links: Link Cited by: §5. Strehl and Ghosh (2002) A. Strehl and J. Ghosh Cluster ensembles — a knowledge reuse framework for combining multiple partitions. In Journal of Machine Learning Research, Vol. 3, p. 583–617. External Links: Link Cited by: §5. Tenenbaum et al. (2000) J. B. Tenenbaum, V. de Silva, and J. C. Langford A global geometric framework for nonlinear dimensionality reduction. Science 290 (5500), p. 2319–2323. External Links: Document Cited by: §5.2, §5. Vanschoren et al. (2013) J. Vanschoren, J. N. van Rijn, B. Bischl, and L. Torgo OpenML: networked science in machine learning. ACM SIGKDD Explorations Newsletter 15 (2), p. 49–60. External Links: Document Cited by: §5. von Luxburg (2007) U. von Luxburg A tutorial on spectral clustering. Statistics and Computing 17 (4), p. 395–416. External Links: Document Cited by: §3.1, §3.1. Wang and Sun (2015) F. Wang and J. Sun Survey on distance metric learning and dimensionality reduction in data mining. Data Min. Knowl. Discov. 29 (2), p. 534–564. External Links: ISSN 1384-5810, Document Cited by: §1, §4.5. Ward (1963) J. H. Ward Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association 58 (301), p. 236–244. External Links: Document Cited by: §5. Weber (2025) M. Weber Geometric machine learning. AI Magazine 46 (1), p. e12210. External Links: Document Cited by: §1. Yu et al. (2015) Y. Yu, T. Wang, and R. J. Samworth A useful variant of the Davis–Kahan theorem for statisticians. Biometrika 102 (2), p. 315–323. External Links: Document Cited by: §3.1. Zhou (2021) Z. Zhou Dimensionality reduction and metric learning. In Machine Learning, External Links: Document Cited by: §1.