The decomposition of the higher-order homology embedding constructed from the k-Laplacian
Year: 2021

Abstract:
The null space of the k-th order Laplacian Lk, known as the k-th homology vector space, encodes the non-trivial topology of a manifold or a network. Understanding the structure of the homology embedding can thus disclose geometric or topological information from the data. The study of the null space embedding of the graph Laplacian L0 has spurred new research and applications, such as spectral clustering algorithms with theoretical guarantees and estimators of the Stochastic Block Model. In this work, we investigate the geometry of the k-th homology embedding and focus on cases reminiscent of spectral clustering. Namely, we analyze the connected sum of manifolds as a perturbation of the direct sum of their homology embeddings. We propose an algorithm to factorize the homology embedding into subspaces corresponding to a manifold’s simplest topological components. The proposed framework is applied to the shortest homologous loop detection problem, a problem known to be NP-hard in general. Our spectral loop detection algorithm complements existing methods of topological data analysis. It scales better than existing methods with no assumptions on the structure of data and is effective on diverse data such as point clouds and images.

heading: 1 Motivation
text: The k-th homology vector space Hk provides rich geometric information on manifolds/networks. For instance, the zeroth, the first, and the second homology vector spaces identify the connected components, the loops, and the cavities in the manifold, respectively. Topological Data Analysis (TDA) [24, 58] (as well as other early works in this field) has found wide use in analyzing biological [31, 47], human behavior [1, 64], or other complex systems [58]. Even though they easily generalize to k ≥ 1, additional efforts are needed to extract topological features (e.g., instances of loops) besides ranks due to the combinatorial complexity of the structures that support them.
Spectral methods based on k-Laplacians (Lk), by contrast, investigate Hk in a linear algebraic manner; abundant geometric information can be extracted from the homology embedding Y (the null space eigenvectors of Lk) of Hk. Analysis of the eigenfunctions (of H0) [14, 37, 40, 49] of the graph Laplacian L0 is pivotal in providing guarantees for spectral clustering and community detection algorithms in well-separated datasets. Recent advances in this field [4, 12, 48] extend the existing spectral algorithms based on L0 to k ≥ 1; however, the theoretical analysis of Y of Hk, unlike spectral clustering, is less developed, in spite of intriguing empirical results by [21]. Here, we put these observations on a formal footing based on the concepts of connected sum and prime decomposition of manifolds (Section 2 and 3). We examine these operations through the lens of the (subspace) perturbation to the homology embedding Y of the discrete k-Laplacian Lk (an estimator of the continuous k-Hodge Laplacian operator ∆k) on finite samples (Section 4). This framework
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
finds applications in, i.e., identifying the shortest homologous loops (Section 5). Lastly, we support our theoretical claims with numerous empirical results from point clouds and images.

heading: 2 Background in Hodge theory and topology
text: Simplicial and cubical complex. An abstract complex is a natural extension of a graph designed to capture higher-order relationships between its vertices. A simplicial k-complex (used when the data are point clouds or networks) is a tuple SCk = (Σ0, · · · ,Σk), with Σℓ being a set of ℓ dimensional simplices, such that every face of a simplex σ ∈ Σℓ is in Σℓ−1 for ℓ ≤ k. As a side note, a graph G = (V,E) is an SC1; and SC2 = (V,E, T ), which is commonly used in edge flow learning [12, 48], is obtained by adding a set of 3-cliques (triangles) T of G. This procedure extends to defining Σℓ as the set of all ℓ-cliques of G, with the resulting complex called a clique complex of the graph G. This complex is also known as a Vietoris-Rips (VR) complex if G is the ϵ-radius neighborhood graph used in the manifold learning literature [12, 14, 53], The cubical k-complex CBk = (K0, · · · ,Kk) is a complex widely used with image data. The difference between this complex and the SCk is that a CBk is a collection of sets of ℓ-cubes, for ℓ < k. Note that we write Σ0 = K0 = V the vertex set and Σ1 = K1 = E the edge set. Σ2 = T and K2 = R are the triangle and rectangle set, respectively. Additionally, we define nℓ = |Σℓ| (or = |Kℓ|) to be the cardinality of the ℓ-dimensional cells and let n = n0 for simplicity. For more information about building various complexes on different datasets please refer to Otter et al. [43].
k-cochain. By choosing an orientation for every k-simplex σk,i ∈ Σk (or Kk), one can define a finite-dimensional vector space Ck (k-cochain space1). An element ωk = ∑ i ωk(σk,i)σk,i ∈ Ck is called a k-cochain; one can further express ωk as ωk = (ωk,1, · · · , ωk,nk)⊤ ∈ Rnk by identifying each σk,i with the standard basis vector ei ∈ Rnk . Functions on nodes and edge flows, for example, are elements of C0 and C1, respectively.
Boundary matrix. The k-th boundary matrix Bk [35] maps a k-cochain of k-cells (simplices/cubes) σk to the (k− 1)-cochain of its faces, i.e., Bk : Ck → Ck−1. Bk ∈ {0,±1}nk−1×nk is a sparse binary matrix, with the sign of the non-zero entries σk−1, σk given by the orientation of σk w.r.t. its face σk−1. Hence, different SC or CB will induce different Bk. For k = 1 on either the SC or CB, the boundary map is the graph incidence matrix, i.e., (B1)[x],[x,y] = 1, (B1)[y],[x,y] = −1, and zero otherwise; for k = 2, each column of B2 contains the orientation of a triangle/rectangle w.r.t. its edges. Specifically, for an SC, (B2)[x,y],[x,y,z] = (B2)[y,z],[x,y,z] = 1, (B2)[x,z],[x,y,z] = −1, and 0 otherwise; for a CB, (B2)[x,y],[x,y,z,w] = (B2)[y,z],[x,y,z,w] = (B2)[z,w],[x,y,z,w] = 1, (B2)[x,w],[x,y,z,w] = −1, and 0 otherwise. Simplex σk+1 is a coface of σk iff σk is a face of σk+1; let coface(σk) be the set of all cofaces of σk. The (k− 1)-th coboundary matrix B⊤k (adjoint of Bk) maps σk−1, as a (k − 1)-cochain, to the k-cochain of coface(σk−1).
k-Hodge Laplacian. Let Wℓ be a diagonal non-negative weight matrix of dimension nℓ, with [Wℓ]σ,σ representing the weight of the ℓ-simplex/cube σ and wℓ ← diag(Wℓ). The weighted k-Hodge Laplacian2 [26] is defined as
Lk = A⊤k Ak +Ak+1A ⊤ k+1, where Aℓ = W −1/2 ℓ−1 BℓW 1/2 ℓ for ℓ = k, k + 1. (1)
The weights capture combinatorial or geometric information and must satisfy the consistency relation wℓ(σℓ) = ∑ σℓ+1∈coface(σℓ) wℓ+1(σℓ+1) (in matrix form: wℓ = |Bℓ+1|wℓ+1) for ℓ = k, k − 1. Hence Ak can be seen as normalized boundary matrix. To determine the weight for the (k + 1)- simplices, one can selected wk+1 to be constant [48] or based on (a product of) pairwise distance kernel (for k = 0, 1) so that the large sample limit exists [12, 14]. The first and second terms of (1) are called respectively the down (Ldownk = A ⊤ k Ak) and up (L up k = Ak+1A ⊤ k+1) Laplacians. For k = 0, the down component disappears and the resulting L0 is the symmetric normalized graph Laplacian used in spectral clustering [55] and Laplacian Eigenmap [5].
1We use chain and cochain interchangeably, see Lim [35] for the distinction between them. 2In this paper, these matrices are also called the k-Laplacian for simplicity.
k-th homology vector space and embedding. The homology vector space Hk is a subspace of Ck (loop space) such that every k-cycle (expressed as a k-cochain) in Hk is not the boundary of any (k+ 1)-cochain. In mathematical terms,Hk := ker(Ak)/im(Ak+1). The rank of the subspace is called the k-th Betti number βk = dim(Hk), which counts the number of “loops” (homology classes) in the SC. Hk is equivalent to the null space of Lk [35, 48]; therefore, a basis ofHk can be obtained by the eigenvectors Y = [y1, · · · ,yβk ] ∈ R
nk×βk of Lk with eigenvalue 0. The homology embedding maps a k-simplex σk to Yσk,: = [y1(σk), · · ·yβk(σk)]
⊤ ∈ Rβk . Note that the basis Y is only identifiable up to a unitary transformation; hence, individual homology embedding coordinates might change with a different basis Y.
Continuous operators on manifolds. The k-cochains are the discrete analogues of k-forms [59]. For k = 1, the following path integral [59] (along the geodesic γ(t) connecting x and y) relates a 1-cochain ω to a 1-form v (vector field): ω([x, y]) = ∫ y x v(γ(t))γ′(t)dt. To estimate a vector field from ω, one can solve a least-squares problem [12], which is the inverse operation of the path integral (e.g., the vector fields in Figure 1 are estimated from Y). Similarly, one can define the differential d and codifferential δ operators which are analog to B⊤k+1 and Bk, respectively. The (continuous) k-Hodge Laplacian operators, which act on k-forms, can be defined for manifolds too, i.e., by ∆k = dk−1δk + δk+1dk. The homology group (the continuous version ofHk) is defined as the null of ∆k. Its elements are harmonic k-forms ζk, computed by solving ∆kζk = 0 with proper boundary conditions; they represent the continuous version of the discrete homology basis Y.
Connected sum and manifold (prime) decomposition. The connected sum [33] of two d dimensional manifoldsM =M1♯M2 is built from removing two d dimensional “disks” from each manifoldM1,M2 and gluing together two manifolds at the boundaries (technical details in [33]). The analog of the connected sum for the abstract complexes will be defined in Section 3. The connected sum is a core operation in topology and is related to the concept of manifold (prime) decomposition. Informally speaking, the prime decomposition aims to factorize a manifoldM into κ smaller building blocks (M =M1♯ · · · ♯Mκ) so that eachMi cannot be further expressed as a connected sum of other manifolds. The well-known classification theorem of surfaces [3] states that any oriented and compact surface is the finite connected sum of manifolds homeomorphic to either a circle S1, a sphere S2, or a torus T2. Classification theorems for d > 2 are currently unknown; fortunately, the uniqueness of the prime decomposition for d = 3 (up to homeomorphism) was shown (Kneser-Milnor theorem [39]). Recently, Bokor et al. [9] (Corollary 2.5) showed the existence of factorizations of manifolds with d ≥ 5, even though they might not be unique. In this paper, we are interested in the following: given finite samples from M, which is a κ-fold connected sum ofMi, can this decomposition be recovered from the discrete homology embedding Y ofM? Namely, we would like to understand how Y relates to the homology embedding of each prime manifoldMi.

heading: 3 Definitions, theoretical/algorithmic aims, and prior works
text: Definitions. We assume that the data X is sampled from a d-dimensional oriented manifold M that can be decomposed into κ prime manifolds (M = M1♯ · · · ♯Mκ). Let Ii be an index set of the data points in X sampled from Mi, for i = 1, . . . κ. Denote by SCk, Lk, Hk(M), and βk the simplicial complex, the k-Laplacian, the k-homology space, and the k-th Betti number of X sampled fromM. Furthermore, let ŜC (i)
k = (Σ̂ (i) 0 , · · · , Σ̂ (i) k ), L̂
(ii)
k ,Hk(Mi), βk(Mi) be the same quantities for manifold Mi (supported on Ii for i ≤ κ). ŜCk and L̂k (without superscript i) are the comparable notations for the disjoint manifolds Mi’s, i.e. ŜC = ∪κi=1ŜC (i) = (Σ̂0, · · · , Σ̂k) with Σ̂ℓ = ∪κi=1Σ̂ (i) ℓ for ℓ ≤ k, and L̂k is a block diagonal matrix with the i-th block being L̂ (ii)
k . Additionally, let Y and Ŷ (both in Rnk×βk ) be the homology basis of Lk and L̂k, respectively. Let Si be the index set of columns of Ŷ corresponding to homology subspaceHk(Mi), with Si∩Sj = ∅ for i ̸= j, |Si| = βk(Mi), and S1∪· · ·∪Sκ = {1, · · · , βk}. Since Ŷ is the homology embedding of a block diagonal matrix L̂k, we choose it so that [Ŷ]σ,m equals the homology embedding of L̂ (ii) k if
σ ∈ Σ̂(i)k with column m ∈ Si and is zero otherwise. Namely, Ŷ lies in the direct sum of subspaces Hk(Mi) for i ≤ κ.
Theoretical aim. We are interested in the geometric properties of the null space eigenvectors Y, and specifically in recovering the homology basis Ŷ of the prime manifolds. Hence, we aim to bound the distance between the spaces spanned by Y and Ŷ. Under a small perturbation, one can provide an analogous argument to the orthogonal cone structure [40, 49] in spectral clustering (the zeroth homology embedding). The main technical challenge is that the connected sum of manifolds is a highly localized perturbation; namely, most cells are not affected at all, while those involved in the gluing process gain or lose O(1) (co)faces. Without properly designing Lk and L̂k, one might get a trivial bound.
Algorithmic aim. We exemplify the algorithmic aim using k = 1, d = 2, and κ = 2, particularly the genus-2 surface shown in Figure 1. The null space basis Y of Lk is only identifiable up to a unitary matrix due to the multiplicity of the zero eigenvalues. For instance, the top and bottom rows of Figure 1 are both valid bases for the harmonic edge flows inH1. However, the basis vector fields in the second row of Figure 1 are more inter-
pretable than those in the top row because Y (the first row) is a linear combination of Z (the second row), with each basis (column in the figure) corresponding to a single homology class (loop). Therefore, here we propose a data-driven approach to obtain the optimal basis Z such that the coupling from other manifolds/subspaces is as weak as possible. Being able to obtain Z from an arbitrary Y can support numerous applications (more in Section 5); however, it is difficult to design a criterion for finding the optimal Z without knowing the geometric structure of the embedding ofHk.
Prior work. The shape of the embedding by the principal eigenvectors of the graph Laplacian L0 is pivotal for showing the guarantees of spectral clustering algorithms for point cloud data or the inference algorithms for the stochastic block model. The analyses used either the matrix perturbation theory [40, 55, 56] or assume a mixture model [49]. For the higher-order k-Laplacian, it is reported empirically that the homology embedding is approximately distributed on the union (directed sum) of subspaces [21]; subspace clustering algorithms [30] were applied to partition edges/triangles under their framework.
On the application side, the eigen-embedding of Lk is used implicitly or explicitly in graph signal processing [27, 44, 45] and in learning algorithms utilizing Hodge decomposition [12, 28, 62]. They use the gradient/curl cochains (which correspond to the images of Ldownk and L up k , respectively) in addition to the homology embedding (ker(Lk)); hence, our framework is intrinsically different from these works.

heading: 4 Main result: connected sum as a matrix perturbation
text: In this section, we analyze the geometric structure of Y by viewing the operation of connected sum through the lens of matrix perturbation theory [51]. We show that, under certain conditions, the homology embedding Y of the joint Laplacian Lk is approximated by Ŷ for the simplices in Σk ∩ Σ̂k. In matrix terms, we show that Y ≈ ŶO (Theorem 1) with O a unitary transformation. We first prepare our assumptions suited for SC built from point clouds. Most of the assumptions (except Assumption 1 for which the connected sum might not be defined) can be extended to the clique complex (for networks) or cubical complex (for images) without too many modifications.
Assumption 1. The point cloud X ∈ Rn×D is sampled from a d-dimensional oriented and compact manifold M ⊆ RD; the homology vector spaces Hk(SC) formed by the simplicial complex constructed from X are isomorphic to the homology groupHk(M) ofM, i.e.,Hk(SC) ≃ Hk(M). Furthermore, assume thatM =M1♯ · · · ♯Mκ, and thatHk(ŜC (i) ) ≃ Hk(Mi) for i = 1, · · · , κ.
This assumption is the minimal assumption needed for the analysis of the embedding of the Lk; it states that any procedure to construct the simplicial complex or weight function for Lk is accepted as long as the isomorphic condition holds. The construction of the SC from the point cloud is out of the scope of this manuscript (see, e.g., Chen et al. [12] for building L1 from X with an analyzable limit and Latschev’s theorem [32] on VR complexes). The last condition requires that the manifoldM can be decomposed; this is most likely true, except for the known hard case ofM with d = 4 discussed in Section 2. Note that Assumption 1 is similar to those used in the multi-manifold clustering [54], where they required the manifold to be constructed by a κ-fold unionM =M1 ∪ · · · ∪Mκ. The goal of our framework is to analyze the homology embedding of k-Laplacian. By contrast, Trillos et al. [54] are interested in identifying each manifold with the notion of clustering (H0). Assumption 1 is for points sampled from manifold M only. To make this assumption hold for networks or images, one can require that the Lk constructed from these two datasets can be roughly factorized into block-diagonal entries. Below we provide two other assumptions that are valid for both SC and CB (with some modifications): the first one controls the eigengap and the second one ensures a small perturbation in the spectral norm of Lk − L̂k. By construction, Lk is positive semidefinite; since we are interested in the stability of its null space, we define, for any matrix L ⪰ 0, the eigengap as the smallest non-zero eigenvalue of L and denote it λmin(L).
Assumption 2. We denote the set of non-intersecting simplices be Nk := Σk ∩ Σ̂k. Let the set of destroyed and created k-simplices during connected sum by Dk and Ck, respectively; they are defined by Ck := Σk\Nk and Dk := Σ̂k\Nk. We have: (1) no k-homology class is created during the connected sum process, i.e., βk(SC) = ∑κ i=1 βk(ŜC (i) ). (2) The eigengaps of LC,Ck and L̂ D,D
k are bounded away from the eigengaps of L (ii) k , i.e., min{λmin(L C,C k ), λmin(L̂
D,D
k )} ≫ min{δ1, · · · , δκ}, where δi is the eigengap of L(ii)k .
The first condition requires that the intersecting simplices Dk ∪ Ck do not create or destroy any k-th homology class; this holds, for instance, when the manifoldM has dimension d > k. Under this condition, we have Hk(M1♯M2) ≃ Hk(M1) ⊕ Hk(M2) [33]. A counterexample for this condition is, e.g., inspecting the cavity space (k = 2) of a genus-2 surface built from gluing two tori together. That is, β2 of a genus-2 surface is 1, while the sum of β2 of two tori is 2. The second condition requires that the principal submatrix of Lk described by the block of Ck ∪ Dk has large eigengap. This happens, e.g., when Ck and Dk are cliques and are contained in small balls.
Assumption 3 (Informal). Let w̃k = |Bk+1[Nk,Nk+1]|wk+1, w̃k−1 = |Bk[:,Nk]|w̃k. For ℓ = k or k − 1, we have maxσ∈Nℓ {wℓ(σ)/w̃ℓ(σ)− 1} ≤ ϵℓ, maxσ∈Nℓ {ŵℓ(σ)/w̃ℓ(σ)− 1} ≤ ϵℓ, and maxσ∈Nℓ {|wℓ(σ)/ŵℓ(σ)− 1|} ≤ ϵ′ℓ. Assumption S1 is the formal version of this assumption.
For k = 1, it states that not too many triangles are being created or destroyed during connected sum. For this assumption to hold, the sampling density in the connected sum region should be smaller than in other regions, i.e., the manifoldM should be sparsely connected (e.g., Figure 2a). Empirically, we observed that the perturbation is small even whenM is not sparsely connected (more discussions in Section 6). Note also that ϵ′ℓ ≪ ϵℓ, for ϵ′ℓ represents the net change in the degree after connected sum. It might be possible to obtain a tighter bound fully by ϵ′ℓ’s, which do not depend on the relative density between the connected sum region and the remaining manifolds; we leave it as future work.
Theorem 1. Let DiffLdownk (DiffL up k ) be the modified difference (defined in Supplement A) of L down k and L̂ down
k (respectively L up k and L̂
up k ). Under Assumptions 1–3 with notations defined as before and λk = k + 2, if ∥∥∥DiffLdownk ∥∥∥2 ≤ [2√ϵ′k + ϵ′k + (1 +√ϵ′k)2 √ϵ′k−1 + 4√ϵk−1]2 λ2k−1 and
∥DiffLupk ∥ 2 ≤ [ 2 √ ϵ′k + ϵ ′ k + 2ϵk + 4 √ ϵk ]2
λ2k, then there exists a unitary matrix O ∈ Rβk×βk such that
∥∥∥YNk,: − ŶNk,:O∥∥∥2 F ≤ 8βk [∥∥∥DiffLdownk ∥∥∥2 + ∥DiffLupk ∥2] min{δ1, · · · , δκ} . (2)
The proof (in Supplement A) is based on bounding the error between Lk and L̂k with L̃k (the Laplacian after removal of k-simplices during connected sum), the use of a variant of the Davis-
Kahan theorem [61], and the bound of the spectral norm of Lk for a simplicial complex, i.e., ∥Lk∥2 ≤ λk = k + 2 [26].
Algorithm 1: Subspace identification Input :SC, k, weights Wk+1
1 Bk,Bk+1 = BOUNDARYMAPS(SC, k) ▷ in Algorithm S2
2 for ℓ = k, k − 1 do 3 Wℓ ← diag{|Bℓ+1|Wℓ+11nℓ+1} 4 Aℓ+1 ←W−1/2ℓ Bℓ+1W 1/2 ℓ+1 5 Lk = A⊤k Ak +Ak+1A ⊤ k+1 6 Y ∈ Rnk×βk ← NULLSPACE(Lk) 7 Z← ICANOPREWHITE(Y)
Return :Independent basis Z
The LHS of (2) contains only the simplices in Nk because Y is not defined on Dk (similar for Ŷ and Ck). Nonetheless, (2) makes sure that the (unbounded) perturbations in the embedding of Ck ∪ Dk do not propagate to the rest of the simplices. The bound in (2) can be extended to CB (Corollary 2) by changing the λk value from (k+2) to 2k+2. The 2k+2 term here is the maximum eigenvalue of the Lk built from any cubical complex (Proposition S3). Corollary 2 (For Lk built from a CB). Under Assumptions 2–3 with DiffLupk as well as DiffL down k defined in Theorem 1 and λk = 2k + 2, there exists a unitary matrix O such that (2) holds.
Subspace Hk(Mi) identification. We propose to approximately recover Ŷ from the columns of the coupled basis Y by blind source separation, as described by Algorithm 1. Specifically, the independent basis Z is obtained by Infomax ICA [6] on Y of Lk, with a modification (Line 7) that preserves the necessary properties of harmonic cochains (i.e., they are divergence-free and curl-free, see also Proposition 3). Algorithm 1 works for CB as well by using the appropriate Bk, Bk+1 construction method (Line 1). Each column of the obtained Z is now approximately supported on a single prime manifoldMi.

heading: 5 Applications: homologous loops detection, clustering, and visualization
text: Algorithm 2: Spectral homologous loop detection Input :Z = [z1, · · · , zβ1 ], V , E, edge distance d
1 for i = 1, · · · , β1 do 2 E+i ← {(s, t) : (s, t) ∈ E and [zi](s,t) > 0} 3 E−i ← {(t, s) : (s, t) ∈ E and [zi](s,t) < 0} 4 τ ← PERCENTILE(|zi|, 1− 1/β1) 5 E×i ← {e ∈ E + i ∪ E − i : |[zi]e| < τ} 6 Ei ← E+i ∪ E − i \E × i 7 Gi ← (V,Ei), with weight of e ∈ Ei being [d]e 8 dmin = inf 9 for e = (t, s0) ∈ Ei do
10 P∗, d∗ ← DIJKSTRA(Gi, from=s0, to=t) ▷ Note that P∗ = [s0, s1, · · · , t] 11 if d∗ < dmin then 12 Ci ← [t, s0, s1, · · · , t]
Return : C1, · · · , Cβ1 Homologous loop detection. In addition to the rank information, one might find it beneficial to extract the shortest cycle of the corresponding Hk generator. This application is found useful in domains including finding minimum energy trajectories in molecular dynamics datasets, trajectory inference in RNA single-cell sequencing [46], and segmenting circular structures in medical images [50]. We propose a spectral shortest homologous loop detection algorithm (Algorithm 2) based on the shortest path algorithm (Dijkstra) as follows: for each dimension i = 1, · · · , β1, the algorithm reverses every edge e having negative [zi]e to generate a weighted digraph Gi = (V,Ei) (Lines 2–4), with the weight of edge e = (j, j′) ∈ Ei equal to the Euclidean distance [d](j,j′) = ∥xj − xj′∥2. The algorithm finds a shortest (in terms of d) loop on this weighted digraph for each i and outputs it as the homologous loop representing the i-th class. We present the following proposition (with the proof in Supplement B) to support Algorithm 2; it implies that if each coordinate of Z extracted from Algorithm 1 corresponds to a homology class, then the detected homologous loop for each homology class is the shortest. Proposition 3. Let zi for i = 1, · · · , β1 be the i-th homology basis that corresponds to the i-th homology class. For every i = 1, · · · , β1, (1) there exist at least one cycle in the digraph Gi such that every vertex v ∈ V can traverse back to itself (reachable); (2) the corresponding cycle will enclose at least one homology class (no short-circuiting).
Since every vertex is reachable from itself, we are guaranteed to find a loop for any starting/ending pair (Lines 9–12). Additionally, there will be no short-circuiting for any loop; each loop we found from Dijkstra is guaranteed to be non-trivial. However, there is one caveat from the second property: even though the i-th loop is non-trivial, it might not always be corresponding to the i-th homology class due to the noise in small [zi]e. Namely, loops that do not represent i-th homology class can be formed with edges e having small [zi]e, resulting in the instability and the (possible) duplication of the identified loops. To address the issue, we propose a heuristic thresholding, by which we keep the n1/β1 edges with the largest absolute value in |[zi]e| (Lines 4–5). We chose to keep n1/β1 by treating each homology class equally, i.e., each class has roughly n1/β1 edges. Unlike Theorem 1 or Algorithm 1, Algorithm 2 cannot be extended to the case when k ≥ 2 because the Dijkstra algorithm is employed. We leave its generalization to extract higher-order cavities as future work.
Compared with the previous approach that finds the shortest loops [16] combinatorially, our approach has better time complexity; specifically, the algorithm by Dey et al. [16] has time complexity O(nn31 + nn21n2), whereas Algorithm 2 runs in time O(n2.37···1 + β21n1 + β1n1n log n). The first, second, and third terms correspond to the time complexity of eigendecomposition of L1, the Infomax ICA, and the Dijkstra algorithm on every digraph Gi, respectively. Note that if the simplicial complex is built from point clouds, the number of triangles n2 may be large; this dependency on n2 makes the algorithm [16] hard to scale. On the other hand, our framework requires that zi are each supported on one homology class; therefore, loops can only be correctly identified using Algorithm 2 if the manifold is sparsely connected (Assumptions 1–3). Additional comparisons between our algorithm and other methods that pose special requirements on the analyzed data can be found in Supplement D.4.2.
Classifying any 2-dimensional manifold. The Betti number β1 of a torus is 2, which is equal to that of two disjoint circles; hence one cannot distinguish these two manifolds only by rank information. Fortunately, they can be categorized using the homology embedding Z. By the classification theorem [3], any 2D surface is the connected sum of circles S1 and tori T1; therefore, Theorem 1 indicates that embedding lies approximately in the directed sum of homology subspace of S1 and/or T2. The homology embedding of S1 is a line since it is in R1. On the other hand, any loop in a torus can be a convex combination of the two homology classes, implying that the intrinsic dimension of the homology embedding is 2. It is hard to obtain Z of any arbitrary torus; we present the homology embedding of the flat m-torus below by expressing the null space basis (1-cochains) as the path integrals of the corresponding harmonic 1-forms [12, 59].
Proposition 4. The envelope of the first homology embedding (1-cochain) induced by the harmonic 1-form on the flat m-torus Tm is an m-dimensional ellipsoid.
The proof (in Supplement B) is straightforward thus is omitted here. Proposition 4 and the classification theorem suggest that the first homology embedding is either a line, a disk, or a combination of the two (with replacement). See an example for the genus-2 surface in Figures 2j and S1.
Note that it is possible to classify any 2-manifold with higher-dimensional homology groups; for instance, one can distinguish T2 from S1♯S1 by the second homology group. However, our approach scales better in computation or memory usage since inspectingH2 needs at least the calculation and storage of tetrahedrons from the neighborhood graph.
Other applications. As pointed out earlier, one can visualize the basis of the harmonic vector fields (of Hk) by overlaying the columns of Y onto the original dataset (Figure 1). Being able to successfully extract a decoupled basis Z increases the interpretability of Hk, as shown in the second row of Figure 1. Theorem 1 also supports the use of subspace clustering algorithm in the higher-order simplex clustering framework [21].

heading: 6 Experiments
text: We demonstrate our approach by computing Y,Z and the shortest loops for five synthetic manifolds: two of them are prime manifolds (TORUS torus, 3-TORUS three-torus) and three (PUNCTPLANE punctured plane with two holes, GENUS-2 genus-2, and TORI-CONCAT concatenation of 4 tori) are factorizable manifolds. Furthermore, five additional real point clouds (ETH and MDA from chemistry, PANCREAS from biology, 3D-GRAPH from 3D modeling, and ISLAND from oceanography) are ana-
lyzed under this framework. For all the point clouds, we build the VR complex SC from the CkNN kernel [8] so that the resulting L1 is sparse and the topological information is preserved. Note that other methods for building an SC from X can also be used as long as Hk is successfully identified (Assumption 1). Lastly, we illustrate the efficacy of our framework to a non-manifold data: RETINA from medical imaging. Please refer to Supplement D for detailed discussions on procedures to generate, preprocess, and download these datasets. All experiments are replicated more than five times with similar results. We perform our analysis on a desktop running Linux with 32GB RAM and an 8-Core 4.20GHz Intel® Core™ i7-7700K CPU; every experiment completes within 3 minutes (1-2 minutes on eigendecomposition of L1, and around 30 seconds on both ICA and Algorithm 2).
Synthetic manifolds. The results for the synthetic manifolds are in Figure 2. Figure 2a (the harmonic embedding of PUNCTPLANE) confirms Theorem 1 that Y is approximately distributed on two subspaces (yellow and red), with each loop parametrizing a single hole (inset of Figure 2a). As discussed previously in Figure 1, the harmonic vector bases (green and blue) are mixtures of the separate subspaces; therefore, these bases have poor interpretability compared with the independent subspace Z identified by Algorithm 1. The shortest loops (Figure 2b) corresponding to z1 (yellow), z2 (red) are obtained by running Dijkstra on the digraphs induced by z1 and z2 separately (Algorithm 2). Figures 2c–2f show the results of the two simple prime manifolds: TORUS and 3-TORUS. The harmonic embeddings of TORUS (Figure 2d) and 3-TORUS (Figure 2f) are a two-dimensional disk and a three-dimensional ellipsoid, respectively; this confirms the conclusion from Proposition 4. The shortest loops obtained from Algorithm 2 for these two datasets are in Figures 2c and 2e, showing that these loops travel around the holes in TORUS (or 3-TORUS). Note that we plot 3-TORUS in the intrinsic coordinate because a three torus can not be embedded in 3D without breaking neighborhood relationships. Three lines in 2e are indeed loops due to the periodic boundary condition, i.e., 0 = 2π, in the intrinsic coordinate. Figures 2h and 2j show the embedding of the coupled harmonic basis (Y) and that corresponding to the independent subspace (Z) obtained by Algorithm 1. Compared with Y, each coordinate of Z corresponds to a subspace, i.e., the left or right handle of GENUS-2, and does not couple with other homology generators. Z is thus a union of two 2D disks, with each disk approximating the harmonic embedding of a torus (see Figure S1 for more detail). Compared with the loops obtained by running Algorithm 2 on Y (Figure 2g), each loop in Figure 2i identified from Z parameterizes the corresponding homology generator without being homologous to other loops. Similar results on TORI-CONCAT are in Figures 2k and 2l, which correspond to the
loops obtained from Y and Z, respectively. The pairwise scatter plots of the eight-dimensional Z (or Y) are in Figure S2 of Supplement D. Note that PUNCTPLANE is an example of a sparsely connected manifold (see the low-density area in the middle), with ϵ1 ≈ 0.035 and ϵ0 ≈ 0.038. Manifolds of other synthetic/real datasets might not be sparsely connected due to the (approximately) constant sampling densities; nevertheless, the perturbations to the subspaces remain small for these datasets.
Small molecule data [13]. Figures 3a–3c and 3d–3f show our analysis on ETH and MDA, respectively. These two small molecule datasets, whose ambient dimensions are D = 102 and D = 98, are suggested to be noisy non-uniformly sampled tori [52]; the harmonic embeddings of these two datasets (Figures 3c and 3f ) confirm this idea. Finding the minimum trajectories corresponding to a specific bond torsion is of interest in chemistry; in these two molecular dynamics systems, this problem can be translated into finding the homologous loops in the point cloud. The homologous loops found by Algorithm 2 overlaid on the first three principal components (PCs) for these two datasets can be found in Figures 3b (for ETH) and 3e (for MDA). The same homologous loops plotted in the bond torsion space (with definition in the insets) based on our prior knowledge are in Figures 3a and 3d. Similar to the discussion for 3-TORUS (Figure 2e), the yellow/red trajectories form loops due to the periodic boundary condition of the bond torsions.
RNA single-cell sequencing data [7]. The trajectory inference methods [46] for analyzing the RNA single-cell sequencing datasets aim to order the cells (points in high-dimensional expression space) along developmental trajectories, which are inferred from the structure of the point clouds. Identifying loops in the dataset can serve as a building block for delineating a correct trajectory, especially for determining cell cycle and cell differentiation. To illustrate the idea, we compute the 1-Laplacian on the CkNN kernel [8] constructed on the UMAP [36] embedding (Algorithm 1). Figure 3g shows the identified loops from Algorithm 2, with the green loop being the cycle of ductal cells and yellow/red loops representing a trifurcation (endocrine cell differentiation).
Additional point cloud datasets. 3D-GRAPH [15] is a 3D model of a Buddha statue with a precomputed triangulation. We treat the 3D model as a point cloud and subsample 3000 farthest points from the original dataset; L1 is obtained from the VR complex of the CkNN kernel. Note that with this small sample size, two smaller loops near the waist of the statue are not detectable. Hence, the number of zero eigenvalues of L1 is 3, with the corresponding homology generators shown in Figure 3h. ISLAND [23], which contains ocean buoys around the Tasman sea, is the other point cloud in
our analysis. The estimated β1 is 3, with the detected loops being the North Island of New Zealand, the South Island of New Zealand, and the main island of Tasmania (Figure 3i).
Non-manifold dataset. Our framework for identifying subspaces is still valid for cubical complexes built from images (by Corollary 2). We demonstrate the idea on RETINA, a medical retinal image [25]. The cubical complex is constructed by intensity thresholding (also called the sub-level set method in TDA [58]) and then applying morphological closing on the binary image to remove small cavities. The weight for every rectangle w2(σ) is set to 1; the estimated null space dimension of the L1 built from CB is β1 = 12, with the identified homologous loops in Figure 3j. The result shows the robustness of the proposed framework even for large β1.

heading: 7 Conclusion and broader impact
text: Our contributions in the emerging field of spectral algorithms for k-Laplacians Lk [12, 21, 35, 48] are summarized as follows. (i) We extend the study of the homology embedding of vertices by the graph Laplacian L0 (spectral clustering) to those of higher-order simplices by Lk. Specifically, the k-th homology embedding can be approximately factorized into parts, with each corresponding to a prime manifold given a small perturbation (small ϵℓ and ϵ′ℓ for ℓ = k, k−1). (ii) The analysis is made possible by expressing the κ-fold connected sum as a matrix perturbation. This convenient property of the homology embedding supports (iii) the use of ICA to identify each decoupled subspace and motivates (iv) the application to the shortest homologous loop detection problem. The proposed framework opens up numerous future directions for us to explore. For instance, one can extend our framework to investigate the generalization of the loosely connected clusters, i.e., when the null spaceHk eigenvalues are not strictly zero; it applies to the case when analyzing the noisy topological structures in the data manifold. Moreover, the connections of our framework to the recent advances in spectral TDA and representation learning can be further explored; they include the persistent spectral method [38, 57] and the disentanglement of representation [63] in generative modeling.
Broader impact. Our analysis provides insight into the structure of the k-th harmonic embedding. This framework can inspire researchers in developing spectral topological data analysis algorithms (e.g., visualization, clustering, tightest higher-order cycles for k ≥ 2 [20, 42]) similar to those that were inaugurated by spectral clustering two decades ago. These applications are especially beneficial to scientists (chemists, biologists, oceanographers, etc.) who use high-dimensional data analysis techniques for studying complex systems. Similar to the limitation of other unsupervised learning algorithms, practitioners without solid understandings of both the analyzed datasets and the used algorithm might draw controversial conclusions (see, e.g., discussions in [2, 41]). Possible approaches to mitigate the negative consequences are to design proper validation and causal inference algorithms for this framework; we leave them as potential directions we will explore.

heading: Acknowledgements
text: The authors acknowledge partial support from the U.S. Department of Energy’s Office of Energy Efficiency and Renewable Energy (EERE) under the Solar Energy Technologies Office Award Number DE-EE0008563 and from the National Science Foundation award DMS 2015272. They thank Alexandre Tkatchenko, Jim Pfaendtner, Stefan Chmiela, and Chris Fu for providing the molecular dynamics data as well as for many hours of brainstorming and advice.
Disclaimer
The views expressed herein do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Review summaries:
This paper presents a method for calculating (or rather *estimating*) the null space of the $k$th Laplacian. This enables the analysis of the geometry of $k$th homology embedding, making it possible to factorise such an embedding into different subspaces corresponding to individual prime factors of a given manifold. The proposed method also turns out to be applicable to improve loop homology generators, leading to geometrically more concise ('localised') representations of such loops.  The new method is showcases on a variety of data sets, demonstrating its overall utility for loop localisation and geometric decompositions. 

In this manuscript, the authors study the decomposition of the homology spaces from Hodge Laplacian matrices, and propose a new algorithm for the detection of the shortest homology generators. In their model, the direct sum of homology embedding or spaces from Hodge Laplacian is approximated by connected sum of manifolds. More specifically, the homology space is decomposed into subspaces with simple topological components. An algorithm is designed to detect shortest homology loops in the structure. The approach to use homology spaces for clustering and detection of shortest loops is very novel and interesting. It expands the area of topological data analysis, and will have a great impact for various data modelling and analysis. 

This paper studies the kernel of the higher-order combinatorial Laplacian, which equivalently amounts to a study of higher-order homology. In particular, given finite samples from a manifold with a ground truth decomposition into a connected sum of smaller ("prime") manifolds, they study the question of when it is possible to computationally recover a "nice" homology basis such that the basis vectors are resolved into independent components contributed from each of the prime manifolds. In this sense, the theoretical aim is quite nice and fundamental. The authors are able to provide a positive answer to their question in the setting where the prime manifolds are connected over a small, sparsely connected region. The proof technique is similar to that for obtaining bounds in spectral graph theory, but is carefully implemented because the problem has lots of block matrices floating around. Finally, the authors propose algorithmic applications of this method to the problem of detecting shortest loops in data, and exemplify their constructions on several point cloud and image datasets. 

This paper studies $k$-th homology vector spaces through the lens of $k$-Laplacian, a generalization of the usual graph Laplacian to higher order combinatorial structures.  It studies in particular the (fairly general) context of decomposable manifold as a sum of prime manifolds, showing (losely speaking) that under some assumptions (sparse connections), the homology space of a manifold can be obtained as the sum of the homology spaces of the corresponding prime manifolds.  Building on Independent Component Analysis technique, it then provides a way to identify the subspaces of the homology basis computed on the whole manifold to recover the basis of the corresponding prime manifolds.  This has applications in particular in homologous loop detection (detecting the shortest generating element in H1), as showcased on both synthetic and real datasets.