Joint Hierarchical Representation Learning of Samples and Features via Informed Tree-Wasserstein Distance
Authors: Ya-Wei Eileen Lin
Ronald R. Coifman
Gal Mishne
Ronen Talmon
Venue: NeurIPS 2025 spotlight ()
Keywords: hierarchical representation learning
unsupervised learning
tree-Wasserstein distance
diffusion geometry
hyperbolic embeddings
Haar wavelet
manifold learning
hyperbolic graph convolutional networks
Primary area: general_machine_learning
TL;DR: We introduce an iterative framework to jointly learn hierarchical representations of both samples and features using tree-Wasserstein distance and data-driven Haar wavelet filters.

Abstract: High-dimensional data often exhibit hierarchical structures in both modes: samples and features. Yet, most existing approaches for hierarchical representation learning consider only one mode at a time. 
In this work, we propose an unsupervised method for jointly learning hierarchical representations of samples and features via Tree-Wasserstein Distance (TWD). 
Our method alternates between the two data modes. It first constructs a tree for one mode, then computes a TWD for the other mode based on that tree, and finally uses the resulting TWD to build the second mode’s tree. By repeatedly alternating through these steps, the method gradually refines both trees and the corresponding TWDs, capturing meaningful hierarchical representations of the data.
We provide a theoretical analysis showing that our method converges.
We show that our method can be integrated into hyperbolic graph convolutional networks as a pre-processing technique, improving performance in link prediction and node classification tasks. 
In addition, our method outperforms baselines in sparse approximation and unsupervised Wasserstein distance learning tasks on word-document and single-cell RNA-sequencing datasets.

Joint Hierarchical Representation Learning of Samples and Features via Informed Tree-Wasserster Distance
Lei Eileen Lin' Ronald R. Coifmann 2 Gal Mishne
Ronen Talmon'
Viterbi Faculty Electrical and Computer Engineering, Technion
Department of Mathematics, Yale University
Halicioglu Data Science Institute, University of California San Diego
Abstract
High-dimensional data often exhibit hierarchical structures in both modes: samples and features. Yet, most existing approaches for hierarchical representation learning consider only one mode at time. In this work, we propose an unsupervised method for jointly learning hierarchical representations of samples and features via Tree- Wasserstein Distance (TWD). Our method alternates between the two data modes. It first constructs tree for one mode, then computes TWD for the other mode based on that tree, and finally uses the resulting TWD to build the second mode's 's tree. By repeatedly alternating through these steps, the method gradually refines both trees and the corresponding TWDs. capturing meaningful hierarchical representations of the data. We provide theoretical analysis showing that our method converges. We show that our method can be integrated into hyperbolic graph convolutional networks as a pre-processing technique, improving performance in link prediction and node classification tasks. In addition, our method outperforms baselines in sparse approximation and unsupervised Wasserstein distance learning tasks on word-document and single-cell RNA-sequencing datasets.
Introduction
High-dimensional data with hierarchical structures are prevalent numerous fields, e.g., gene expression [1-3], image analysis [4-6], neuroscience [7, 8), and citation networks [9-11]. Therefore, finding meaningful hierarchical representations for these data an important task that has attracted considerable attention [12-20]. While hierarchical structure is often assumed to exist in either the samples or the features [13, 21-23] (i.c., one of the data modes of matrix), many real-world datasets exhibit hicrarchical structure in both modes. For example, word-documeet data [24-26] commonly contain hierarchies within features (e.g., related keywords) and within samples (e.g.. document topies). Another example is recommendation systems [27, 28), where items (features) could be arranged into taxonomies or product categories, and users (samples) can be represented by multiple levels of 30].
behavioral or demographic relations [0).
An emerging approach to represent hierarchical data is based on embedding in hyperbolic space [13, 14, 21-23, 31J. However, most existing hyperbolic representation learning methods focus only on one mode of a data matrix [13, 21, 22]. either the rows (samples) or the columns (features). A straightforward way handle hicrarchies in both modes is to learn separate hierarchical represen- tation for each. We postulate that a joint approach facilitates significant advantages and investigate how the hierarchical structure of features conversely, how the sample hierarchy can improve the discovery the feature hierarchy.
improve the discovery of the sample hierarchy, and,
To this end, we propose to jointly learn the hierarchical representation of features and samples using Tree- Wasserstein Distance (TWD) [32] in an iterative manner. Our approach departs from
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Word
Dooners TAG PaRt
Aor TeD Mat
Inital TWD
Rofina Doc TWO Muste
Umii Converge
Dec. Hear Wavsiet
Word Hour Raveet
Docemen Toee
Winl Tre
Figure I: Overview of learning hierarchical representations across samples and features using TWD. Consider word-document data matrix. We construct initial tree for one data mode (e.g., words). This tree is then used to compute the TWD in the other mode (e.g.. documents). The newly computed TWD informs a tree update of that mode, and the updated tree is subsequently used to compute the TWD in the cross-mode. This alternating procedure continues iteratively, refining both trees.
standard uses of TWD, where trees are typically built for computational efficiency in approximating Wasserstein distances with Euclidean ground metrics [32-36]. Instead, we use trees learned from hierarchical data through diffusion geometry and hyperbolic embedding [37], allowing the hierarchical structure of one mode to be incorporated in the computation of a distance of the other. We start by learning an initial tree of one mode (either features or samples). Then, this learned tree is incorporated into the computation of the TWD of the other mode. The computed, informed TWD is then used to infer a tree of that mode, i.c., sample TWD is used to update the sample tree and feature TWD is used update the feature tree, shown in Fig. 1. This procedure is repeated, alternating between the two modes, iteratively refining both trees and recomputing the corresponding TWDs. We provide theoretical guarantees that this iterative process converges and that fixed point exists.
Based on the learned trees that represent the hierarchical structures of the data, we extend the iterative algorithm by incorporating additional adaptive filtering step. This filtering is implemented using Haar bases [38-41], which are constructed based on the learned trees. More concretely, at each iteration, we view each sample (resp. feature) as a signal supported on the feature tree (resp. sample tree) [42]. We then apply Haar wavelet filters [38, 43] induced by the feature tree (resp. sample tree) to the samples (resp. features), illustrated in Fig. 1. Assuming the trees reflect the intrinsic hicrarchical structure, applying the resulting data-driven wavelet filters can remove noise and other nuisance components [44]. We show that integrating the filters into the iterative process also leads to convergence. Empirically, we demonstrate that this approach improves hierarchical representations in terms of sparse approximation, and the resulting TWD leads to superior performance in document and single-cell classification.
We further demonstrate the practical benefit of our hierarchical representations by using them to initialize hyperbolic graph convolutional networks (HGCNs) [21, 45). By incorporating our data- preprocessing step, observe improved performance in
driven hierarchical representation link prediction and node classification tasks compared to standard HGCNs. This highlights both the compatibility and effectiveness of our approach for hyperbolic wraph-basce models.
Our contribution. (1) We present an iterative framework for jointly learning hierarchical represen- tations of features and samples using TWD. (2) We further enhance the hierarchical representations using Haar wavelet filters constructed from the learned trees. (3) We show that our method achieves superior performance in sparse approximation, as well as document and single-cell (scRNA-seq) classification. (4) We also show that the proposed hierarchical representation can be used to initialize HGCNs, improving link prediction and node classification on hierarchical graph data.
2 Related work
Tree- Wrre-Wssssrstinn distance. Tree- Wasserstein distance (TWD) [32-36] was introduced to mitigate the high computational cost of the Wasserstein distance [46], which is a powerful tool to compare sample distributions while taking into account feature relationship [47-49]. Most TWD methods involve two steps: (i) constructing feature tree and (ii) using this tree to compute sample TWD matrix. While the tree is typically built to approximate the Euclidean ground metric in TWD literature.

---

[37] recently proposed tree construction via hyperbolic embeddings [50] and diffusion geometry [5 I], enabling the tree metric to reflect geodesic distances on a latent unobserved tree underlying features.
Co-Manifold learning. Co-manifold learning aims to jointly recover the geometry of samples and features by treating their relationships as mutually informative, i.e., the feature manifold informs the sample manifold, and vice versa. It has been applied in joint embedding [52-55] and dimensionality reduction [7, 56, 57], assuming that both rows and columns of the data matrix lie on smooth, low- dimensional manifolds. While effective in capturing geometric structures, these methods often rely on Riemannian manifolds with non- non-negtivve curvature assumptions [55, 58], which may not be suited well with the negative curvature often associated with hierarchical data.
Unsupervised ground metric learning. Wasserstein singular vectors (WSV) [59] were introduced to jointly compute Wasserstein distances across samples and features, where the Wasserstein distance in one mode (e.g., samples) acts as the ground metric for the other (e.g., features). However, WSV is computationally expensive. To address this, Tree-WSV [60] was recently proposed, using tree-based approximations of the Wasserstein distance to reduce computational complexity. While these methods alternate between computing distances in manner similar to ours, neither WSV nor Tree-W WSV was designed to learn the hierarchical representations of the data. In contrast, our method explicitly learns the hierarchical representation of samples and features, further enhanced by wavelet filtering and its integration into HGCNs. These aspects were not explored in the WSV or Tree-Wo WSV frameworks.
3 Background
Tree construction using hyperbolic and diffusion geometry. Consider a set of high-dimensional points 2 = {zg E [R" " and let M E Rm be a suitable pairwise distance matrix between these points. In a recent work [37], a binary tree was derived from hyperbolic embedding obtained through multiscale diffusion densities built from M [50], employing the concept of lowest common ancestors [61] within this embedding space. The resulting tree, denoted T (M), has m leaves corresponding to the m points. Details on the embedding and tree construction are in App. A.
Tree- rres-assernttin distance. Consider tree T (V, E, A) with Nier leaves, where V is the vertex set with N nodes, E is the edge set, and A E N is the edge weight matrix. The tree distance dr is the sum of weights of the edges on the shortest path between any two nodes T. Let (v) be the set of nodes in the subtree rooted at v E V. Each u E V has a unique parent v, with edge weight w du v) [62]. Given distributions Pl: P2 E R giai supported on T, the TWD [32] is defined by
TW(P1, P2, T) = Ewev we EmETroo) (pi(u) - pp(())
(1)
Haar wavelet. Consider a complete binary tree B with m leaves. Let E = -et 1,..., L denote the levels in the tree, where f 1 the root level and é L is the level. Let T(e, s) be the set of all leaves in the subtree, whose root is the s-th node of the tree at level f, where = 1.....N, Nr and N2 is the number of nodes in the ê-level. At level f, subtree T(e, splits into two sub-subtrees f(l + 1, $1) and T(E 1, s2). A zero-mean Haar wavelet Br.- E R has non-zero values only the indices corresponding to leaves in the sub-subtrees and is piecewise constant on each of them (see App. for an illustration and further details) [40, 39, 433. The set of these Haar wavelets, along with a constant vector, is complete and forms an orthonormal Haar basis, denoted by B E Rm with each column corresponding a Haar wavelet (basis vector). Any vector a RT can be expanded in this Haar basis as a = Σe, Be.e, where ag = (a, Br,,) is the expansion coefficient.
4 Proposed method
Problem setting. Given data matrix X E I with n rows (samples) and m columns (features), we denote X,, and X., as the i-th sample and the j-th feature, respectively. Our goal is to learn hierarchical representations for both samples and features. We model the hierarchical structure of the data by constructing rooted weighted trees as follows: a sample tree T, with n leaves, where each leaf represents one sample, and a feature tree T, with leaves, where each leaf represents one feature.
3

---

Algorithm Iterative Joint Hierarchical Representation Learning via Free-Wasserstein Distance
Input: Data matrix X E 17 Me E R T and M, E _=*" Tr. Ve >
Output: Trees T(W?)) and T(W?). and TWDs w
()) and w!
I + 0, X, fr, X and X. (c, X /
Initialization
: (oX,;T((M))
D Construct initial feature tree and compute initial sample TWD
( - ((x:: (MM,))
P Construct initial sample tree and compute initial feature TWD
repeat
w f- + T(W!)) and w + (x:::T((.))
Iterative update
/-1+1 until convergence
Placing data points as leaves follows the established line of works in statisties [63-65], manifold learning [52-55, 7, 43], and TWD literature [32-36]. We propose an iterative scheme that alternates between learning the sample tree and learning the feature tree.
4.1 Iterative scheme for joint hierarchical representation learning
To jointly learn hierarchical representations for samples and features, we use TWD as a means to facilitate the relationships between samples and features. We begin by constructing tree for one data mode; without loss of generality, we first construct an initial feature tree. Given pairwise distance matrix M. E mg over the m features, we build complete binary tree T(M.) with m leaves [37], where each leaf corresponds to a feature (see App. for details). This feature tree then serves as the hierarchical ground metric for computing TWD between the samples:
W (0) (i, = TW. rv. T(M.)) + 7-S(r; -r))
where r; = XT | is the i-th normalized sample, is norm regularize based on the snowflake 1,37
penalty [57], and 7/r O. For brevity, we write the resulting sample pairwise TWD matrix
( (X,; ((M))) R=x""
(3)
where 8, = [r....., rn] E Raxm and & denotes a function that computes the pairwise TWD
between rows of the matrix X, using the tree T TMM.) defined over features. A similar analogous construction applied the other mode (the samples). Given an initial pairwise distance matrix M, R=x= over the n samples, we construct an initial sample tree T(M,) with n leaves and use it compute TWD between the features:
w(O) (j.j') TW cy, T(M,) + Tel(e, c.y).
(4)
where cj = X..// || and 7e > 0. The resulting feature pairwise TWD is then written as
( (8x: T(M,)) gm
(5)
where X, -c,,,.... Cm mxxx These matrices of the initial sample and feature pairwise TWDs in Eq. (3) and Eq. (5) serve as the starting point of the proposed iterative scheme.
Our iterative scheme alternates between the two data modes in coordinate descent manner [66, 67]. At iteration = 0, the initial pairwise distance matrices M, and M. are used to construct the corresponding sample and feature trees, respectively. These trees are then used compute the initial TWDs in Eq. (3) and Eq. (5). For all subsequent iterations / > 1, the learned feature TWD from the previous step is used to construct new feature tree, which is then used compute the updated sample TWD. The same process is applied to the other mode: the sample TWD is used to construct new sample tree, which turn used compute the subsequent seaseqquetfeaturre TWD. Formally, using the notation introduced above, the update steps at + iteration are given by
W = (x,,; ((W!!)))))))))))) ()) B=*a
W = 4 T(W()) Rm
(6)
This alternating scheme, outlined Alg. allows the hierarchical representation in one mode to iteratively inform and refine the representation in the other mode.

---

The sequences w and : generated by Alg. I have at least one limit point, and all limit points are fixed points if Tr To > 0.
Theorem 1.
The proof is in App. C. Thm. implies the existence of a limit point to which the proposed iterative scheme converges. Alg. 1 can be implemented in practice without regularization, i.e., 7 = To 0. However, the conditions Yr, Te > 0 are necessary for Thm. 1.
We remark that while other TWD methods [32-36] could in principle be incorporated into iterative framework, we chose the tree construction method [37] for two main reasons. First, we empirically observe that computing sample and feature TWDs with these alternative TWD methods often fails to converge. Second, their use as cross-mode tree references would yield trees whose tree distances approximate the Wasserstein distance. However, the Wasserstein metric is not inherently a hierarchical metric. As a result, the trees derived from such approximations do not represent the hierarchy present in the data. In contrast, constructing hierarchical representation using [37] yields tree whose geodesic (shortest path) distances reflect the hierarchical relationship underlying the data [68, 51]. As shown empirically in Sec. 6, the trees obtained by Alg. 1 yield meaningful hierarchical representations of features and samples across benchmarks from multiple domains, leading to improved performance compared using other TWD methods within our iterative scheme. In addition, we show that the trees and the corresponding TWDs are progressively refined throughout the iterations.
4.2 Haar wavelet filtering
To improve the refinement of trees across iterations, we apply filtering step at every iteration. Specifically, our filters are constructed from Haar wavelets [38-41]. was shown that Haar wavelets can be derived adaptively from trees [43]. Here, we propose to build the Haar wavelets from the trees inferred through the itcrative process. Viewing each sample fees.ffaturee as a signal supported on the feature tree (resp. sample tree) [42] allows us apply data-driven Haar wavelet filters. Since by construction the filters reflect the intrinsic hicrarchical structure of the data, they enhance this meaningful representation while suppressing noise and other nuisance components.
Given a feature tree T(Me), we construct a Haar basis B. E Rm associated with the tree (see Sec. 3 and App. A for details). Subsequently, each sample can be expanded in this Haar basis, and we denote O;; = B.) E RRm as the vector of the expansion coefficients of the i-th sample. To define wavelet filter, we select a subset of the Haar basis vectors in B, as follows. For each coefficient index j. we compute the aggregate I,1 norm Ex= lox, and sort these values in descending order. Then, we sequentially add the corresponding indices to Q until the cumulative contribution 1/2 = Eyee Ei (4) exceeds threshold 0. > O. Let B Rmxd denote the matrix consisting of the d basis vectors in Q. The filtering step. which yields the filtered samples, is defined
T (X:; ((M.)) (X E mm
where Vi denotes wavelet filtering operator applied to the rows of the matrix X, using the Haar wavelet induced by the tree T(M,) defined over features. An analogous wavelet filter can be constructed from sample tree T(M,)) and applied to the features:
(Z (M,)) = (ZB,),, E RR
(8)
where Z- = XT, and B, E Rn consists of the top d' basis vectors selected using threshold 8, > 0 on cumulative coefficient magnitude.
We can apply this Tree Haar Wavelet filtering in our joint iterative scheme to update the trees as follows. Given initial inputs ( = X, Z = T and using the notation introduced above, at each iteration 2 O, filter the data using the Haar basis vectors derived from the trees:
x (I+)) ( (1) ((W.))) E m
(1+1) v ( " ((W(")) E R
(9)
where W (0) M, and w!! M. for iteration I = 0. The filtered data is normalized into discrete
histograms! Then we refine the trees and TWDs based on the normalized Haar-filtered data:
W(+1) क X (1+1), ((((()))
E Rx
W([+1]
क (+1)) (T((())) E Rmm
(10)
where X (1+1) and X (++1) are the column-normalized matrices of X ((+1) and (++)), respectively.
We summarize this iterative learning scheme with the Haar wavelet filters in Alg. 2.
The resulting filtered data may contain negative values. To represent histogram, we subtract the vector with minimum value.
5

---

Algorithm 2 Haar Filtering over Alternating Tree Refinement
Input:
Data matrix X 17 MM
R Raxm and M, R"*", 7/ and Ye: thresholds de and 0,
Output: Trees T(W," and T(W?). and TWDs W
()) and W
t O, (
X. Z (0) XT,
M, and ! M.
Initialization
repeat
X(1-1)
(x T( ((W )) and z V (T(.))
Tree haar wavelet filtering
X
(+1))
(x /
x(111)
(1+1) (z 1 l1
Cj
w
t 1 (X, 11+ :(W!)) and w
1 (Xx (1+1),
((W!)))
D Iterative update
1+1+1 until convergence
Theorem 2.
The sequences W and w generated by Alg. 2 have at least one limit point, and all limit points are fixed points if yr To > 0.
The proof of Thm. 2 is in App. C. Same as Thm. 1, Thm. 2 implies the existence of a limit point to which the proposed iterative scheme with the addition the Haar wavelet filters converges. We argue that the resulting data-driven wavelet filters attenuate noise and other nuisance components [43], and therefore improve the quality of the learned trees. As shown in Sec. 6, this filtering step contributes to more informative hierarchical representations underlying the digensinesio data, resulting in superior performance on various tasks. We note that while constructing the filters using L1 based selection criterion [69] is effective, alternative filtering strategies [70-75] can be considered, depending on the downstream tasks. We leave this extension for future work.
We conclude this section with few remarks. First, while our iterative procedure is broadly applicable to various TWD methods [32-36], the theoretical results and the ability to obtain meaningful hierar- chical representations rely our choice of the specific tree construction method [37]. Second, in each iteration, our method can be computed in O + m.2 [37, 76] (see App. E). contrast, a naïve computation requires (mn"² + m² log m nmn + n log n),, making our method more efficient. Third, similar to WSV [59] and Tree-WST ground metric learning technique [77]. However, its specific focus learning hicrarchical repre- sentations for both samples and features distinguishes it from these methods and leads to superior empirical performance for hierarchical data (see Sec. 6). For further comparison with Tree- Tree-WSV, refer to App. B and App. G.3. Finally, we remark that either the sample or the feature hierarchical structure can be provided as a prior and used for initialization of our method (see Sec. 5).
[60], our approach can be viewed as an unsupervised
5 Incorporating learned hierarchical representation within HGCNs
In Sec. 4, we introduced joint hierarchical representation learning framework using TWD. Typically, the initial pairwise distance matrices M, and M. are data-driven, e.g., using standard metrics such as the Euclidean distance or the cosine similarity. An advantage of our method is that it can incorporate prior knowledge when a hierarchical structure available for one of the modes. Without loss of generality, we consider scenarios where hierarchical structure over the samples is known and represented by a graph H ([n), E, A), with [n] = {l.....,., n... To integrate this prior, we initialize M, using the shortest-path distances dæ induced by H. This initialization introduces structured prior in one mode while allowing the hierarchy in the other mode to be learned jointly via TWD.
Incorporating such prior is particularly relevant for hierarchical graph data [9, 21, 78], where the node hierarchy is provided, while the feature structure remains implicit. We demonstrate the compatibility of our approach with hyperbolic graph convolutional networks (HGCNs) [21, 45]. Specifically, during the neighborhood aggregation step in HGCNs, we replace the predefined hierar- chical graph with the sample tree inferred by our method after convergence. This learned tree guides the aggregation process, enabling hierarchy-aware message passing that reflects structured relations among samples and across features. As shown in Sec. 6, incorporating method pro-processing

---

Table 1: The L1 norm of the Haar expansion coefficients. Values reported in the format samples / features (the lowest bold and the second lowest underlined).
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Co-Quadtree
26.6/
24,
64.3/ 108.7
87.3/ 102.7
157.0/ 242.4
1616.6/ 77.3
Co-tteeeeee
27./31.5
69.1 20.9
73.7/97.5
88.1/110.1
173.7/ 240.0
1688.8/ 81.7
Co-WCTWD
2
57.6/ 17.1
63.3/817
77.4/664
136.4/202.9
1308.1 68.44 68.4
Co-WQIWD
25.2/32.1
56.9/
61.3/100.1
11.7/111.0
135.0/ 224.8
1324.8/674
Co-UltraTree
37.6/32.1
69.8/ 18.8
76.0/
86.3/133.3
155.4/ 226.6
1450.0/82.2
Co-TSWD-1
26.0/
70.1 15.0
69.8/91.6
100.4/11..2
/
1716.1/700
Co-TSWD-5
30.5
58.9 16.0
65.7/98.0
83.3/102.5
170.2 234.7
/ 70.7
Co-TSWD-10
21.1/34.5
53.6
14.3/
7
141.4 229.7
1373.3 /
Co-SWCTWD
29.5
/9.0
/9.6
95.9/104 104.6
2
1.55.0
89.7
Co-SWOTWD
33.5/
57.6/ 13,9
61.3 86.0
77.0/98.9
14181 209.6
1369.8 68.0
Co-MST-TWD
30.4/4.1
69.1
77.5 109.3
97.0/ 126.4
2
1755.3 83.7
CO-TR-TWD
/4.8
/9.1
/9.1
82.2 / 102.3
124.6/ 238.6
926.8/ 68.8
Co-HHC-TWD
22.5/22.7
51.51 13.6
/93./
74.1 110.1
/
1148.8 79.0
Co-gHHC-TWW
20.9/10.6
65.0/ 16.7
72.8 112.9
88.5 115.5
162.71 240.7
1612.8 70.8
Co- TWD
22.5/22.1
51.9/ 13.4
88.4/8133
73.1 /92.8
133.7/ 202.0
1294.6 6 78.1
QUE
19.8/198
$7.8 10.4
77.6/
88.8 72.0
93.1/ 173.3
906.4 58.5
Tree-WSV
23.6/24.9
54.3 18.2
65.4
84.2 106.3
/
1637.4 73.2
Alg.
12.9/7.0
25.4/3.7
40.4 24.6
51.0/ 30.1
64.4/ 110.8
511.3/ 50.1
Alg. 2
10.1/4.
22,4/3.4
37,2/ 19.4
46.6 26.6
50.9 / 93.7
489.4/ /5.6
step improves performance on link prediction (LP) and node classification (NC) for hierarchical graph datasets. While our focus here on HGCNs due to their improved performance NC and LP, the same initialization strategy can be adopted in other hyperbolic architectures, such as hyperbolic neural networks (HNNs) [79, 80]. Further technical details on this integration are provided in App. D.
Experimental results
We evaluate our methods on sparse approximation and unsupervised Wasserstein distance learning with word-document and scRNA-seq benchmarks. Additionally, we examine the integration of our methods into HGCNs for hierarchical graph data on LP and NC tasks. The implementation details. including hyperparameters, are reported in App. E. Additional experiments, e.g., empirical convergence, ablation study, runtime analysis, and co-clustering performance, are presented in App. F.
6.1 Evaluating hierarchical representations via sparse approximation
We first demonstrate the advantages of the learned trees from Alg. 1 and Alg. for sparse approx- imation tasks on the data matrix. The quality of the feature tree (and similarly, the sample tree) is evaluated by the L1 norm of their expansion coefficients across all samples (and the fcatures, respectively) [81-83]. A lower L1 norm indicates a more efficient (sparser) representation of the data using fewer significant Haar coefficients, thus indicating the learned tree structures better reflect the hierarchical information of the data [43, 53]. We test four word-document datasets
[25]: BBCSPORT, TWITTER, CLASSIC, and AMAZON, and two scRNA-seq datasets [84]: ZEISEL and CBMC. Both types of data exhibit hierarchical structures in their features and samples [85, 86]. In document data, words (features) form semantic hierarchies (e.g., animal - mammal - dog), while documents (samples) follow topic-subtopic structures (e.g., science - biology - genomies). In scRNA-seq data, genes (features) are organized by functional relationships such as biological pathways and gene ontologies (c.g., immune response genes - cytokine genes - specific interlcukins), and cells (samples) follow developmental or taxonomic hicrarchies (c.g., hematopoictic stem cell - progenitor cell + mature blood cell types), which are widely modeled as hierarchical in computational biology [87]. Additional information about these datasets can be found in App. E.
We propose to learn hierarchical representations for both samples and features, where each informs the other through TWD. To the best of our knowledge, this approach to hierarchical representation learning has not been previously explored. demonstrate its effectiveness, we compare Alg. and Alg. 2 against existing TWD- bWDD- methods as follows. For a fair comparison, we adapt each baseline TWD method our iterative setting. Specifically, we begin by using the competing method to compute the sample TWD matrix, from which the sample tree is constructed to approximate the corresponding
7

---

Table 2: Document and single-cell classification accuracy.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Co-Quadtree
96.2+0.4
69.6+0.3
95.9=0.2
89,4+0.2
81.7=100
80.740.3
Co-Flowtrce
95.7+0.9
71.5+0.77
95.6=0.5
91.4+0.4
84.3=0.7
83.0+1.2
Co-WCTWD
93.2+1.2
70.242.11
94.7-2.6
87.4110
82.5-2.9
79.412.1
Co-WOTWD
1.
70.7+22.2
9..
88.2+2..1
82.3+3.11
80.5
Co-UltraTree
9..
70.1
93.6=2.0
86.5+288
85.8=1.1
84.6+1.3
Co-TSWD-1
88.211.4
70.411.2
94.710.9
86.110.55
80.2-1.4
73.2.11.0
Co-TSWD--S
88.7+1.7
71.0+1.5
96.7+0.88
91.5+0.4
82.0-0.9
75.4+0.7
Co-TSWD-10
89.2+1.11
71.4188
95.5=0.2
91.820.7
83.8=0.5
77.2+0..9
Co-SWCTWD
93.512.4
70.5.1.0
9.44113
90.7.11.5
82.7.1.7
79.0.10.9
Co-SWQIWD
96.2111.2
72.412.1
96.0=1.1
90.612233
82.4-1.4
813111
Co-MST-TWD
88.7+-2.4
68.4+3.3
91.3++2.9
87.1+1.4
80.1-2.8
76.5+
Co-TR-TWD
89.5+112
70.9+1.7
93,4-2.2
89.5+1.4
80.7=0.8
78.5+0.9
Co-HIH-TTW
86.112..1
70.1.11.3
93.611.5
88.510.5
83.2-1.4
77.610.8
Co-gHHC-TWD
84.0++200
70.4+1.6
90.7+1.7
87,2+1.9
79.9-1.4
84.2+1.2
Co-UItaFit-TW
86.8+0.9
70.9
91.9=1.0
89.942.0
83,7=2.9
79.1
QUE
84.7
72.4.00.6
91.9
91.6110.9
83.6
82,5.199
WSV
85.9+1.0
714+1.
92.6=0.7
89.0+1.5
81.6-244
17.
Tree- WSV
16,3+
712+1.9
92,4+1.0
88.7+1.9
82.0-2.9
76.4+2.4
Alg.
96.7....3
74.1
97.310.2
94.0.10.4
0.04
86.7.10.5
Alg.
97.3+0055
76.7+00.7
97.6=0.1
94,2+0022
94.0=0.6
93,3++0.7
tree-based Wasserstein ground metric. This sample tree is then used to compute the feature TWD matrix, which in turn defines the feature tree. each iteration, the same baseline method used to compute the TWDs for samples and features. The competing TWD methods include: Quadtree [32], Flowtree [88], TSWD [35], UltraTree [36], weighted cluster TWD (WCTWD), weighted
Quadtree TWD (WQTWD) [34], their sliced variants SWCTWD and SWQTWD, MST [89]. Tree Representation (TR) [22), gradient-based hierarchical clustering (HC) in hyperbolic space (gHHC) [90], gradient-based Ultrametric Fitting (UltraFit) [91], and HC by hyperbolic Dasgupta's cost (HHC) [92]. See App. for details on these methods. We use the prefix "co-" " to denote the adaptation cach method to iterative framework. In addition, include comparisons with co-manifold learning that involves trees induced by a diffusion embedding (QUE) [52], and Tree-WSV [60], which learns unsupervised ground metries based tree approximation of WSV [59].
Tab. reports the L, norm of the Haar coeffi- cients across all samples and all features, respec- tively. We report the value after convergence for our methods, and for baselines, we either report it after convergence or, convergence is not achieved, the obtained value after 25 itera- tions. We see that our methods provide a more efficient representation, showing that the spar- sity and quality of the trees produced by our methods are superior and outperform baselines. Fig. 2 shows the L1 norm of the Haar expansion coefficients obtained by the proposed methods across iterations on ZEİSEL dataset, where we observe that the L1 norm is iteratively reduced and reaches convergence. Note that this is not the objective we are minimizing but conse- quence of our method of learning well hierarchical representations of the data. Notably, Alg. consistently achieves better sparse approximation performance than Alg. We attribute this provement to the wavelet filtering step, which jointly considers the data and the structure and further improves the quality of the learned hierarchies each iteration. One might ask whether wavelet filters could similarly benefit other TWD- baselines. In App. F, we test this by applying wavelets using trees constructed from various TWD baselines. The results show that it does not consistently improve the quality of the tree representations. We argue that this is because the trees in these methods are primarily designed to approximate the Wasserstein distance as the ground metric, rather than to represent the hierarchical structure of the data. Therefore, these trees not faithful hierarchical
Zeisel Sample Tree
Zeisel Feature Tree
160
Aig
200
Aig.1
140
120
1.40
120
60
5 10 15
5
1.0
15
20
of iteration
of iteration
Figure 2: The L1 norm of the Haar coefficients from the sample tree (left) and the feature tree (right) during the sparse approximation task across iterations on the ZEISEL dataset.
8

---

Table 3: ROC AUC for LP, F1 score for the DISEASE dataset, and accuracy for the AIRPORT,
PUBMED, and CORA datasets for NC tasks (the highest in bold and the second highest underlined).
DISEASE
AIRPORT
PUBMED
CORA
LP
NC
LP
NC
LP
NC
NC
EUC
59,8+2.
32,5
92.0+0.0
60.913344
83.3=0.1
4872=0.7
82.5+0.3
23.8
HYP
63.540.6
45,5+3.3
94.5+0.0
70.2+0.1
87.5=0.1
68.5=0.3
87.6=0.2
22.0+11.5
EUC-MIXED
49.6+1.1
35.2+3.4
91.5+0.1
68.3+2.3
86.0=1.3
63.0-0.3
84.4+0.2
46.1+0.44
HYP-MIXED
55.1
56.911.5
93.340.0
69.610.11
83.8-033
73.9=0.2
85.6.0.5
45.9100.3
MLP
72.610.6
28.812..
89.810.5
68.6100.6
84.1
72,4=022
83.1
$1.5110
HNNs
75.110.3
41.011.8
90.810.2
80.51
94.910.1
69.8-0.4
89.0.1.11
54.6110.4
GCN
64.7+0.5
69.7+0.44
89.3+0.4
81,4+0.6
91.1+0.5
78. 78.1-0.2
90.4+0.2
81.3+0.3
GAT
69.8+0.3
70.440.4
90.5+0.3
81.5+0.3
91.2-0.1
79.0=0.3
93.7=0.1
83.0+:0.7
GRAPUSAGE
65,9+0.3
69.1+0.6
90.4+0.55
82.1+0.5
86.2=100
77.4-2.2
85.5-0.6
77.9+2.44
SGC
65.11002
69.51002
89.810.3
80.610.1
94.110.0
78.9-0.0
91.5110.1
81.010.1
HGCNs
90.8+0.3
74,5+0.9
96.4+..1
90.6+0.2
96.3+0.0
80.3-0.3
92.9++011
79.9+0,2
H2H-GCN
97.0++.33
88.6+1.7
96.4+..1
0..3+
9.9
0.59=0.5
95.0++0.0
82.8+0.4
HGCN-AIE
93.210.66
87.910.7
93.7110.2
89.910.4
94.110.7
81.7=0.2
93.1=0.11
82.910.3
HGCN-AIg. 2
98.4100.4
89.410.3
97.210.1
92.110033
97.2=0.2
83.6-0.4
96.910.3
83,910.2
representations of the data, and thus, applying wavelet filters that depend on both data and tree structure fails to improve hierarchical representation learning.
6.2 Document and single-cell classification using learned TWD
We further demonstrate the effectiveness of the TWDs obtained by our methods through document and cell classification tasks. We compare our results with the same competing methods used in Sec. 6.1, and additionally include WSV [59] that learns unsupervised ground metrics as a baseline. Classification is performed using &NN based on the obtained distances, with voss-valio over five trials. Each trial randomly splits the dataset into 70% training and 30% testing sets.
Tab. shows the document and single-cell classification accuracy. The accuracy of our methods is based on the TWDs after convergence. For the baselines, the accuracy is reported either based on the distance obtained after convergence or the distance after 25 iterations if convergence is not achieved. We see that our methods outperform the baselines by large margin. This indicates that the TWDs, learned through our iterative scheme, effectively capture the interplay of the hierarchical structures between rows and columns. In addition, we observe that the classification accuracy improves with each iteration our methods (see Fig. 9 in App. F). While the competing methods also show marginal improvement through the iterative procedure, our methods consistently achieve better performance. Our approaches exhibit fast convergence, typically within 10- 14 iterations in all the tested datasets.
6.3 Link prediction and node classification for hierarchical graph data
Finally, we show the utility of our methods as a pre- pre-processing step for HGCNs [21, 25], evaluated on LP and NC tasks. We adhere the experimental setups and baselines used in these works maintain consistency. For LP task, we use Fermi-Dirac decoder [93, 13] to compute probability scores for edges. Then, the networks performance of I.P is assessed by measuring the area under the ROC curve (AUC). For NC task, we employ a centroid-based classification method [94], where softmax classifiers and cross-entropy loss functions are utilized. Additionally, an LP regularization objective is integrated into the NC task [21, 45J. The NC task evaluated using the Fl score for binary-class datasets and accuracy for multi-class datasets. We test four datasets [9, 95, 21], including CORA, PUBMED, DISEASE and AIRPOT. Descriptions of these datasets and their splits are included in App. E. We compare our methods with two shallow methods: Euclidean embedding (EUC) and Poincaré embedding (HYP) [13]. We also include comparisons with the concatenation shallow embeddings and node features, denoted as EUC-MIXED and HYP-MIXED. Furthermore, we include multi-layer perceptron (MLP) and its hyperbolic extension, HNNs [79], as well as four GNNs: GCN [96], GAT [97], GRAPHSAGE [98], and SGC [99]. Lastly, we include HGCNs [21] and H2H-GCN [45].
trained by minimizing cross-entropy loss with negative sampling. The

---

Tab. 3 shows the performance of integrating our methods with HGCNs on the LP and NC tasks compared the competing methods. We repeat the random split process 10 times and report the average performance and standard deviation. Our methods consistently outperform the competing baselines across both tasks. Similar to Tab. and Tab. 2, we observe in Tab. that using wavelet filters demonstrates superior performance by significant margin. This indicates that it effectively represents the hierarchical structure of the graph data, improving the expressiveness of both GNNs and hyperbolic embeddings. A natural question is whether wavelet filters alone (i.e., without our iterative scheme) could benefit HGCNs, We investigate it in App. F and find that such integration does not yield comparable outcomes. To our knowledge, with TWD and apply them to HGCNs. We note that our methods could have been integrated within HGCNs in other ways. However, shown in Tab. the straightforward integration already delivers favorable results. Thus, we opt to explore more complex integration techniques for future work.
work is the first to incorporate wavelets
Conclusions
This work introduces an iterative framework for jointly learning hierarchical representations of both samples and features using TWD. The proposed method begins by constructing a tree for one mode (either features or samples), which is then used to compute TWD and infer the tree construction for the other mode. The process alternates between modes, with each tree informing the inference of the other through pairwise TWD computations. To further improve the quality of the tree representations, we apply wavelet filters derived from the learned trees to the data at each iteration, which effectively suppress noise and filter out nuisance components. We show theoretically that the procedure converges and empirically that the trees and TWDs are refined across the iterations. Specifically, empirical evaluations on word-document and scRNA-seq datasets show that the resulting tree representations and TWDs lead to meaningful hierarchical representations. We further demonstrate that the proposed method can serve as a preprocessing step for HGCNs applied to hierarchical graph data, improving performance in hierarchical graph-based learning problems on link prediction and node classification.
Limitations and future work. One limitation of our approach lies in the use of trees to represent data geometry. While this representation aligns with our assumption that the data exhibits underlying hierarchical structures, it may not generalize well data supported on more complex geometries, e.g., spherical manifolds [100], spaces with mixed curvatures [101], asymmetric data [102-104], or general graphs [105]. In future work, we plan to explore more flexible geometric representation methods that can accommodate a broader class of data geometries, while utilizing the proposed iterative procedure between samples and features. We also plan to explore a differentiable variant of TWD, such as the soft TWD [106]. which could enable integration of our iterative process into neural architectures. We will further incorporate supervision or tasc-specific signals into the learning process. Finally, we plan to extend the approach to multi-way data (e.g., tensor-valued inputs).
Acknowledgments
We thank the anonymous reviewers for their insightful feedback. We thank Gal Maman for proofread- ing this manuscript. The work of YEL and RT was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 802735-ERC-DIFFOO GM was supported by NSF CCF-2217058 and CCF-2403452. The work of RC was supported by AFOSR MURI 052721 and DARPA-PA-23-03.
The work of
References
[1] Amos Tanay and Aviv Regev. Scaling single-cell genomics from phenomenology to mechanism. Nature, 541 (7637)3))3331 2017.
[2] Riccardo Bellazzi, Andrea Codegoni, Stefano Gualandi, Giovanna Nicora, and Eleonora Vercesi. The gene mover's distance: Single-cell similarity via optimal transport. arXiv preprint arXiv:2102.01216 2021.
[3] Malte Luecken and Fabian J Theis. Current best practices in single-cell RNA-seq analysis: tutorial. Molecular systems biology. 15(6))88746, 2019.
10

---

[4] Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian. Image denoising by sparse 3-d transform-domatt domain collaborative filtering. IEEE Transactions on image processing. 1668):2080-2095, 2007.
[5] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Deep image prior. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 9446-9454, 2018.
[6] Pascal Mettes, Mina Ghadimi Atigh, Martin Keller-Ressel, Jeffrey Gu, and Serena Yeung. Hyperbolic deep learning in computer vision: survey. International Journal of Computer Vision, pages 1-25, 2024.
[7] Gal Mishne, Ronen Talmon, Ron Meir, Jackie Schiller, Maria Lavzin, Uri Dubin, and Ronald R Coifman. Hierarchical discovery. IEEE Journal of Selected Topics in Signal Processing, 10(7):1238-13 1253, 2016.
analysis for neuronal structure and activity pattern
[8] Siyuan Gao, Gal Mishne, and Dustin Scheinost. Poincaré embedding reveals edge-based functional networks of the brain. In International Conference Medical Image Computing and Computer-Assister Intervention, pages 448 457. Springer, 2020.
[9] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi- Rad. Collective classification in network data. Al magazine, 29(3):93 93, 2008
[10] Elizabeth A Leicht and Mark EJ Newman. Community structure directed networks. Physical review letters, 100(11):11 188703, 2008
[11] Aaron Clauset, Cristopher Moore, and Mark EJ Newman. Hierarchical structure and the prediction of missing links in networks. Nature, 453(7191):98-001 2008.
[12] Shuchen Wu, Noćmi Élteto, Ishita Dasgupta, and Eric Schulz. Learning structure from the ground up-hherarchical representation learning by chunking. Advances in Neural Information Processing Systems, 35:36706 36721,, 2022.
[13] Maximillian Nickel and Douwe Kiela. Poincaré embeddings for learning hierarchical repre- sentations. Advances in Neural Information Processing Systems, 30, 2017.
[14] Maximillian Nickel and Douwe Kiela, Learning continuous hierarchies in the Lorentz model of hyperholic geometry. In International Conference Machine Learning, pages 3779-3788. PMLR, 2018
[15] Karan Desai, Maximilian Nickel, Tanmay Rajpurohit, Justin Johnson, and Shanmukha Ra- makrishna Vedantam. Hyperbolic image-text representations. In International Conference on Machine Learning, pages 7694-7731. PMLR, 2023.
[16) Aleksandr Ermolov, Leyla Mirvakhabova, Valentin Khrulkov, Nicu Sebe, and Ivan Oseledets. Hyperbolic vision transformers: Combining improvements in metric learning. Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 7409-7419, 2022.
[17] Wei Peng. Tuomas Varanka, Abdelrahman Mostafa, Henglin Shi, and Guoying Zhao. Hyper- bolic deep neural networks: A survey. IEEE Transactions on pattern analysis and machine intelligence, 14112):10033-10044 2021.
[18] Cheng Tan, Yijie Zhang, Zhangyang Gao, Bozhen Hu, Siyuan Li, Zicheng Liu, and Stan Z Li. Rdesign: hierarchical data-cfficient representation learning for tertiary structure-basee rma design. arXiv preprint arXiv:2301. 10774, 2023.
[19] Haochen Chen, Bryan Perozzi, Yifan Hu, and Steven Skiena. Harp: Hierarchical representation learning for networks. In Proceedings of the AAAI conference artificial intelligence, volume 32, 2018
[20] Gal Mishne, Zhengchao Wan, Yusu Wang, and Sheng Yang. The numerical stability of hyperbolic representation learning. International Conference on Machine Learning, pages 22949. PMLR, 2023.

---

[21] Ines Chami, Zhitao Ying, Christopher Ré, and Jure Leskovec. Hyperbolic graph convolutional neural networks. Advances in Neural Information Processing Systems, 32, 2019.
[22] Rishi Sonthalia and Anna Gilbert. Tree! I am no Tree! I am a low dimensional hyperbolic embedding. Advances in Neural Information Processing Systems, 33:845 856, 2020.
[23] Menglin Yang. Min Zhou, Rex Ying, Yankai Chen, and Irwin King. Hyperbolic representation learning: Revisiting and advancing. In International Conference on Machine Learning, pages 39639-39659, PMLR, 2023.
[24] Jonathan Chang and David M Blei. Hierarchical relational models for document networks. 2010.
[25] Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger. From word embeddings to document distances. In International Conference on Machine Learning, pages 957-966 PMLR, 2015.
[26] Gao Huang, Chuan Guo, Matt J Kusner, Yu Sun, Fei Sha, and Kilian Weinberger. Supervised word mover's distance. Advances in Neural Information Processing Systems, 29, 2016.
[27] Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recom- mender systems. Computer, 42(8)330-37, 2009.
[28] Francesco Ricci, Lior Rokach, and Bracha Shapira. Introduction to recommender systems handbook. In Recommender systems handbook, pages 1-35. Springer, 2010.
[29] Jizhe Wang. Pipei Huang. Huan Zhao, Zhibo Zhang, Binqiang Zhao, and Dik Lun Lee. Billion- scale commodity embedding for e-commerce recommendation in alibaba. In Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pages 839-848, 2018.
[30] Yulong Gu, Zhuoye Ding, Shuaiqiang Wang, and Dawei Yin. Hierarchical user profiling for e-commerce recommender systems. In Proceedings of the 13th international conference on web search and data mining, pages 223-231, 2020.
[31] Simon Weher, Bar Zöngür, Nikita Araslanov, and Daniel Cremers. Flattening the parent bias: Hierarchical semantic segmentation in the Poincaré ball. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 28223-28232, 2024.
[32] Piotr Indyk and Nitin Thaper. Fast image retrieval via embeddings. 3rd international workshop on statistical and computational theories of vision, 2, 2003.
[33] Steven N Evans and Frederick A Matsen. The phylogenetic Miniorovich-Rubinstte metric for environmental sequence samples. Journal of the Royal Statistical Society Series B: Statistical Methodology, 74(3):569-592, 2012.
[34] Makoto Yamada, Yuki Takezawa, Ryoma Sato, Han Bao, Zornitsa Kozareva, and Sujith Ravi. Approximating I-Wasserstein distance with trees. Transactions on Machine Learning Research, 2022.
[35] Tam Le, Makoto Yamada, Kenji Fukumizu, and Marco Cuturi. Tree-sliced variants of Wasser- stein distances. Advances in Neural Information Processing Systems, 32, 2019.
[36] Samantha Chen, Puoya Tabaghi, and Yusu Wang. Learning ultrametric trees for optimal trans- port regression. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pages 20657-20665, 2024.
[37] Ya-Wei Bileen Lin, Ronald R. Coifman, Gal Mishne, and Ronen Talmon. Tree- Wasserstein distance for high dimensional data with latent feature hierarchy. In International Conference on Learning Representations, 2025.
[38] Stephane Mallat. Multiresolution approximations and wavelet orthonormal bases of LR).. Transactions of the American mathematical society, 315 1989.
[39] Stéphane Mallat. A wavelet tour of signal processing. Elsevier, 1999.
12

---

[40] Alfred Haar. Zur theorie der orthogonalen funktionensysteme. Mathematische Anmalen. 71(1):38 53, 1911.
[41] Ingrid Daubechies. The wavelet transform, time-frequency localization and signal analysis. IEEE transactions on information theory, 36(5):911-1005, 1990.
[42] Antonio Ortega, Pascal Frossard, Jelena Kovašević, José MF Moura, and Pierre Vandergheynst. Graph signal processing: Overview, challenges, and applications. Proceedings of the IEEE, (6665):08-828, 2018.
[43] Matan Gavish, Boaz Nadler, and Ronald R Coifmann Multiscale wavelets trees, graphs and high dimensional data: Theory and applications to semi supervised learning. In JCML, 2010.
[44] Minh N Do and Martin Vetterli. The contourlet transform: an efficient directional multires- olution image representation. IEEE Transactions on image processing, 14(12):2091-2106, 2005,
[45] Jindou Dai, Yuwei Wu, Zhi Gao, and Yunde Jia. A peperbooli-to-hyperroe graph convolu- tional network. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 154 163, 2021.
[46] Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009.
[47] Gaspard Monge. Mémoire sur la théorie des déblais des remblais. Mem. Math. Phys. Acad. Royale Sci., pages 666-704, 1781.
[48] Leonid V Kantorovich. On the translocation of masses. In Dokl. Akad. Nauk. USSR (NS), volume 37, pages 199-201, 1942.
[49] Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends in Machine Learning, 11(5-6):355-66 607, 2019.
[50] Ya- fa- Eileen Lin, Ronald R. Coifman, Gal Mishne, and Ronen Talmon. Hyperbolic diffusion embedding and distance for hierarchical representation learning. In International Conference on Machine Learning, pages 21003-21025. PMLR, 2023.
[51] Ronald R Coifman and Stéphane Lafon. Diffusion maps. Applied and computational harmonic analysis, 21(1):5-30, 2006.
[52] Jerrod Isaac Ankenman. Geometry and analysis of dual networks questionnaires. Yale University, 2014.
[53] Matan Gavish and Ronald R Coifman. Sampling, denoising and compression of matrices by coherent matrix organization. Applied and Computational Harmonic Analysis, 33(3):354-369, 2012.
[54] Or Yair, Ronen Talmon, Ronald Coifman, and Ioannis G Kevrekidis. Reconstruction of normal forms by learning informed observation geometries from data. Proceedings of the National Academy of Sciences, 11 E185 2017.
[55] William Leeb and Ronald Coifman. Hölder- norms and their duals on spaces with semigroups, with applications to earth mover's distance. Journal of Fourier Analysis and Applications, 22(4))910-953, 2016.
[56] Nauman Shahid, Nathanael Perraudin, Vassilis Kalofolias, Gilles Puy, and Pierre Van- dergheynst. Fast robust pca on graphs. IEEE Journal of Selected Topics in Signal Processing, 10(4):740-756, 2016.
[57] Gal Mishne, Eric Chi, and Ronald Coifman. Co-manifold learning with missing data. In International Conference on Machine Learming, pages 4605-4614, PMLR, 2019.
[58] Charles Fefferman, Sanjoy Mitter, and Hariharan Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 2994):983-1049, 2016.
13

---

[59] Geert-Jan Huizing, Laura Cantini, and Gabriel Peyré. Unsupervised ground metric learning using Wasserstein singular vectors. International Conference on Machine Learning, pages 9429-9443. PMLR, 2022.
[60] Kira Michaela Düsterwald, Samo Hromadka, and Makoto Yamada. Fast unsupervised ground metric learning with tree- Wasserstein distance. In International Conference on Learning Representations, 2025.
[61] Dingkang Wang and Yusu Wang. An improved cost function for hierarchical cluster trees. arXiv preprint 18Xiiv:8112001 2018.
[62] John Adrian Bondy and Uppaluri Siva Ramachandra Murty. Graph theory. Springer Publishing Company, Incorporated, 2008
[63] Eric C Chi, Genevera I Allen, and Richard G Baraniuk. Convex biclustering, Biometrics, 73(1):10-19 2017.
[64] Eric C Chi and Kenneth Lange. Splitting methods for convex clustering. Journal of Computa- tional and Graphical Statistics, 24(4):994-0013, 2015.
[65] Eric C Chi and Stefan Steinerberger. Recovering trees with convex clustering. SIAM Journal on Mathematics of Data Science, 1(3):383-077, 2019.
[66] Stephen J Wright. Coordinate descent algorithms. Mathematical programming. 151(1)33-3, 2015.
[67] Jerome H Friedman, Trevor Hastie, and Rob Tibshirani. Regularization paths for generalized linear models via coordinate descent. Journal of statistical software, 33:1-22, 2010.
[68] Joe Kileel, Amit Moscovich, Nathan Zelesko, and Amit Singer. Manifold learning with arbitrary norms. Journal of Fourier Analysis and Applications, 27(5):82, 2021.
[69] Ronald R Coifman and M Victor Wickerhauser. Entropy-based algorithms for best basis selection. IEEE Transactions on information theory, 38(2):113-718, 1992.
[70] Albert Cohen, Ingrid Daubechies, and J-C Feauveau. Biorthogonal bases of compactly supported wavelets. Communications on pure and applied mathematics, 45(5):485-560, 1992.
[71] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval. Wavelets on graphs via spectral graph theory. Applied and Computational Harmonic Analysis, 30(2): 2011. 1.
[72] David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst. The emerging field signal processing on graphs: Extending high-dimensional data analysis networks and other irregular domains. IEEE signal processing magazine, 30(3):83-99) 98, 2013.
[73] Nicolas Tremblay and Pierre Borgnat. Graph wavelets for multiscale community mining. IEEE Transactions Signal Processing, 62(20):5227-5239, 2014.
[74] Nora Leonardi and Dimitri Van De Ville. Tight wavelet frames multislice graphs. IEEE Transactions on Signal Processing, 61(13):335733367 2013.
[75] Xiaowen Dong, Dorina Thanou, Pascal Frossard, and Pierre Vandergheynst. Learning laplacian matrix in smooth graph signal representations. IEEE Transactions on Signal Processing. (M(23):6160-617) 2016.
[76] Chao Shen and Hau-Tieng Wu. Scalability and robustness of spectral embedding: Landmark diffusion is all you need. Information and Inference: A Journal of the IMA, 11(4):1527-1555 2022.
[77] Marco Cuturi and David Avis. Ground metric learning. The Journal of Machine Learning Research, 15(1):533-564 2014.
[78] Roy M Anderson and Robert M May. Infectious diseases of humans: dynamics and control. Oxford university press, 1991.
14

---

[79] Octavian Ganea, Gary Bécigneul, and Thomas Hofmann. Hyperbolic neural networks. Ad- vances in Neural Information Processing Systems, 31, 2018.
[80] Ryohei Shimizu, Yusuke Mukuta, and Tatsuya Harada. Hyperbolic neural networks++. arXiv preprini arXiv:2006.088. 0822000000000,,.. 2020.
[81] Jean-Luc Starck, Fionn Murtagh, and Jalal M Fadili. Sparse image and signal processing: wavelets, curvelets, morphological diversity. Cambridge university press, 2010.
[82] Albert Cohen, Wolfgang Dahmen, Ingrid Daubechies, and Ronald Vore. Tree approximation and optimal encoding. Applied and Computational Harmonic Analysis, 11(2):12-2226, 2001.
[83] Emmanuel Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on information theory, 52(())489-51) -09899955990000 2006.
[84] Bianca Dumitrascu, Soledad Villar, Dustin G Mixon, and Barbara E Engelhardt. Optimal marker gene selection for cell type discrimination in single cell analyses. Nature communica- tions, 12(1):1-8, 2021.
[85] David M Blei, Thomas L Griffiths, and Michael Jordan. The nested chinese restaurant process and bayesian nonparametric inference of topic hierarchies. Journal of the ACM (JACM). 57(2):1-30 2010.
[86] Junyue Cao, Malte Spielmann, Xiaojie Qiu, Xingfan Huang, Daniel M Ibrahim, Andrew J Hill, Fan Zhang. Stefan Mundlos, Lena Christiansen, Frank J Steemers, et The single-cell transcriptional landscape of mammalian organogenesis. Nature, 566(7745):49-502, 2019,
[87] Lihong Peng. Xiongfei Tian, Geng Tian, Junlin Xu, Xin Huang, Yanbin Weng, Jialiang Yang, and Liqian Zhou, Single-cell ma-seq clustering: datasets, models, and algorithms. RNA biology, 17(6))765-783, 2020.
[88] Arturs Backurs, Yihe Dong, Piotr Indyk, Ilya Razenshteyn, and Tal Wagner. Scalable nearest neighbor search for optimal transport. In International Conference on Machine Learning, pages 497-506. PMLR, 2020.
[89] Robert Clay Prim. Shortest connection networks and some generalizations. The Bell System Technical Journal, 1666):1389-1401, 1957.
[90] Nicholas Monath, Manzil Zaheer, Daniel Silva, Andrew McCallum, and Amr Ahmed. Gradient- based hierarchical clustering using continuous representations of trees in hyperbolic space. In Proceedings of the 25th ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pages 714-722, 2019.
[91] Giovanni Chierchia and Benjamin Perret, Ultrametric fitting by gradient descent. Advances in neural information processing systems, 32, 2019.
[92] Ines Chami, Albert Gu, Vaggos Chatziafratis, and Christopher Ré. From trees to continuous embeddings and back: Hyperbolic hierarchical clustering. Advances in Neural Information Processing Systems,
2020.
[93] Dmitri Krioukov, Fragkiskos Papadopoulos, Maksim Kitsak, Amin Vahdat, and Marián Boguná. Hyperbolic geometry of complex networks. Physical Review E-Statistical, Nonlinear, and Soft Matter Physics, 82(3):036106, 2010.
[94] Qi Liu, Maximilian Nickel, and Douwe Kiela. Hyperbolic graph neural networks. Advances in Neural Information Processing Systems, 32, 2019.
[95] Zhilin Yang. William Cohen, and Ruslan Salakhudinov. Revisiting semi-supervised learning with graph embeddings. In International conference on machine learning, pages 40-48, PMLR, 2016.
[96] Thomas Kipf and Max Welling, Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations, 2016.
15

---

[97] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. arXiv preprint arXiv:1710.10003, 2017.
[98] Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. Advances in Neural Information Processing Systems, 30, 2017.
[99] Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger. Simplifying graph convolutional networks. In International conference on machine learning. pages 6861-6871. PMLR, 2019.
[100] Clément Bonet, Paul Berg, Nicolas Courty, François Septier, Lucas Drumetz, and Minh-Tan Pham. Spherical sliced-wael wasser-esssrt arXiv preprint arXiv:2206.08788) 08780, 2022.
[101] Albert Gu, Frederic Sala, Beliz Gunel, and Christopher Ré. Learning mixed-curvatie repre- sentations in product spaces. International conference on learning representations, 2018.
[102] Thomas Dagès, Simon Weber, Ya-Wei Eileen Lin, Ronen Talmon, Daniel Cremers, Michael Lindenbaum, Alfred M Bruckstein, and Ron Kimmel. Finsler multi-dimension scaling: Manifold learning for asymmetric dimensionality reduction and embedding, In Proceedings of the Computer Vision and Pattern Recognition Conference, pages 25842-2588 2025.
[103] Thomas Dagès, Michael Lindenbaum, and Alfred M Bruckstein. Metric convolutions: A unifying theory to adaptive convolutions. arXiv preprint arXiv:2406.05400, 2024.
[104] Simon Weber, Thomas Dagès, Maolin Gao, and Daniel Cremers. Finsler-Laplace-Belll Belltrami operators with application shape analysis. In Proceedings of the IEEE/CVF Conference Computer Vision and Patlern Recognition, pages 3131-3140, 2024.
[105] Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine, 34(4):18 42, 2017.
International Conference on Machine Leurning, pages 10086-10095. PMLR, 2021. [106] Yuki Takezawa, Ryoma Sato, and Makoto Yamada. Supervised tree- Wasserstein distance. In
[107] Mikhael Gromov. Hyperbolic groups. In Essays in group theory, pages 75-263. Springer, 1987.
[108] Thomas H Cormen, Charles Leiserson, Ronald L Rivest, and Clifford Stein. Introduction to algorithms. MIT press, 2022.
[109] Martin R Bridson and André Haefliger. Metric spaces of non-positive curvature, volume 319. Springer Science & Business Media, 2013.
[110] John G Ratcliffe, S Axler, and KA Ribet. Foundations of hyperbolic manifolds, volume 149. Springer, 1994.
[111] FeiWei Fileen Lin, Yuval Kluger, and Ronen Talmon. Hyperbolic Procrustes analysis using Riemannian geometry. Advances in Neural Information Processing Systems, $4:5959-59771 2021.
[112] Alan F Beardon. The geometry of discrete groups, volume 91. Springer Science & Business Media, 2012.
[113] Amit Singer. From graph to manifold Laplacian: The convergence rate. Applied and Compula- tional Harmonic Analysis, 21(1): 12-3144,, 2006.
[114] Mikhail Belkin and Partha Niyogi. Towards theoretical foundation for Laplacian- manifold methods. Journal of Computer and System Sciences, 14(8):1289-1308 2008
[115] Elad Lavi, Amir Bourvine, Ya-Wei Eileen Lin, and Ronen Talmon. Hyperbolic distance based on emd and diffusion for hyperspectral imaging. In ICASSP 2025-2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1-5. IEEE, 2025.
16

---

[116] Maurice Fréchet. Les éléments aléatoires de nature quelconque dans un espace distancié. In Annales de l'institut Henri Poincaré, volume 10, pages 215-310, 1948.
[117] Gregory Beylkin, Ronald Coifman, and Vladimir Rokhlin. Fast wavelet transforms and numerical algorithms Communications on pure and applied mathematics, 44(2):41-183,,
1991.
[118] Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich. Displacement interpolation using Lagrangian mass transport. In Proceedings of the 2011 SIGGRAPH Asia conference, pages 1-12, 2011.
[119] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in Neural Information Processing Systems, 26, 2013.
[120] Lenaic Chizat, Pierre Roussillon, Flavien Léger, François- Saavirr Vialard, and Gabriel Peyré. Faster Wasserstein distance estimation with the Sinkhorn divergence. Advances in Neural Information Processing Systems, 33:2257-2269, 2020.
[121] Johannes Gasteiger, Marten Lienen, and Stephan Günnemann. Scalable optimal transport in high dimensions for graph distances, embedding alignment, and more. In International Conference on Machine Learning, pages 5616-5627. PMLR, 2021.
[122] Max Sommerfeld, Jörn Schrieber, Yoav Zemel, and Axel Munk. Optimal transport: Fast probabilistic approximation with exact solvers. Journal of Machine Learning Research, 20(105):1-23, 2019.
[123] Jason Altschuler, Jonathan Viles-Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. Advances in neural information processing systems, 30, 2017.
[124] David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams. Convolutional networks on graphs for learning molecular fingerprints. Advances in neural information processing systems, 28, 2015.
[125] Justin Gilmer. Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George Dahl. Neural message passing for quantum chemistry. In International conference on machine learning, pages 1263-1272. PMLR, 2017.
[126] Xiao Wang, Houye Ji, Chuan Shi, Bai Wang, Yanfang Ye, Peng Cui, and Philip Yu. Het- erogeneous graph attention network. In The world wide web conference, pages 2022-2032, 2019.
[127] Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec. Graph convolutional neural networks for web-scale recommender systems. In Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pages 974 983, 2018.
[128] Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh Pappu, Karl Leswing, and Vijay Pande. Moleculenet: benchmark for molecular machine learning. Chemical science, 9():51135330 2018.
[129] Gal Mishne, Ronen Talmon, Israel Cohen, Ronald R Coifman, and Yuval Kluger. Data-driven tree transforms and metrics. IEEE transactions on signal and information processing over networks, 4(3):451-466) 2017.
[130] Julien Rabin, Gabriel Peyré, Julie Delon, and Mare Bernot. Wasserstein barycenter and its application to texture mixing. In Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, I, Ein-Gedi, Israel, May 29-June 2011, Revised Springer, 2012.
Selected Papers 3, pages
[131] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and radon wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51:22 45, 2015.
17

---

[132] Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde. Gener- alized sliced wasserstein distances. Advances neural information processing systems, 32. 2019.
[133] Clément Bonet, Benoit Malézieux, Alain Rakotomamonjy, Lucas Drumetz, Thomas Moreau, Matthieu Kowalski, and Nicolas Courty. dliced-wasserstel on symmetric positive definite matrices for m/eeg signals. In International Conference on Machine Learning, pages 2777- 2805, PMLR, 2023.
[134] Clément Bonet, Laetitia Chapel, Lucas Drumetz, and Nicolas Courty. Hyperbolic sliced- Wasserstein via geodesic and horospherical projections. In Topological, Algebraic and Geo- metric Learning Workshops 2023, pages 334-370. PMLR, 2023.
[135] Matthieu Heitz, Nicolas Bonneel, David Coeurjolly, Marco Cuturi, and Gabriel Peyré. Ground metrie learning on graphs. Journal of Mathematical Imaging and Vision, 63:89-107, 2021.
[136] Gloria Zen, Elisa Ricci, and Nicu Sebe. Simultancous ground metric learning and matrix factorization with earth mover's distance. In 2014 22nd International Conference on Pattern Recognition, pages 3690-3695, IEEE, 2014.
[137] Fan Wang and Leonidas Guibas. Supervised earth mover's distance learning and its computer vision applications. In European Conference on Computer Vision, pages 442 455. Springer, 2012.
[138] Ruilin Li, Xiaojing Ye, Haomin Zhou, and Hongyuan Zha. Learning to match via inverse optimal transport. Journal of machine learning research, 20(80)1 2019.
[139] Andrew M Stuart and Marie-Therese Wolfram. Inverse optimal transport. SIAM Journal on Applied Mathematics, (01)):59-6199 2020.
[140] Heinz H Bauschke, Patrick L. Combettes, Heinz H Bauschke, and Patrick L Combettes. Correction to: convex analysis and monotone operator theory in Hilbert spaces. Springer, 2017.
[141] Alexander Grigoryan. Heat kernel and analysis on manifolds, volume 47. American Mathe- matical Soc., 2009.
[142] Yiding Zhang, Xiao Wang. Chuan Shi, Xunqiang Jiang, and Yanfang Ye. Hyperbolic graph attention network. IEEE Transactions on Big Data, 1600-171,, 2021.
[143] fa-Wei Eileen Lin, Yuval Kluger, and Ronen Talmon. Hyperbolic diffusion Procrustes analysis for intrinsic representation of hierarchical data sets. ICASSP 2024-2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 6325-6329. IEEE, 2024.
[144] Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeff Dean. Distributed repre- sentations of words and phrases and their compositionality. Advances in Neural Information Processing Systems, 26, 2013.
[145] Amit Zeisel, Ana B Muñoz-Manchado Simone Codeluppi, Peter Lönnerberg, Gioele La Manno, Anna Juréus, Sueli Marques, Hermany Munguba, Liqun He, Christer Betsholtz, et al. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science, 347(6226): 138-11121 2015.
[146] Marlon Stoeckius, Christoph Hafemeister, William Stephenson, Brian Houck-Loomen Pratip Chattopadhyay, Harold Swerdlow, Rahul Satija, and Peter Smibert. Simultancous epitope and transcriptome measurement in single cells. Nature methods, 14(9):855-868, 2017.
[147] Jingcheng Du, Peilin Jia, Yulin Dai, Cui Tao, Zhongming Zhao, and Degui Zhi. Gene2vec: distributed representation of genes based on co-expression. BMC genomics, 20:7-15, 2019.
[148] Ofir Lindenbaum, Moshe Salhov, Arie Yeredor, and Amir Averbuch. Gaussian bandwidth selection for manifold learning and classification. Data mining and knowledge discovery, 34:1676-1712, 2020.
18

---

[149] Pablo A Jaskowiak, Ricardo JGB Campello, and Ivan G Costa. On the selection of appropriate distances for gene expression data clustering, In BMC bioinformatics, volume 15, pages 1-177 Springer, 2014.
[150] Tom Kenter and Maarten De Rijke. Short text similarity with word embeddings. In Proceedings of the 24th ACM international on conference on information and knowledge management, pages 1411 1411-1420, 2015.
[151] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference knowledge discovery & data mining, pages 2623-2631, 2019.
[152] Lihi Zelnik-Manor and Pietro Perona. Self-tuning spectral clustering. Advances in Neural Information Processing Systems, 17, 2004.
[153] Xiucai Ding and Hau-Tieng Wu. Impact of signal-t-noosee ratio and bandwidth graph Lapla- cian spectrum from high-dimensional noisy point cloud. arXiv preprint arXiv:2011.10725, 2020.
[154] Ronald R Coifmaa and Matan Gavish, Harmonic analysis of digital data bases. Wavelets and Multiscale Analysis: Theory and Applications, pages 161-197, 2011.
[155] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. Pot: Python optimal transport. The Journal of Machine Learning Research, 22(1):3571-3577 2021.
[156] Chakib Fettal, Lazhar Labiod, and Mohamed Nadif. Efficient and effective optimal transport- based biclustering. In Advances in Neural Information Processing Systems, 2022.
[157] Bruno Cesar Feltes, Eduardo Bassani Chandelier, Bruno lochins Grisci, and Marcio Domn Cumida: an extensively curated microarray database for benchmarking and testing of machine learning approaches in cancer research. Journal of Computational Biology, 2664):376-386, 2019.
[158] John A Hartigan and Manchek A Wong, Algorithm AS 136: A k-means clustering algorithm. Journal of the royal statistical society. series e (applied statistics), 28(1):100-108 1979.
[159] Charlotte Laclau, levgen Redko, Basarab Matei, Younes Bennani, and Vincent Brault. Co- clustering through optimal transport. International Conference on Machine Learning, pages 1955-1964. PMLR, 2017.
[160] Vayer Titouan, Ievgen Redko, Rémi Flamary, and Nicolas Courty. Co-optimal transport. Advances in Neural Information Processing Systems, 33:17559-
2020.
[161] Deng Cai, Xiaofei He, Xiaoyun Wu, and Jiawei Han. Non-negative matrix factorization on manifold. In 2008 eighth IEEE international conference on data mining, pages 63-72. IEEE, 2008.
[162] Lawrence Hubert and Phipps Arabie. Comparing partitions. Journal of classification, 2(1):193- 218, 1985.
[163] Noga Alon, Richard M Karp, David Peleg, and Douglas West. A graph-theoretic game and its application to the k-server problem. SIAM Journal on Computing, 24(1):78-100, 1995.
[164] Facundo Mémoli. Gromov-wasserstein distances and the metric approach to object matching. Foundations of computational mathematics, 11(4):417-477, 2011.
[165] Ching- Yao Chuang and Stefanie Jegelka. Tree mover's distance: Bridging graph metrics and stability of graph neural networks. Advances Neural Information Processing Systems, 35:2944-20 2957, 2022.
19

---

[166] Alexander Grossmann and Jean Morlet. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM journal on mathematical analysis, 15(4):723-736
1984.
[167] Nhat Khang Ngo, Truong Son Hy, and Risi Kondor. Multiresolution graph transformers and wavelet positional encoding for learning long- rrange and hierarchical structures. The Journal of Chemical Physics, 159(3), 2023.
[168] Xiaowen Dong, Dorina Thanou, Laura Toni, Michael Bronstein, and Pascal Frossard. Graph signal processing for machine learning: A review and new perspectives. IEEE Signal processing magazine, (7(6):117-12) 2020.
[169] Robert Grone, Russell Merris, and Viakalathur Shankar Sunder. The laplacian spectrum of graph. SIAM Journal on matrix analysis and applications, 11(2):218-22 238, 1990.
[170] Fan RK Chung. Spectral graph theory, volume 92. American Mathematical Soc., 1997.
[171] Petre Stoica, Randolph L Moses, et al. Spectral analysis of signals, volume 452. Citeseer, 2005.
[172] Aliaksei Sandryhaila and José MF Moura. Discrete signal processing on graphs: Graph fourier transform. In 2013 IEEE International Conference Acoustics, Speech and Signal Processing, pages 6167-6170. IEEE, 2013.
[173] Wei Hu, Jiahao Pang, Xianming Liu, Dong Tian, Chia-Wen Lin, and Anthony Vetro. Graph signal processing for geometric data and beyond: Theory and applications. IEEE Transactions on Multimedia, 14:3961-3977 2021.
20

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We introduce iterative framework for jointly learning hierarchical repre- sentations samples and features using the tree- ree-Wasserstein distance. Our method alter- nates between tree construction and distance computation across modes and incorporates data-driven Haar wavelet filters to represent meaningful hierarchical structure. The main contributions summarized in Sec. 1. Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The paper discusses the limitations of the work in Sec. 7, where we outline the constraints of the proposed method.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model vell-specification asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken low lighting. Or speech-to-texx system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledger in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
21

---

3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Right after each stated theoretical result, we indicate that the proof and the theoretical analysis can be found in App. C.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: App. E provides detailed descriptions of the proposed method. The method can be re-implemented, and the results are reproducible.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
also be provided via detailed
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility.
22

---

In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: We provide detailed descriptions of method's implementation in App. E. Our code is available at https: / /gi thub. wwwccomoooa///c http//aavveeeeeeeeee lin weaee-iii/quusiioooo htt///uuustooniirree
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips cc/
Souiiee/CooeSSummtt awwsonoooiccyyy olicy) for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
/ //nips. wcwppblpiiicccccccccccccccccccba htt//gudde/////e Pilees/oodSuubbsse ssionsolicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: We provide detailed information on data splits, hyperparameter settings, initial used distance, and the criteria used for their selection. The type of classifier and all the training details are also described. This information is App. E.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Sec. 6, we assess learned tree quality by measuring the L1 norm of the Haar coefficients. For document and single-cell classification, we report both the mean and standard deviation of the classification accuracy. For hierarchical graph data, we evaluate node classification and link prediction tasks: the link prediction performance is measured
23

---

using ROC AUC, while the node classification performance is reported using the FI score on the DISEASE dataset and accuracy on the AIRPORT, PUBMED, and CORA datasets.
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: Detailed descriptions of the type of computer resources and time of execution are provided in App. E and App. F.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: :tttpp:://euurre /ttp/eeurpp..// cc
Answer: [Yes]
Justification:
We respect the NeurIPS code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts

---

Question:
Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper presents work whose goal is to advance the field of unsupervised representation learning. The paper has no foreseeable societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but encourage authors to take this into account and make best. faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We provide credit to the original papers of all external data, code, and models used in this work, ensuring that their contributions are acknowledged and their licensing
terms are followed.
25

---

Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Ye]]
Justification:
The new assets introduced in the paper are well documented. Detailed descrip- tions our method's implementation are provided App. E.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: No crowdsourcing nor research with human subjects was involved in this work.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
26

---

Answer:
[NA]
Justification:
No crowdsourcing nor research with human subjects was involved in this work.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs were used for editing purposes (e.g., grammar, spelling, word choice), which does impact the scientific rigor or originality of the work.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non- con-standard components.
Please refer to our LLM policy (https for what should or should not be described.
27

---

Supplementary Material
This supplemental matcrial is organized as follows:
Appendix A contains the additional background to the proposed method. We begin by reviewing the fundamentals of tree structures. Next, we present key concepts from hyperbolic and diffusion geometry, which are used to construct trees from high-dimensional data. We then describe the construction of a Haar wavelet basis on given tree and introduce the Wasserstein distance and tree- Wasserstein distance. Lastly, we briefly discuss graph convolutional networks and their hyperbolic counterpart.
Appendix B provides additional related work relevant to our method. It includes manifold learning based on the Wasserstein distance metric, further discussion of existing TWD methods, the sliced- Wasserstein distance, and additional background on co-manifold learning and unsupervised ground metric learning.
Appendix C contains the proofs of Thm. and Thm. 2, along with the supporting lemmas.
Appendix D describes how the learned hierarchical representations obtained by our method can be incorporated into hyperbolic graph convolutional networks as initialization step.
Appendix E provides details about the experiments presented in Sec. 6. We first describe the datasets used and their statistics. We then report the initial pairwise distances, scaling factors, Haar wavelet thresholds, and other hyperparameters used in our experiments. The norm regularization terms applied in the iterative learning scheme are also specified. We explain how our method is scaled to handle large datasets. We provide details the experimental setups for document classification and cell-type classification, well as for link prediction and node classification on hierarchical graph data.
Appendix F presents additional experimental results supporting our method. We first show a synthetic toy example simulating a video recommendation system for visualization. We then demonstrate the empirical uniqueness of the learned hierarchical representations under strong regularization. An ablation study evaluates the role of the iterative joint learning scheme in classification performance and contrasts it with applying wavelet filtering only after convergence. We further assess the effect of integrating wavelet filtering with alternative TWD variants, well as with HGCNs without the iterative updates. We provide details on the sparse approximation analysis using Haar coefficients across iterations. as well report the classification performance over the iterative learning scheme. Additionally, we include runtime analysis and explore the use of alternative regularizers within the iterative learning scheme. We include additional experiments co-clustering tasks and incorporating fixed hierarchical distances in the iterative learning scheme for HGCNs. We evaluate the effectiiveness of our method as preprocessing step across different neural network backbones. We further include a discussion and comparison with other tree distance metrics. Finally, we present the visualization of different iterations of the learned trees and how they evolve in our method on the toy example.
Appendix G provides additional notes on the motivation for incorporating wavelet filters into our iterative learning scheme. We explain why wavelet filters are considered over Laplacian- based filtering in our context. We clarify the objective of our method and highlight the key differences between our approach and unsupervised ground metric learning with trees (Tree-WSV). We discuss the interpretability of TWDs.
A Additional background
We include supplementary background in the appendix.
A.1
Graph and tree
Trees are fundamental structures in graph theory. A tree a type of graph characterized by hierarchical organization. Below, outline key concepts and definitions of trees and their associated metrics.
28

---

Graph and shortest path metric. A graph G = (V, A) consists of a set of nodes V, edges E, and a edge weight matrix A IR Rmmx For any two nodes j. j' the shortest path metric da.j"")) represents the minimum sum of edge weights along the path connecting node j and node j'.
Tree metric. A metric d: V
R is a tree metric there exists a tree T E, such that d(j.j') corresponds to the shortest path metric du (i.j") on T. Notably, tree metrics exhibit O-hyperbolicity [107], making them well-suited for representing hierarchical data.
Binary trees. Binary trees are a specific class of trees where each internal node has exactly two child nodes. In a balanced binary tree, internal nodes (except the root) have a degree of 3, while leaf nodes have a degree of 1. A rooted and balanced binary tree is a further refinement where one internal node, designated as the root, has degree of 2, and all other internal nodes maintain a degree of 3. Note that flexible tree can be represented as binary tree through transformation, e.g., left-child right-sibling [108].
Lemma A.1. Any tree metric metric derived from a tree graph where the distance between two nodes is the length of the unique path connecting them) can be isometrically embedded into an & space [109].
A.2 Hyperbolic geometry
hyperbolic geometry is a Non-Euclidean geometry with a constant negative curvature. It is widely used in modeling structures with hierarchical tree-like relationships [110]. There are four models commonly used to describe hyperbolic spaces. The Poincaré disk model maps the entire hyperbolic plane inside unit disk, with geodesics represented as ares orthogonal to the boundary straight lines through the center. The Poincaré half-plane model uses the upper half of the Euclidean plane, where geodesics semicireles orthogonal the horizontal axis vertical lines. The Klein model represents hyperbolic space a unit disk but sacrifices angular accuracy, making it more suitable for some specific calculations. The hyperboloid model relies on higher-dimensional Lorentzian inner product, offering computational efficiency.
In this work, we consider two equivalent models of hyperbolic space [110]: the Poincaré half- space and the hyperboloid model. The Poincaré half-space was used for embedding that reveals the hierarchical structure underlying nigh-dimensional data [50], while the hyperboloid model is advantageous for its simple foosed-form Riemannian operations [14, 111] used in hyperbolic graph convolutional networks [21, 45].
The d-dimensional Poincaré half-space defined as
Hd (a E R
| [a)) >}
(11)
with the Riemannian metric tensor ds² = (da² 2 (1) + +... da' [112]. Let dL denote the hyperboloid manifold in dimensions, defined by
Ld = ( # (b, b) 1.,,,,,1))))) 0}.
(12)
where 2 is the Minkowski inner product (b, b) =b= L. Let ToLd be the tangent
space at point E I.d, given by I = (v E R l v)c = 0). We denote ||v = v(v, v) as the norm of v E T.Ld. For two different points b,, b2 E Ld and 0 + v Tb, 14, the exponential and logarithmic maps of L are respectively given by
Expo, (v) cosh([|v]lc))) + sinb(//llle////llllll
(13)
Logb, (bz) cosh 1 (n) - r/bi)/ 2/2
(14)
where 1) - (bb,, ba) c- The parallel transport (PT) of vector v E Ld along the geodesic path from by E Ld to by E Ld is given by
bbbb,-b. = (b2 Ab1, v)2 (b, bz),
(15)
X+1
where À == - (b,, bz) while keeping the metric tensor unchanged. Due to the equivalence between the Poincaré half-space and the Lorentz model, there exists diffcomorphism P. Hdd + Ld that
29

---

maps points from the Poincaré half-space E Hd to the Lorentz model b E L by
b P(a)= (1+ ||c|| 2c(1).... 2c(d + 1))
(16)
|||
where
(2a(1), ,.... 2a(d), ||a| 1)
|a| 1 + 22(d (d + 1) +1
A.3 Diffusion geometry
Difftusion geometry [51] is a mathematical framework that analyzes high-dimensional data by capturing intrinsic geometric structures (i.c., manifolds with non-negative curvature) through diffusion processes. It is rooted in the study of diffusion propagation on graphs, manifolds, or general data spaces, where the spread of information or heat over time reflects the underlying connectivity and geometry. Diffusion geometry represents data as nodes in graph and models their relationships using diffusion operator. By simulating diffusion processes over this structure, the method identifies meaning ful relationships based on proximity and connectivity in the diffusion space. We introduce the construction of the diffusion operator and its desired property below.
in Consider a set of sigh-dimensionall points Z = {z, E R" Jj=1 lying on a low-dimensional manifold. Let K = (apMMMM2// E Rm be an affinity matrix, where M E mm is suitable pairwise is the Hadamard power, and
distance matrix between the points { parameter. Note that the matrix K can be interpreted an undirected weighted graph G (Z, where Z is the node set and K represents the edge weights.
the scale
The diffusion operator [51] is then constructed by P = KD-1 where D is the diagonal degree matrix with entries (j ) = E, K(j,1). We remark that a density normalization affinity matrix be considered to mitigate the effects of non-uniform data sampling [51]. Note that the diffusion operator P is aolumn-sttcchsteie
allowing it to be used as a transition probability matrix of a Markov chain on the graph G. Specifically, the vector D P' is the propagated density after diffusion time E I of a density 6, initially concentrated at point j.
The diffusion operator has been demonstrated to have favorable convergence [51]. As - and c c- 0, the operator pi/e converges pointwise to the Neumann heat kernel associated with the underlying manifold at the limits. This convergence indicates that the diffusion operator can be viewed as discrete approximation of the continuous heat kernel, thereby effectively capturing the geometric structure of the underlying manifold in finite-dimensional setting [51, 113, 114).
A.4 Tree construction based hyperholic and diffusion geometries
Recently, the work [50] introduced a tree construction method that recovers the latent hierarchical structure underlying high-dimensional data 2 based on hyperbolic embeddings and diffusion geome- try. Specifically, this latent hierarchical structure can be modeled by a weighted tree T, which can be viewed a discretization of an underlying manifold This manifold J is assumed to be a complete, simply connected Riemannian manifold with negative curvature, embedded in a high-dimensional ambient space IR", and equipped with geodesic distance dg.
Consider the diffusion operator constructed as in App. A.3, the multi-scale diffusion densities 1 p2_* 8, are considered with dyadic diffusion time steps = 2 for k E z [50, 115], and are embedded into a Poincaré half-space by
(j,k) +y yブ [11 22/2-22
11
(17)
where 1 2012 is the Hadamard root of w.. Note that diffusion geometry effective in recovering the underlying manifold [51]; however, the study [50] demonstrated that considering propagated densities with diffusion times on dyadic grids can capture the hidden hierarchical structures by incorporating local information from exponentially expanding neighborhoods around each point. The multi-scale hyperbolic embedding is a function Embedding Z + M defined by
Embedding(z.) = [ ,y²))
(18)
30

---

Algorithm 3 Tree Construction Using Hyperbolic and Diffusion Geometry [37]
Input: Pairwise distance matrix M e R=x" the maximal scale K, and the scale parameter
Output: A binary tree B with m leaf nodes
function T(M) K < exp( -222/)) D t diag K), P KD-1
for k 10.1...., K] and j E [m] do
2 -P22 8j
(1) /
y)
B < leaves (j) j [m])
for j / E [m] do
for k (0. 1..1.... K} do
1//1//2 1g. 0//2) T ,2//2-2
proj(y) V VyJ)
Onj'" (K+1) proj (v} V vygg) proj(y) v vyf))
S {(j sorted by az,r)
for (j.j') j') S do
if and j' are not in the same subtree Zj < internal node for node j, Ip + internal node for node j'
add internal node for I, and Ly, and assign the geodesic edge weight
return B
where M m++ H x mH of (K 1) elements. The geodesic distance induced in M is the & distance on the product manifold M, given by
dmll.g') = Σ 2sinh (aanh-((--t (2-h/2211 || - v1111)).
(19)
Theorem A.1 (Theorem [50]). For sufficiently large K, the embedding distance is bilipschitz equivalent to the underlying tree distance, i.e., d.M = dr.
The multi-scale hyperbolic embedding is used to construct binary tree ((M) with leaves [37], where each leaf corresponds to a data point in Z. Pairs of points are merged based on the Riemannian mean [116] of the radius of the geodesic connecting the hyperbolic embeddings y and y for
k {O..... K}. Specifically, the k-th hyperbolic embedding, the orthogonal projection the m + 1-axis in Rm++1 which is the radius of the geodesic of y and y is given by
proj(y} v vyt)
101/2
pgg 02/2) T 24/2-22
(20)
2
The Riemannian mean of the orthogonal projection on the + 1-axis in 1m++1 across all K + 1 hyperbolic embeddings has a closed form
h IO...., 0, a, where ajs-
proj(y)) v vy.) projoy y
(21)
The values aj.j' serve as hierarchical linkage scores and guide the merging process in constructing the binary tree B. Tree edge weights are then assigned using the & distance on the manifold M, which corresponds to the geodesic distance. The tree construction based on hyperbolic embedding and diffusion densities summarized in Alg.
Theorem A.2 (Theorem [37].)- For sufficiently large K and m, the tree metric d in Alg. 3 is bilipschitz equivalent to the underlying tree metric dr, where the ground truth latent tree and dr is the hidden tree metric between features defined as the length of the shortest path on T.
A.5 Tree-Based wavelet
31

---

Haar wavelets can be constructed directly and adaptively from trees straightforward manner [43]. Without loss of generality, we focus here the Haar basis induced by the feature tree, where the resulting wavelet filter is subsequently applied to the samples. Symmetrically, one can construct a Haar basis from the sample tree and apply the corresponding wavelet filter to the features.
<=2
=
911
Consider a scalar function f X + data matrix. Let S = {flif X - }R be the space of all functions on the dataset. For a complete binary tree T(M)) with m leaves, let & = I,..., L. denote the levels the tree, where e leaf level with leaves. Define ((f,e) E) be the of all leaves in the E-th subtree of ((M)) at level é. Define St as the space of features that are constant across all subtrees at level f, and let 1x be a constant function on X with the value We have the following hierarchy
defined on the
the root level and L T.e is the
Figure 3: An illustration Haar basis induced by a binary tree.
S Span(1x); She = S, and S C...CS.. (22)
Therefore, the space S can be decomposed
A L.-1
S
Q¹)
D St
(23)
where Q" is the orthogonal complement of St. As the tree T(M)) is full binary tree, the Haar-like basis constructed from the is essentially the standard Haar basis [40]. denoted by {,er, R")
At level f, a subtree (e, 8) splits into two sub-subtrees (e + 1, 81) and Pe + 1, 82). A zero- mean Haar wavelet Br,. E R has non-zero values only at the indices corresponding leaves in the sub-subtrees and is piecewise constant on each of them. The set of these Haar wavelets, along with constant vector Bo, complete and forms an orthonormal Haar basis. We collect all of these basis vectors as columns of a matrix, denoted by B E m Fig. 3 illustrates a Haar basis induced by binary tree.
Proposition A.1 (Function Smoothness and Coefficient Decay [43]). Let Be,c be the Haar basis function supported on T(E,s), If the fanction f is Lipschitz continuous, then for some C >> O,
KF, (r,Br,))sl S (CdTr,ss(T(eee))
(24)
where (T(M))((( is the tree distance between the internal node rooted at e-subiree to the leaves.
Prop. A. indicates that the smoothness of the samples (X,.) leads to exponential decay rate of their wavelet coefficients as a function of the tree level. Consequently, the wavelet coefficients can serve as measure to assess the quality of the tree representation T(M.).
Proposition A.2 (L, Sparsity [43]). Consider a Haar basis Be, where O C S and such that [Be (i)l < 1/ |0 1/2 Let f= Eo aeße and assume Ee |Bal I C. For any O, the approxi- mation f Σ Ejo|en aghe, then
||-
Σ ( S SCVR. JE[m)
(25)
Prop. A.2 implies that with L, approximation, is sufficient to estimate the coefficients for function approximation. In the Haar domain, estimating these coefficients be achieved using the fast wavelet transform [117].
Note that when two trees share the same branching, assigning different weights to their leaves changes the inner product that defines orthogonality. The inner product changes, each wavelet rebalances the weighted masses of its children, and coefficient magnitudes scale with the square roots of those
32

---

masses, Because the basis is built by successive weighted contrasts, altering the weights cascades through all levels: coarse vectors that once averaged equal-sized parts may now emphasize one side, while fine-scale vectors adjust in the opposite direction to keep the basis orthonormal. Consequently, although the supports of the basis vectors coincide, the weighted tree produces a non-equivalent orthonormal basis and, in turn, different expansions for the data.
A.6 Wasserstein distance
The Wasserstein distance [46] measures the discrepancy between probability distributions over given metric space. It is related to optimal transport (OT) theory [47, 48], where the goal is to quantify the cost of transforming one distribution into another.
Consider two probability distributions μ and y defined on a metric space (X, d), where is the ground metric. The Wasserstein distance is defined as
f Jxxx d(z, d)) y).
OT d,
inf TEI((a.)
(26)
where [(,,,) v) denotes the set of all joint probability measures on X' x & with marginals μı and v, i.e.,
(A x) = ((A), B r(B),
(27)
for all measurable sets A, B C X. It represents the minimal "cost" to transport the mass of μ to which is also known as "earth mover's distance" [49, 43].
In the discrete setting, the Wasserstein distance is commonly applied to two discrete probability distributions μ Σ 1402, and v/ Σラ= Py where 82 denotes the Dirac delta function at I, and ( { are the support points of μl and v, respectively. The distributions satisfy Σ M4 = プド= Vj = In this casc, the Wasserstein distance can be formulated as the solution to the following linear programming problem:
OT(H,, v, min ΣΣ dr
(28)
DET(UM) 1=1
where r Ruxm is the transport matrix satisfying the marginal constraints:
Σ Ta Mi: I Vi.j.
(29)
When computing the Wasserstein distance between discrete probability distributions, the compu- tational bottleneck often lies in constructing and processing the ground pairwise distance matrix. This matrix encodes the distances between every pair of discrete support points in the distributions, requiring ((m²) storage and O(m3 log m) computation [49]. This complexity arises due to solving linear program for the optimal transport problem, which scales poorly with the number of points m. Consequently, applying optimal transport directly becomes infeasible for large-scale datasets [1 18).
Several approximations have been proposed to reduce this computational complexity. The Sinkhorn distance [1 19, 120] introduces an entropy regularization term to the objective function, enabling the use of iterative matrix scaling algorithms in quadratic. Graph-based methods exploit sparsity in the transport graph simplify the problem structure, while sampling-based approaches approximate distributions using a subset of support points, reducing computational demand [121-123].
A.7 Tree- Wasserstein distance
The Wasserstein distance requires solving a linear programme whose computational cost is quadratic in the number of support points. A widely used strategy to mitigate this computation cost is to approximate the ground metric with tree metric, thereby defining the tree- Wasserstein distance (TWD) [32]. Because transport on a tree admits a closed-form solution, TWD can be computed in (m)) time for tree with m leaves, which makes attractive for large-scale applications.
Let T (V, F, A) be tree with Nica leaves. The tree distance dr V x V + R is the sum of weights of the edges the shortest path between any two nodes on T. Let Tr(w) ) be the set of nodes
33

---

in the subtree of the tree T rooted at the node E V. For any node E V. there exists a unique parent s.t. the weight is defined by the distance, i.e., Wa dr(u, v) [62]. Given two discrete distributions Pl; P2 E NR supported the tree T, the TWD is formally defined by
TW P2, = Evev We EueTr(e) (pi(u) (anu))p(uu))
(30)
where inner sum represents the net mass that must cross the edge (v, parent weighting by Wer which accumulates the transport cost across the tree
Given tree 7 and probability measures supported on the tree T, the TWD computing using T' is the Wasserstein distance the tree metric dr [35, 34), i.e.,
O'((2 P2, dT) TW P2, T).
(31)
However, is important to note that most of the TWD methods [32-36] are designed approximate Wasserstein distance with an Euclidean ground metric. Specifically, their goal is
min (0T(PI,P2,de P2. (P1 P2. T2)2? T
(32)
where dE denotes the Euclidean metric. Therefore, the tree construction in these methods is used to approximate the Euclidean ground metric. This design choice, while useful for accelerating computation, inherently biases the tree structure toward approximating Euclidean distances rather than representing the hierarchical structure of the data, which is the primary focus of our work.
A.8 Graph convolutional networks
Graph Convolutional Networks (GCNs) [96] have gained significant attention in graph machine learning, where nodes graph are typically assumed in Euclidean spaces. By generalizing con- volutional operations to graphs, GCNs effectively capture the dependencies between nodes. This capability enables high accuracy in tasks such as node classification, link prediction, and graph classification, making GCNs useful in applications like social networks [98], molecular structures [124, 125]. knowledge graphs [126], recommendation systems [127], and drug discovery [128]. A brief overview of the GCN framework is presented below.
Consider a graph G (V, A) with node features (x; E J where V is the vertex set containing nodes, E C V x V is the edge set, and A E IR=x" is the adjacency matrix. Each node i V is associated with feature vector x; E T representing m-dimensional node attributes. In each layer of GCN message passing, the graph convolution can be divided into two steps: feature transformation and neighborhood aggregation. Specifically, the feature transform defined by
h J h (8-1)
(8)
(33)
where h, (0)
is the initial feature, and ( is learnable weight matrix. The neighborhood aggregation then updates the representation by
(2) h
(6) h
Σ +E[[S]]
wis
wish! h(9
(34)
where [[i]] denotes the neighbors of node i, wis a weight associated with the edge between nodes and 8, and (-)) is non-linear activation function.
By stacking multiple such layers, GCNs propagate information through the graph, enabling each node representation to integrate signals from multi-hop neighborhoods. The feature transformation step learns task-specific embeddings, while the aggregation step incorporates structural information from the graph topology.
A.9 Hyperbolic graph convolutional networks
Hyperbolic Graph Convolutional Networks (HGCNs) [21, 45] generalize Graph Convolutional Networks (GCNs) [96] to hyperbolic spaces, using the inductive bias of hyperbolic geometry to better capture hierarchical structures in graph data. By embedding nodes in negatively curved space,
34

---

HGCNs provide improved representational capacity for graphs with hierarchies, while preserving the scalability of standard GCNs. This extension requires adapting the steps in GCNs, including feature transformation and neighborhood aggregation, to conform with the geometry of hyperbolic space. HGCNs have shown empirical advantages across a range of tasks, including social network analysis, recommendation systems, and biological network modeling, demonstrating the effectiveness of geometric inductive biases deep learning on non-Euclidean domains. We briefly review the HGCNs framework below.
HGCNs extend GCNs to hyperbolic spaces (e.g.. Lorentz model) by embedding node features in hyperbolic geometry and redefining core operations to respect its geometric structure. Below, we outline the key components of the HGCN architecture.
Mapping Euclidean node features to hyperbolic representations. HGCNs [21] begin by mapping input Euclidean node features Lorentz model using the exponential map at the origin. Specifically, given a Euclidean feature vector X; E I,m,, its hyperbolic embedding is initialized as
H!(() ExPo (10. x 1') E L.
(35)
where the Exp is the exponential map.
Hyperbolic feature transformation. To perform feature transformations in hyperbolic space, HGCNs [21] use the logarithmic and exponential maps to move between the manifold and the tangent space at the origin. Let W E Rm x be learnable weight matrix. The hyperbolic equivalent of matrix multiplication defined as
W x # Expo W (Logo ")).
(36)
To incorporate a bias term b E Rmn HGCNs use parallel transport and define hyperbolic bias addition
as
b Expx, (PT (b)),
(37)
where PT is the parallel transport operator. Then, the full hyperbolic feature transformation is defined as
11 (w()))
x (-1).") + (0)
(38)
Hyperbolic neighborhood aggregation. Neighborhood aggregation HGCNs performed in the langent space and leverages hyperbolic attention to account for the hicrarchical relationships within the graph. Given hyperbolic embeddings x and x both are first mapped to the tangent space the origin. The attention weights are computed via an Euclidean Multi-Layer Perceptron (MLP) applied the concatenation of the tangent vectors
Wiy softmax (MTP(Togoox")) [[ooo (xf ))),
(39)
where [[:]] denotes the neighbors of node i, and || indicates concatenation. Aggregation is then performed via a hyperbolic weighted average in the tangent space of x
AGG(x = Expxx"
Σ jElN
WiyLoga (xf*)
(40)
Then, a non-linear activation applied in hyperbolic space to complete the layer's update.
Remark A.1. In the hyperbolic feature transformation step of HGCNs [21], the standard approach involves mapping points between the hyperbolic manifold and the Euclidean tangent space via loga- rithmic and exponential maps. This allows feature transformations to be performed using Euclidean operations, such as linear transformations followed by bias addition. However, this design only partially uses the advantage of the geometry of hyperbolic space, the actual transformation takes place in the flat tangent space rather than directly on the curved manifold. To exploit the structure of hyperbolic space, [45] proposed fully hyperbolic feature transformation. Specifically, they introduce a Lorentz linear transformation that operates directly in the hyperbolic space
1, 0
x (). [[0 (()x((1))) (11 s.t. (
and (J =I
(41)
0,
35

---

Here, J(t is an orthogonal matrix, while the first coordinate (the time-like component in Lorentz space) remains unchanged. This formulation preserves the hyperbolic structure and allows feature transformations to be conducted entirely within the manifold. When integrating our method into the HGCN framework, we adopt this syperbolietoo-hyperbooll feature transformation to maintain geometric consistency.
B Extended related works
Manifold learning with Wasserstein distance. Manifold learning is nonlinear dimensionality reduction approach designed for nigh-dimensional data that intrinsically lie on a lwwer-dimensional manifold [58, 51, 1 114]. Most classical methods derive manifold representations from discrete data samples by constructing graph, where nodes correspond to data points and edges reflect pairwise similarities. These similarities are often defined using Gaussian affinity kernel based on the Euclidean norm between the data samples. Recent studies have demonstrated the benefits of using the Wasserstein distance more geometry-aware alternative to Fuclidean distances in manifold learning frameworks [43, 52], where the Wasserstein distance incorporates the feature relationship. This direction has shown promising results in applications such as matrix organization [129], analysis of neuronal activity patterns [7]. and molecular shape space modeling [68]. Building on these ideas, our work explores the integration of the diffusion operator [51] with the TWD to jointly learn hierarchical representations for both samples and features.
rree-asserrtten distance. The Tree- Wasserstein Distance (TWD) offers a computationally ef- ficient alternative to the Wasserstein distance by approximating the ground cost metric with a tree structure. It is designed to compare probability distributions by quantifying the amount of mass that must be transported between them, where the transport is governed by distances along a tree. By approximating the original Euclidean ground metric with tree metric, TWD enables faster computation while maintaining a close approximation to the true Wasserstein distance. This trade-off between efficiency and fidelity has been validated in range of studies 323. where TWD has been shown to provide a reliable and scalable surrogate for optimal transport in Euclidean settings. For instance, the Quadtree [32] recursively partitions the ambient space into hypercubes build random trees, which are then used to calculate the TWD. Flowtree [88] refines this approach by focusing on optimal flow and its cost within the ground metric. Similarly, the Tree-Sliced Wasserstein Distance (TSWD) [35] further improves robustness by averaging TWD values computed over multiple ran- domly sampled trees, with variants such as TSWD-1, TSWD-5, and TSWD-10 corresponding the number of sampled trees used. Recent advancements include WQTWD and WCTWD [34]. which employ Quadtree or clustering- bbased tree structures and optimize tree weights to approximate the Wasserstein distance. UltraTree [36] introduces an ultrametric tree by minimizing OT regression cost, aiming to approximate the Wasserstein distance while respecting the original metric space.
Most recently, a tree construction method [37] in TWD literature that differs fundamentally from existing TWD approaches [32-35]. Whereas conventional TWD methods focus how TWD can be close to the true Wasserstein distance with the Euclidean ground metric, this method [37] aims to efficiently compute Wasserstein distance with a ground metric that recovers latent hicrarchical structure underlying hierarchical high-dimensional features. Unlike prior TWD methods that treat tree construction as a heuristic or approximation step to Euclidean distances, this approach [37] treats the tree as a geometric object that reflects intrinsic hierarchical relationships among features. In particular, the constructed tree is not merely a proxy for an ambient Euclidean metric, but rather data-driven structure whose tree distances approximate geodesic paths on latent, unobserved hierarchical metric space. This allows the resulting TWD to encode meaningful relationships in settings where the data are governed by hidden hierarchies. In our work, we adopt this tree construction method [37], as it provides a principled and geometry-aware foundation for learning hierarchical representations of both samples and features.
Sliced- Wasserstein distance. Another widely used approach for efficiently approximating the Wasserstein distance is the -aseedd-Wssrsttein Distance (SWD) [130, 131]. Instead of solving the high-dimensional optimal transport problem directly, SWD projects the input distributions onto one- dimensional subspaces along random directions, where the Wasserstein distance admits a closed-form solution. Averaging the results over many such slices provides meaningful measure of similarity while retaining much of the geometric structure inherent in the distributions. This approach is
36

---

advantageous in high- tigh-dimensional settings, where computational efficiency is critical, making it popular applications such as generative modeling [132], domain adaptation [133], and statistical inference [134].
Co-Manifold learning. Co-manifold learning aims to simultaneously uncover the manifolds of both samples and features in a data matrix by treating their relationships mutually informative. That is, the geometry of the features informed by the samples, and vice versa. This joint modeling approach has been widely explored in tasks such as joint embedding [52-55] and dimensionality reduction [7, 56, 57], under the assumption that both samples and features lie on smooth, low-dimensional manifolds. By iteratively refining representations across both modes, co-manifold methods aim to capture the manifolds of both samples and features. These approaches typically assume that the manifolds involved are Riemannian and characterized by non-negative curvature [55, 58]. However, such assumptions may he limiting in scenarios where the data exhibit hierarchical organization, which is more naturally modeled in negatively curved spaces such as hyperbolic geometry.
In our work, we adopt the co-manifold learning perspective, using the geometry of the samples to inform the structure of the features and vice versa. However, unlike existing co-manifold approaches that primarily focus on smooth manifolds with non-negative curvature, our focus lies in capturing the hierarchical nature of both samples and features. By modeling these hierarchies explicitly through trees and employing the corresponding TWDs for cross-mode inference, we address a distinct setting where the data is nigh-dimensional with hicrarchical structures.
Unsupervised ground metric learning. Ground metric learning [77] focuses on learning a distance function that defines the cost matrix in the optimal transport problem, thereby capturing meaningful relationships between data elements. The choice of ground metric plays a critical role in downstream tasks [135, 136] such as clustering, rransport-based inference, and dimensionality reduction. Many existing approaches rely on prior knowledge of the cost matrix structure or access to label information guide the learning process [137-139]. In many cases, however, labeled data is unavailable, making it necessary to develop unsupervised methods that rely solely on the data itself.
Unsupervised ground metric learning seeks to infer distance metric that captures meaningful pairwise relationships directly from the data, without relying on external labels or prior knowledge. A recent work in this direction is the Wasserstein Singular Vectors (WSV) method [59], which jointly computes Wasserstein distances over samples and features by using the Wasserstein distance matrix in one domain (e.g., samples) the ground metric for computing distances in the other (e.g.. features). However, the method is computationally demanding due to the repeated evaluation of high-dimensional Wasserstein distances. To address this issue, the recent work proposed Tre [60], which uses tree-based approximations of the Wasserstein distance to reduce computational cost. Both WSV and Tree-WSV adopt alternating scheme that iteratively updates distances across rows and columns, a strategy that aligns conceptually with our framework.
When working with hierarchical data, which is the central focus of our work, the Wasserstein metric is not inherently a hierarchical metric. As a result, methods such as WSV and Tree-WSY, which rely respectively on the Wasserstein metric and its tree-based approximation, are limited in their ability to represent hierarchical structures underlying samples and features. While Tree-WS WSV introduces tree metrics to improve computational efficiency, the trees constructed approximate Wasserstein distances, not to represent the hierarchical metric space.
In contrast, our method is explicitly designed to represent the hierarchical structures for both rows and columns. The key distinction lies in our use of the tree construction method [37] in the TWD: we construct the tree which employs hyperbolic embeddings and diffusion densities to reflect meaningful hierarchical relationships among data points. We then integrate this specific tree-based TWD into the construction of difffusinn operator, whose diffusion densities are used to generate hyperbolic embeddings that encode the joint structure of both samples and features. In doing so, we not only model each hierarchical structure but also capture their mutual influence through an altermating learning process.
Although WSV and Tree-WSY share similar alternating framework, their focus remains on metric approximation for ground metric learning rather than representation learning. contrast, our approach is centered on jointly learning hierarchical representations for samples and features and enhancing
37

---

them through wavelet filtering and integration within HGCNs. This provides additional insights into the hierarchical data that are not present in the WSV or Tree-WSV frameworks.
C
Proofss
We present the proofs supporting our theoretical analysis in Sec. 4. To aid in the proofs, we include several separately numbered propositions and lemmas.
Notation. Let Hm C CR C be of pairwise tree distance matrices for points. That is, given H E Hm. the following properties satisfy: (i) there is a weighted tree T = (V, E, A) such that H(j,)' dr (i-j') V j.J' [m)], (ii) H(j.) = 0 V [m]. (iii) (() S Jj) (2) V J2. Jo [m]] and (v) H(j.))) H(''.j) V j.j [m]. Let W Wm-an. C RT*= represent the set of pairwise OT' distance matrices for n points, where the ground metric is given by the tree distance matrix H Hmn. with points. That is, given W E Wm--nn the following properties satisfy: (i) W(i, = OT Pe: H) = ( Hv. T) for any two probability distributions
pig py Rm supported on the tree T. (ii) Wii,)) 0 V [n], (iii) W(i,, i2) S W(i,, ig)
W(1,, i2) for any probability distributions Mig. E R supported on the tree T, and (v)
W (i,,"") W(i, for any two probability distributions μi, pa IRm supported on the tree T.
C.1 Proof of Theorem
Theorem 1. The sequences w and w! all limit points are fixed points if "ro Tic 0.
generated by Alg. have at least one limit point, and
Proof. Consider the set S., (He E Hml Hal S 0g and H.(j.J") 0. Vj * j'), the Wasserstein distance between samples using a hierarchical ground metric H. is bounded by [59]
0. |r. 1l s OT(r,, ry, H.) S 2
fg ||:: rell1
2
where r, X |I. Similarly, consider S,, = (H, E Hall |H S 0, and Hi f') 2 0, vi # i'), the Wasserstein distance between features using H, as the hierarchical ground metric is bounded by
Cy H,) 0 èي cy < OT(ē,, cy. 2
0r é -crll1
where Cj x | 2l1. Therefore, the updated TWDs in Eq. (6), where we consider norm regularizer ((-), can respectively bounded by
0=
0
_e ||: relh + + -rv) S (( < - + +755(r - -re)
2
and
0, |c, - + volce- cyr) (..')
0, cyll -acc,,ycyll
|c,-
2
2
Moreover, from Prop. C.1 and Lemma C.2, the pairwise tree distances H, and H. constructed via Eq. (19) are bounded by the Wasserstein metries:
@oT(..H.) Ho(j.5") I HH)) and (........) Hi, S @oT(..H..)
Let
W Wm P, { (w!))
w!!)
I H.) and W(O(t,)) ēoT(..... H.) Vi vitt]
and
Pm = {w!) E W Wm..o w
S .o.(.....) and (i.') > HoT(.....) H.) Vj+j].
The product set
K. = (Sn x Sm) (Pn. x Pm))
38

---

is closed, bounded, and convex; hence compact in finite dimensions. Each update step of the iterative scheme remains in K:
w ob (XX;; ((W!))
w + (Xx:: T(W()))
contained in a compact set. We define the alternating step as the self-map
F(W,. We) = ((W)), 2(X (
Note that the tree decoder function T is marginal-separatee there exists e > 0 such that for every iteration and every distance pair (d1, (8,, 12)
W (11,1) w (1,1)) 2
this e-subset of K, the decoded binary tree is locally constant and thus continuous. Therefore, F is continuous K. and F is non-expansive [140]. Non-expansiveness, combined with bounded iterates, implies asymptotic regularity. By Opial's ' lemma, every cluster point a fixed point. Since the sequence is confined to the compact set IC, least one limit point exists. Therefore, every limit point of the iteration a fixed point, and the sequences generated by Alg. admit at least one convergent subsequence whose limit is fixed point of the alternating scheme.
Proposition C.1. Let |I | be norm Ih and let X = (x+ C TR be set of data points. Suppose the diffasion operator is constructed using a Gaussian kernel defined with respect II-IIR Then, for all + i' the embedding distance di (i) in Eg. (19) admits the bounds
@R S dd S QR.
(42)
where eR, @R Ryo are constants determined by the norm ||-1
Proof. Let CR. Or 0 be constants such that
CR Xx x x S CR.
Vi Y.
This ensures that the Gaussian affinity matrix
|x] - X 11
K(i,2') exp
bounded as
CR
K(i,t)) S exp )
cxp
Since the resulting graph is fully connected and the node degrees are uniformly bounded below (i.e., dig(i) 2 dmin), the continuous-tim timm diffusion operator Pt exp(-tL) satisfies the following heat kernel bound [141]:
Cr el
ps ('
exp
for some constant > 0, where L is the Laplacian matrix. Similarly, using the spectral decomposition of L and the exponential decay of eigenvalues, can upper bound
P (i < C2
CR
exp
et
for constant 02 0.
Using the bounds above, we obtain:
(v exp
|| [PP((.)) (P S S/n
(-cR)
CR
Cr el
exp
39

---

and, by Jensen's ' inequality
1/2
PP(,;:) :))
2 V2
11-1
1-
2
exp
C
Cl
exp
(-
C (-CR
c
ct
and similarly, an upper bound
Cr ))
P'((.:) - (' S
C2 n
exp
CR
Ci
exp
n
et
The embedding distance in the hyperbolic spaces is constructed from multiscale Hellinger-type distances of the form
K
dm(i,t) Σ 2inh
22
Vpz - (V,
k=0
Substituting into the definition of dM and use the monotonicity of we obtain:
dm(i,d') - Σ 2sinh-1 2 k/241
Vn. (V% exp
(-
CR et
91
Cr (-Cr (-el et
exp
=: eR-
For the lower bound, we use the inequality in Lemma C.1 and based the result [55, Lemma 3]:
1/2
di 2
Σ 2 k/2
K
2/22
/2 11 1211111
C2
exp
CR
Ci
exp
CR ct
)
=:
@R-
n
ct
n
k=0
These define constants eR and ΩR that depend only on the ambient norm and the Gaussian kernel.
Lemma C.1. Let p, q be two probability distributiinss on X. For a constant E Zzo. we have
2sinh -1 (2-k/2++ volle) 2 9l11.
(43)
where 1/p- - V4l12 is the unnormalized Hellinger distance between and Lemma C.2. Fix norm || || on T. and data cloud X {x,);=1 TR.. Let dg denote the tree metric output by Alg. 3. Then there exist positive constants U &r depending only on || l|re and the kernel bandwidth, such that
ØR < de(i,f) < d'R.
Vi
Thm. A. shows that the embedding distance in the hyperbolic spaces d is bi-l bi-Lipschitz equivalent to the intrinsic tree distance dr constructed from it: dr ~ d.v- Thm. A.2 in turn relates the decoded tree metric de produced by Alg. 3 dr, that dp dr. Combining these relations gives dæ = d.. Prop. C.I provides explicit upper and lower bounds on d, expressed through ØR and eR- Transferring those bounds through the bi-Lipschitz chain yields the stated inequalities for dg.
C.2 Proof of Theorem 2
Theorem 2. The sequences (!()) and W! ()) generated by Alg. 2 have at least one limit point, and all limit points are fixed points if Tr, Te O.
Proof. The discrete Haar transform is an orthonormal linear operator H. At each step the algorithm keeps only a subset of Haar coefficients chosen by an -based rule and applies soft-thresholding with thresholds 9_ and 0,. Soft-threshoding Sa(-) is 1-Lipschitz and non-expansive. Therefore each filter
V(:T(W
- So. p(0) H,
40

---

is a composition of operators whose spectral norms do not exceed 1, where p() is the projection (or masking) operator acting on the vector of Haar coefficients at the 1-th iteration. analogously for the feature mode. Therefore, for every 4,
V(X") ( (T(:) -
x 1
11
Because ( and ( finite, all subsequent filtered data and, through the convex updates that solve the regularized least-squares sub-prot uub-problems, all matrices w W remain in a ball of finite radius, which respectively form sets that is convex, closed, and bounded, therefore compact.
((W.") is firmly non-expansive, their composition is non-expansive, and it yields asymptotic regularity [140]. Together with boundedness, the asymptotic regularity activates Opial's lemma, which guarantees that every cluster point of the sequence is
Because each of ( (); T (((?") and
fixed point. The sequence ! and W! lies in the compact set. By Bolzano- Weierstrass it has at least convergent sequence.
D Incorporating learned hierarchical representation within HGCNs
We introduce an iterative learning scheme for joint hicrarchical representation learning based on TWD, where the initial pairwise distance matrices M, and M, are typically computed using standard data-driven metries, such as Euclidean distance the distance based on cosine similarity. A key advantage of our approach is its ability to incorporate prior structural knowledge when a hierarchical structure available for one of the modes. In such cases, the distance metric of the known hierarchy can be used to initialize the corresponding distance matrix, allowing the iterative scheme to refine the structure of the other mode. We demonstrate the applicability of this formulation to hierarchical graph data and show how it can be integrated into HGCNs [21, 45).
In hierarchical graph data, the input is modeled as hierarchical graph [n] .I.....71 represents the nodes, and each node is associated with an m-dimensional attribute vector. We denote the collection of node features as data matrix X E R=x=, consistent with the high-dimensional data matrix setting in our framework. In hierarchical graph data, the node hierarchy (corresponding the sample hierarchy in the data matrix) known, while the feature structure remains implicit. Recent works have proposed developing graph convolutional networks in hyperbolic space [21, 94, 45, 142], where node features are projected from Euclidean to hyperbolic space.
(mn), E, where
We show that the hierarchical representations learned by our method can be used as a preprocessing step for HGCNs. Specifically, first derive multi-scale hyperbolic embeddings from the data matrix X and use them as hyperbolic node features. During the neighborhood aggregation step, we replace the predetined hierarchical graph with the sample tree obtained from our iterative learning scheme after convergence for neighborhood inference. It enables hicrarchy-aware message passing, where the learned tree reflects structured relations among samples and features in a data-driven manner. We provide technical details of this integration below.
D.1 Hyperbolic feature transformation
It is important note that the used tree construction [37] decodes binary tree from multi-scale hyperbolic embeddings (see App. A). As result, after convergence, the produced hyperbolic embeddings can be directly used as initial hyperbolic node representations in HGCNs. For each node i, the multi-scale hyperbolic embedding is denoted by
(44)
where each y lies in a Poincaré half-space H forming a point in the product manifold of hyperbolic spaces H x... x H This representation serves as hyperbolic transformation of the original Euclidean node features and provides a coometry-awar initialization for downstream hyperbolic graph learning tasks.
The Poincaré half-space H=*+1 advantageous for representing the exponentially increasing scale of diffusion densities in hyperbolic space [50]. However, due to the computational advantages of
41

---

performing Riemannian operations in the Lorentz model IL", and the known equivalence between the Poincaré and Lorentz models with diffeomorphism [110, 143], we adopt the Lorentz model for feature transformation [21, 45]. The hyperbolic feature transformation is performed by first transforming points from the Poincaré half-space to the Lorentz model, followed by applying the Lorentz linear transformation [45] and mapping back to the Poincaré half-space.
To enable efficient hyperbolic feature transformation, we first map each point y from the Poincaré half-space to the Lorentz model using
y = P(y!) = (1 15 (.,, 1))
(45)
||
where
5 24 (n),
1 + 1)
The hyperbolic feature transformation is then applied via Lorentz linear map [45]
1, 0 O, j
h = J s.t. J.
and JJ
(46)
This transformation preserves the hyperbolic geometry (i.e., manifold-preserving)) and corresponds to linear feature map in hyperbolic space. After the transformation, the features h are projected back to the Poincaré half-space, completing the integration of our learned embeddings into HGCNs.
D.2 Hyperbolic neighborhood aggregation
Since the hyperbolic feature transformation is manifold-presss preserving, the neighborhood aggre- gation in hyperbolic space as a weighted Riemannian mean neighboring points, followed by a non-linear activation function. Specifically, during the aggregation step HGCNs, replace the predefined hierarchical graph H with the sample tree inferred by our method after convergence for neighborhood inference. This learned tree serves hierarchy-aware structure that guides message passing, capturing data-driven relations among samples and across features.
After mapping the transformed features back to the Poincaré half-space, the hyperbolic neighborhood aggregation for node at layer l and scale k is given by
(wish(o))) SETAII
Σ
(47)
where [[i]] denotes the neighborhood of node i, was aggregation weights, and is non-linear activation function. In contrast to the Lorentz model formulation [21, 45], applying non-linea activations in the Poincaré half-space does violate manifold constraints, the operations remain in the same space. Notably, the weighted sum within the activation function corresponds to weighted Riemannian mean in the Poincaré half-space model.
The integrated hyperbolic feature transformation and aggregation are conducted directly on the manifold [45]. This contrasts with the earlier method [21]. which relies projecting points to the tangent space for linear operations and then mapping back. By operating entirely within the hyperbolic manifold, integrated method maintains geometric fidelity throughout the learning process in HGCNs.
E Details on experimental study
We provide a detailed description of the experimental setups and additional implementation details for the experimental study in Sec. 6. The experiments conducted NVIDIA DGX A100.
E.1 Datasets
We evaluate our methods on word-document data, single-cell RNA-sequencing (scRNA-seq data), and hierarchical graph datasets. We report the details of the dataset statistics in Tab. 4.
42

---

Table 4: Summary of dataset statistics for word-document, single-cell RNA-sequeneinn (scRNA-seq),
and hierarchical graph benchmarks.
Dataset
# Samples / # Nodes
# Classes
# Edges
# Features
BBCSPORT
517
13.243 BOW
TWTTTER
2,176
3
6,344 BOW
CLASSIC
4,965
24,277 BOW
AMAZON
5,600
4
42,063 BOW
ZHISEL
3,005
47
4,000 Genes
CBMC
8,617
56
500 Genes
DISEASE
1,044
2
1,043
1,000
AIRPORT
3,188
4
18,631
4
PUBMED
19,717
88,651
500
CORA
2,708
7
5,429
1,433
Word-document benchmarks. We evaluate method four standard word-document bench- marks commonly used in the word mover's distance [25] and TWD literature [34-37]: (i) the BBCSPORT dataset, consisting of 13,243 bags of words (BOW) and 517 articles categorized into five sports types, (ii) the TWITTER dataset, comprising 6,344 BOW and 2, 2,676 tweets classified into three types of sentiment, (iii) the CLASSIC dataset, including 24,277 BOW and 4,965 academic papers from four publishers, and (iv) the AMAZON dataset, containing 42,063 BOW and 5,600 reviews of four products. The pre-trained Word2Vec embeddings [144] are considered as the word embedding vector, trained on the Google News dataset, which includes approximately 3 million words and phrases. Vee represents these words and phrases vectors in 30000 The document types serve as labels for classification tasks in Sec. 6.2.
Single-cell RNA -equeeqqunni data. Two scRNA -eqq datasets [84] are considered:(i) the ZEISEL
dataset: From the mouse cortex and hippocampus [145], comprising 4,000 gene markers and 3,005 single cells, and (ii) the CBMC dataset: from a cord blood mononuclear cell study [146], consisting of 500 gene markers and 8,617 single cells. We used the divisive biclustering method [145] to obtain 47 classes for Zeisel and 56 classes for the CBMC. The Gene2Vec [147] is used as the gene embedding vectors [2]. Cell types used classification labels in the classification experiments presented in Sec. 6.2.
Hierarchical graph datasets. For hierarchical graph data, we consider the following datasets: (i) the DISEASE dataset: constructed by simulating the SIR disease spreading model [78], where node labels indicate infection status and node features indicate susceptibility to the disease [21], (ii) the AIRPORT dataset: flight network dataset where nodes represent airports and edges represent airline routes. The label of node indicates the population of the country where the airport is located [21]. (iii) the CORA dataset: citation network containing 2,708 nodes, 5,429 edges, and 1,433 features per node, with papers classified into seven machine learning categories [9], and (iv) the PUBMED dataset: Another citation network with 19,717 nodes, 44.338 edges, and 500 features per node, encompassing three classes of medicine publications [9]. The node labels are used for node classification experiments in Sec. 6.3.
E.2 Hyperparameters settings
We describe the components required to initialize and configure our method, including the choice of initial distance metrics, the hyperparameter configuration, and the norm regularization.
Initial distance metric. The initial pairwise distance matrices M, E Raxn and M. E mm are typically data-driven and constructed from the observations using standard similarity measures such as Euclidean or cosine distances [148]. These initial distances are then used to define Gaus- sian kernel for the diffusion operator, as detailed in App. A. For word-document and scRNA-seq datasets, the initial distance matrices, M. E mm (feature distances) and M, E R (sam- ple distances), are computed using cosine similarity in the ambient space, following prior work [149, 150]. Empirically, cosine-based distances yield more stable and reliable tree representations
43

---

compared to using Euclidean distances as initial distance, which were found to be less robust and led to decreased performance in our experiments. While different initial distance metrics lead to different tree representations, we observe that the resulting trees from Alg. and Alg. 2 converge to unique representations when the regularization parameters "rr and . are sufficiently large (see App. F). However, we emphasize that such large regularization does not necessarily produce the most meaningful hierarchical representations. When "/r and To are too large, the regularization terms dominate the distance computation in the iterative learning scheme, leading to trees that may lack meaningful hierarchical structure.
Hyperparameter configuration. We explore broad range of hyperparameters accommodate variations across tasks. For hyperparameter optimization, we conduct a grid search for each dataset using Optuna [151]. The Gaussian kernel scale is selected from the set {0.1, 1, 2, 5, 10} X X. where X denotes the median of the pairwise distances [152, 153J. The number of hyperbolic components in the product manifold set with the range K E {0. 1, ....,.................. 19). The regularization parameters "/r and Te are chosen from the set { (10-5 3, 5, x 10 10,00,22 x 10 10-1,5 x 10-1, ,,,, 5, 10 5 x 10 10 The thresholds Ug and 0, the filtering step are set as follows. At the first iteration the total L1 norms of the Haar coefficients across all samples and all features are computed, de- noted by T7. and Th., respectively. Thresholds 1), and 1,, are then respectively selected from the sets [0.1n,, 0.2 ..., 0.9ne) and (0.1 0. 0.9nr) As the optimal hyperparameter configura- tion varies across tasks and there is no universal choice that works for all cases, Optuna [151] is used to systematically explore the parameter space and efficiently identify task-specific configurations.
0).
Norm regularizer.
We adopt a norm regularizer based on the snowflake penalty [57], given by
<(1: - r--== 1 I: 1ll
(48)
VE 10-6
This smoothness- prreeesinppoooingggggo regularizer increases monotonically over [0, oo) and imposes stronger penalties on small differences, thereby encouraging local regularity in the representations. While alternative regularization strategies, such as entropic regularization 19], could be considered, we demonstrate in App. F that they are less effective.
E.3 Scaling to large datasets
In the tree construction [37]. one key step involves building diffusion operators [51] between both samples and features. Although such construction can be computationally intensive, recent devel- opments in diffusion geometry literature have introduced various techniques to significantly reduce the runtime and space complexity of this process. In particular, one can use the diffusion landmark approach [76], which enhances scalability by reducing the complexity from O(n²) to O( 122)) for samples and from O(m²) Omm2+2c) for features, where e <<1 represents the proportional size of the landmark set relative to the original dataset. We briefly describe the diffusion landmark approach below for the setup of the sample diffusion operator and note that it can be seamlessly applied the feature diffusing operator.
Given a set of samples X' = [x; S R in an ambient space R", let = = (x]] C C R be the landmark set of A, where X' and the original set is denoted by = exp( - /e), where M(i, i') represents a suitable distance between the sample x; E X' and the sample x E X', and e > 0 is the scale parameter. Let D be diagonal matrix, where (i,, 8. KK 1,, and 1, [l, ..... 1] E Ra. The landmark- matrix Y KK T R has an eigen-structure similar to the diffusion operator, which can be computed by applying SVD to the matrix D - 11/2 = V. To construct the diffusion operator on X. one can use the landmark set X" and its eigenvectors.
log, (n') 1. The affinity matrix between the landmark set
In our method, for large datasets, we fix the landmark set size to = 0.1, which reduces the computational complexity of tree construction at each iteration to O for samples and 0mm!2) for features. Additionally, the computation of the TWD is linear, resulting in an overall per-iteration complexity of 1 + -22)). In contrast, a naive approach without diffusion landmarks leads to significantly higher complexity: O(mre) for tree construction and O(n² log n) and O(m log m) for computing Wasserstein distances over rows and columns, respectively. The total per-iteration cost
44

---

becomes Omm?² + nm3 + n² log + m² logm). As shown in Sec. 6, our method scales efficiently and can handle datasets with thousands of samples and features.
E.4 Experimental setup
We describe the experimental setups for sparse approximation to evaluate the hierarchical representa- tions, for document and single-cell classification using the learned TWD, and for link prediction and node classification on hierarchical graph data.
Sparse approximation setup. We assess the learned hierarchical representations using sparse approximation criterion. Specifically, we compute the L, norm of the Haar coefficients across all samples and features. Let B, and B, denote the Haar bases induced by the feature tree and sample tree after convergence, respectively. We measure sparsity by computing:
12 ΣΣ (c), GT),
1 S ΣΣ (2,.B.)
(49)
m
These quantities reflect the average sparsity of the transformed data in the respective Haar bases. A lower L1 norm indicates more compact representation, where the signal is concentrated on fewer significant coefficients. This suggests that the learned tree structures more meaningfully represent the hierarchical organization of the data [43, 53].
Document and single-cell classification setup. For document and single-cell classification tasks (i.e.. sample classification). we use the sample TWD matrix to perform classification. The dataset split by partitioning the TWD matrix into 70% training samples and 30% testing samples. We apply a k-nearest neighbors classifier with & E (1.3.... 19}. The random split is repeated five times [25, 26), and we report the best average classification accuracy across these runs.
Link prediction and node classification setup. We follow experimental setups [21,45] to maintain consistency. For the link prediction task, the set of edges is split into 85% for training, 5% for validation, and 10% for testing. A Fermi-Dirac decoder [93, 13] is used to compute probability scores for edges, and the model is trained using enosopy loss with negative sampling. The link prediction performance is assessed by the area under the ROC curve (AUC). For the node classification task, dataset-specific splits are applied [21, 45]: 70%/15% 15% for AIRPORT, 30% 10%/60% for DISEASE, and the standard setup of 20 training examples per class for CORA and PUBMED. The node classification is performed using a centroid- based method [94], where each class is associated with prototype and predictions are made using a softmax classifier trained with cross- coos-py loss. Additionally, an LP-based regularization objective is incorporated during training [21, 45J. The NC task is evaluated using the Fl score for binary-class datasets and accuracy for multi-class datasets: F1 score for DISEASE, and accuracy for AIRPORT, CORA, and PUBMED.
F Additional experimental results
We present additional experimental results supporting our method, including toy problem, ablation studies, empirical convergence, runtime analysis, and co-clustering performance.
F. Toy example: synthetic video recommendation system
We consider toy example motivated by a video recommendation system, where samples correspond users and features correspond videos. Both users and videos are assumed follow hierarchical structures, and representations are generated through a probabilistic diffusion process over these trees. The video hierarchy is based on content categories. The root node represents all videos and branches into three main types: fiction, documentary, and animation, as illustrated in Fig. 4 (right). Each category further divides into subgenres. For example, fiction splits into action, drama, and sci-fi; documentary into biography and historical documentary: and animation into family and comedy. We generate video embeddings [wr) IRd using a hierarchical probabilistic model. The root node embedding is drawn from a Gaussian prior, and embeddings for child nodes are recursively sampled from Gaussian centered at their parent's is embedding:
yoot N0 ()) welw, N N(wp,o.11
45

---

Molia
Doo Dackmboory
Geaviiiuuus Twrding Serily Suleing- kanal
00d
Dar
the
the
Die
Lier
Vidoo
Viloo
valoo
Vilo
Valoo
Valoo
Moo 1X
Vido 14
Vioo
Vido
Figure 4: Hierarchical structure used the toy video recommendation example. The right tree depicts the user hierarchy based device type and viewing context. The left tree represents the feature hierarchy of videos, branching by genre and subgenre (e.g., fiction - action, drama, sci-fi). Users and videos are colored according to their first-level subcategory.
Alg.1
Alg.
User
Video
User
Video
@bokkp auTs
Sel Ing-farna
D.oaing Basty Chuad
Haray
Cieedy
Actes
Figure 5: Learned tree representations for the toy video recommendation example using Alg. 1 and Alg. 2. Users and videos colored according to their first-level subcategory.
This recursive diffusion continues until all leaf nodes (i.e., individual videos) receive their embeddings. As result, videos belonging to the same (sub)category remain close in the space, reflecting their semantic proximity. A similar process is used to generate user embeddings {s, ) C Rd based on a hierarchical structure defined by user device and consumption context, shown Fig. 4 (left). The user embeddings are sampled as:
Statt (O ofI), Nspp: I)
Given user embeddings {s and movie embeddings (w, w the interaction matrix Y is defined based their proximity:
Yu = ||s w1ll2 + Esj. where N N o ).
Note that other functions of the distance can also be considered, e.g., clastic potential operator [154]. In our synthetic experiment, we set a. = 1, ar = 0.6, op = 0.5, 01 = 0.25, muise- and d = 30. Finally, to simulate an unstructured observation, we apply random permutations [129, 54] to the rows and columns of to obtain the observed data matrix X.
Fig. 5 presents the learned hierarchical representations for users and videos in the example of toy video recommendation. The trees are obtained by Alg. and Alg. 2 after convergence. Leaves are colored according to their first-level subcategory labels. The hierarchical representations learned by Alg. 2 exhibit a clearer structure and a more distinct separation of subcategories, closely aligning with the ground-truie hierarchies shown in Fig. 4. It highlights the benefit of incorporating wavelet filters into the joint hierarchical representation learning process. This visual evidence supports the effectiveness of filtering in learning meaningful hierarchical representations, We remark that in the absence of noise, we observe that both algorithms produce comparable results. In the presence of noise, the tree learned by Alg. deviates slightly from the ground-trun -tuth hierarchy, resulting in less accurate representation (see arrow in Fig. 5). In contrast, Alg. 2 effectively attenuates the noise, accurate and coherent hierarchical structure, shown in Fig.
leading to
Fig. 6 shows the trees produced when alternative TWD constructions are integrated into our iterative learning scheme. In contrast to the clear hierarchical structures produced by approach, these
46

---

Co-Quadtree & Co-WQTWD & Co-SWQIWD
Co-WCTWD & Co-SWCTWD
User
Video
User
Vidco
Figure 6: Learned tree representations for the toy video recommendation example using alternative TWD methods in our iterative learning scheme. Users and videos are colored according to their first-level subcategory.
trees fail to represent the hicrarchical structures the samples and the features. This result further illustrates that when using alternative TWD methods in the iterative learning scheme, their trees designed for hierarchical representation learning and therefore produce disorganized structures that not meaningful for visual comparison. In addition, present quantitative results in Tab. 5 below, reporting the performance on sparse approximation and classification accuracy. The experimental setup follows that of Sec. 6. For the classification tasks, we use the first-level subcategory labels as the ground truth labels. Tab. 5 shows that our method consistently outperforms baselines, validating its effectiveness in representing meaningful hierarchical structures.
Table 5: Quantitative results for the toy example. As in Sec. 6, sparse approximation is measured by the L1 norm of the Haar coefficient expansion, reported the format samples features (lowest values in bold, second-lowest values underlined). Classification accuracy for users and vidcos is reported in the same format, with the highest values in bold and the second-highen highest underlined.
Sparse Approximation
Classification Accuracy
Co-Quadtree Co-Flowtree
9.2
90.1 88, 4+0.4
9.1
89.2.1077 6:1..5
Co-WCTWD
933
82.6:110/ 86.410.8
Co- CO-WQWWD
/9.7
83.1+0.8 85.9=0.7
Co-Ultrarree
8,9/9.0
88.51144 82.310.9
Co-TSWD-1
10.3/114
85.61114 88.2211.1
Co-TSWD-5
1
87.240.88 01.+-088/449400. 84,9+0.9
Co-TSWD-10
9.7
89.7.00. 87.4.0.8
Co-SWCTWD
8.2889
894410/ 88.310.4
Co-SWQTWD
8.4/866
88.6+0.2/ 88.4+0.9
Co-MST-TWD
1
82.1.1.77 8..6.1
Co-TR-TWD
9.8/
92.6.10.4 / 810.7
Co-HHC-TWD
103/114
89.2+-2.0 88 1.7
Co-gIIIIC-
10.1/ 10.9
88.4....6 87.0.1
To-UIraFie-WD
9.7 10.3
90.1.1.3 89 19.22
QUE
8.6
93.4+ 92.5+00.5
WSV
92.5.10.8 91.6..0.6
Tree- WSV
8.4/8. 8.7
93.910.66 92.710.88
Alg. I
3.2/4.0
0.3 95.
Alg. 2
2.5/3.1 3.1
99.4+0.2/ 99.6+0.11
F.2 Empirical uniqueness of hierarchical representations under strong regularization
For sufficiently large "/ and Te. the resulting TWDs are unique. This means that even if the initial distances are random, our methods will converge to unique trees and TWDs. Fig. 7 illustrates that using different random initializations yields identical final expansion coefficients for the Zeisel dataset. However, we would like to point out that such large regularization does not necessarily lead to the best hierarchical representation. When 7, and Ne are too large, the regularization terms dominate the distance computation in the iterative learning scheme, leading to trees that may lack meaningful hierarchical structure. This is evident when comparing these results to those in Fig. Sec. 6, where moderate regularization (tuned by Optuna [151]) yields more informative trees. The purpose of this experiment therefore to demonstrate empirical uniqueness under strong regularization, rather than to advocate for such settings in practice.
47

---

Zeisel Sample Tree
Zeisel Feature Tree
250
280
Alg.1 rd
Alg.1
225
Alg.1 rd
260
Alg.1 rd 3
200
Alg.1 rd
240
Alg.1 rd
175
rd r
220
Alg.2
Alg.2 rd
Alg.2
150
Alg.2 rd
200
Alg.2 rd
125
180
160
100
140
120
10 15
20
5
10
20
# of iteration
# of iteration
Figure 7: Illustration of empirical uniqueness under strong regularization for ZEISEL datasets. The L, norm of the Haar coefficients from the sample tree and the feature tree during the sparse approximation task across iterations the ZEISEL dataset. Despite different random distance initializations, the final expansion coefficients are identical when using sufficiently large 7/r and Ye. indicating convergence unique hierarchical representations.
F,3
Ablation study: effectiveness of iterative learning scheme in classification tasks
To demonstrate the effectiveness of our iterative procedure in improving the joint hierarchical representation learning for both samples and features, we evaluate its impact on the learned TWDs in downstream classification tasks. Specifically, and single-cell datasets using both the initial TWDs, i.e. computed without any iterative interaction between the sample and feature trees, and the TWDs obtained after applying our full iterative learning scheme. Tab. 6 reports the classification results using the initial TWDs, serving as a non-iterative baseline, alongside the results from iterative method reported in Tab. 2. The comparison shows a consistent improvement in classification accuracy when using the TWDs refined through iterative process. This suggests that alternating updates between sample and feature trees help better inference for the data, regardless of the specific TWD method used. Importantly, our proposed methods Alg. 1 and Alg. not only improve the quality of the learned TWDs but also achieve better performance than all baseline methods. These results highlight the advantage of our approach in jointly refining hierarchical structures across data modes, validating the effectiveness of the proposed iterative learning framework for hierarchical data.
compare the classification performance
document
F4
Ablation study: incorporating Wavelet filtering with alternative TWD methods
The key difference between Alg. I and Alg. 2 is the application of filtering step in each iteration of Alg. which aims suppress noise and other nuisance components. This raises the question of whether wavelet filtering could also benefit other TWD-based baseline methods when considering in our iterative learning scheme. To explore this, we apply Haar wavelet filtering using trees constructed from various TWD baselines and report the resulting L1 norms of the Haar coefficients across all samples and features in Tab. For each baseline, we append ". "-Wveelet" indicate the variant with wavelet filtering applied. We observe that incorporating wavelet filtering into these baselines does not consistently lead to improved hierarchical representations. We argue that this is because the trees used in these baselines are primarily optimized to approximate Wasserstein distances as ground metries, rather than to represent the hicrarchical structure of the data. As such, they do not serve as meaningful hierarchical representations, since the Wasserstein metric is inherently not hierarchical metric. Consequently, applying wavelet filters, whose effectiveness depends on the data and the tree structure, does not improve, and may even hinder, hierarchical representation learning.
We further evaluate the effect of wavelet filtering on downstream classification tasks for document and single-cell datasets. Tab. 8 reports the classification accuracies obtained by applying wavelet filtering within the iterative learning scheme using alternative TWD methods. Similarly, for each baseline, we append "- "Wavelet" to indicate the variant with wavelet filtering applied. We see that the results show inconsistent improvements across datasets and domains, with any gains being marginal at best. In contrast, when wavelet filtering is applied within proposed method, the improvement
48

---

Table 6: Document and single-cell classification accuracy using the initial TWD (without iterative refinement) and the updated TWD (with iterative refinement). Results marked with * are taken from the work of Lin et al. [37]. Arrows (1) indicate improvements achieved through the iterative learning process. The best performance is marked in bold, and the second-best is underlined.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Quadtree
95.5+0.5*
69.6+0.8*
95,9++0.4
89.3++.3'
80.1+1.2"
80.6+0066
Co-Quadtree
96.2:00.4 (1)
69.620.3
95.910.2
89.4=0.2
81.7+100
80.7=0.3
Flowtree
95.3+1.1
70,2+0.9*
94.4+0.6'
90.1+0.3'
81.7+009"
81,8+0,9
Co-Flowtree
95.7100 (f)
71.540.7 (t)
95.620.5 (1)
914400 (1)
84.340.7 (1)
83.0 (1)
WCTWD
92,6=2.1*
69.1+2.6*
93.7+2.9"
18,2=14'
81.344.9
78.443.3'
Co-WCTWD
93.2.11.2 (f)
70.2.12.1 (f)
94.7.2.6 (1)
87,411.0
82.5119 (t)
79.412.1 (1)
WQTWD
94.321.7*
69,4122,4
94,613.22
87,4
80.913.5*
79.113.0"
Co- Co-WQTWD
9..7+ (1)
70.7+22.
1.3
88.2- (1)
82.3+3.1
80.5+ (1)
UltraTree
93.111..5
68.113.2*
92.311.9
86.213.11
83.911..'
82.311.6'
Co-UhraTree
9..3+1.4
70.1+288
93.6+2.0
86.5+288 (1)
85.8+1.1
1.3 (1)
TSWD-1
87.6+1.9
69.8+1.3*
94.5=+0.*
85.5=0.6'
79.6+18"
72.6+18'
Co-TSWD-1
88.2.1144
(0.41..2
94.7_0.9
86.110.5
80.2.11.4)
73.211.00
TSWD-5
88.1+1.3
70,5+1.1
95,940.4
90.8-0.1"
81,3+144
74,911.1
Co-TSWD-5
88.711.77 (1)
71.011.5
9..7
91.5
82.010.9
75.410.7
TSWD-10
88.610.9*
70.7+11.
95,911.66
91.110,5*
83,21088
76,510..7
Co-TSWD-10
89.2+1.1 (1)
71)++18
95.5+0.2
91.8-0.77 (1)
8388+0.5
77.2+0.9 (1)
SWCTWD
92811..
70.211.2*
94.1
90.2.1..2*
81.913.1"
78.3
Co-SWCTWD
(3,5+2.4
1..
94,4=1.3 (1)
90.7=1.5
82.7+1.7 (t)
79.0=0.9 (1)
SWQTWD 94,5+110"
7..6+1.9"
95.4++2.0*
1.18+1.1'
80.7+2..5
7..8+2.5'
Co-SWQrWD
96.2+1.2
72.422.1
96.0=1111 (1)
90.6=2.3
82.4+144 (1)
81.3+111 (1)
MST-TWD
88.411.'
68.2+1.9
90.0+3.1'
86,421.2'
80.143.1"
76.24225.
Co-MST-TWD
88.7+224
68.4+33.3 (1)
91.3+ (1)
87.1=114
80.1+:2.8
3..5+
TR-TWD
89.210.9*
70.2.1077
92.910.8*
88.711.1
80.310.7
78411..2
Co-TR-TWD
8..5+1.2
70.9+1.7 (t)
93.4-222
89.5=1.4
80.7+0.8 (1)
78.5+0.9
HHC-TWD
8.8
70.410.4*
93.4:10.8
88.5
82.310.7
77.31111
Co-HHC-TWD
86.1+2.1 (1)
93.6=1.5 (1)
(1)
88.5+0.5
83,2+114
77.6=0.8 (1)
gHHC-TWD
83.242.4
69.9+1.8
90.3+2.22
86.9+2200*
79.4+199"
73.6+166
Co-gHHC-TWD
84.0+2.0 (t)
70.4+1.6
90.7=1.7 (t)
87.2=199 (1)
79.9+14
84.2+1.2
UhrDFit-TW
84.9+1.4"
69.5+1.2"
91.6+009*
87,4+1.6'
81.9+3.3"
77.8+1.2
Co-UlraFic-
86.8+0.9 (1)
1..9+1.
91.9+100 (1)
2..
83.7429
79.1+188
WMD
95.410.7*
71.310.6*
97.210.1*
92.610.3*
GMD
84.2+0.7"
81,4+007'
HD-TWD
96.1+
73.4+0.2*
96.9++..2
93.1+0.4*
89.1+
84.3++0.3*
Alg. 1
96.7+0033 (t)
74.1
97.3=02
94.0+0.4
90.1
86.7=0.5
Alg. 2
97.310.5
76.7.1077
97.6
94.2_0.2
94.0.10.6
93.310.7
is significantly more pronounced. These findings are consistent with the sparse approximation results in Tab. reinforcing the observation that wavelet filtering effectively improves hierarchical representation learning only when applied to trees that meaningfully represent the data hierarchy.
F.5 Ablation study: integrating wavelet filters to HGCNs without iterative learning scheme
In Sec. we demonstrate that our method can serve as a preprocessing step for hierarchical graph data and be integrated into HGCNs. The implementation details are provided in App. D. Here, examine whether wavelet filtering alone, without our iterative learning scheme, can also improve the performance of HGCNs.
important to note that Haar wavelet filters are induced by a tree structure [43]. However, in the case of hierarchical graph data, the structure is not always a tree. To quantify how closely graph resembles a tree, one can use the notion of 5-hyperbolicity [107]. tree has whereas hierarchical graphs that deviate from tree-like geometry exhibit larger d d-hyperbolicity values,
49

---

Table 7: The L1 norm of the Haar coefficients across all samples and features when applying wavelet filtering induced by trees constructed from various TWD-based baselines. Lower values indicate sparser representations. Arrows (1) denote cases where wavelet filtering leads to improvement (i.e., reduction in L, norm) compared to the unfiltered baseline. Values are reported in the format samples/ fealures (the lowest in bold and the second lowest underlined).
BBCSPORT
TWITTER
CLASSIC
AMAZON 87.3
ZEISEL
CBMC
So-Quddrree
27.6
64.3/ 108.7
1616.6/
Quaddree-aarccee
29.4/29.7
75.8/35.4
65.9 120.8
(1) (1)/0663
173.9/214.7 (T)
1792.4/ 89.2
Co-Flowtrce
27,/31.
69.1/209
73.7
/
173.7 /2000
1688.8 81.7
Co-Flowtreo-Wavee
26.3 31,9
62.9/22.7
74.8 /100
(7) 109.2
180,4 245.1
1537 (T) 79955
Co-CCTWD
26.5/26.7
57.6/77.
63.3/81.7
77.4/964
232.4/20299
1308.1 66
Cs-WCTWD-WWee
28.1/28.4
58.5/2077
67.1 /84.3
80.97
139.6/211.7
1412.3/794
Co-WTTWD
25.2/32.1
56.97 30.2
61.37 100.
74.7/ 111.0
135.0/ 224.8
1324.8 1
Co-WQrWD-Wan
27.1/40.3
33.1
68.71 121.3
80.9 0 114.7
123.2/
1417.5 /7.6
Co-UltnyTree
37.6/32.1
69.8/ 18.8
76.0/ 125.5
86.3 133.3
155.4/ /26.
1450.0
Co-UlenTTree
42.5/28.3 <ft)
75.1/25.9
/72.3/257.8
1524.6
Co-TSWD-1
26.0
70.17
69.8/916
100.4/1112
179.5/2115
1716.1 79.5
Co W-II-SSWD--aaas
29.7/34.5
76.2/ (4)8
73.9/98.4
103.6/ 129.1
/239.4
1788.5 /94.
Co-TSWI-5
305/24.7
58./ 160
65.7/9%0
X3.3 102.3
170.2 234.7
1560.2/
Co-TSWD-S-Naveler
33.1/28.0
64.2/ 19.1
66.7
90.5
196.7
1798.4/93.1
Co-TSWD-10
21.1/34.3
22.6/23.6
59.37 14088
74.5 135.8
4414//2997
1373.3/81.2
Co NO-TSW--JO-Waael
24.5/37.9
23.0/
63.47 151.2
73.9 /1/9 114..7
153.6/2324 /232.4
1425.6/ 79.2 (1)
Co-SWCTWD
35,0/29.5
24.7
79.8/93.
95.9/ 104.6
170,4/215.9
1595.07 89.7
CoCTWC
36.1/314
72.4/33.9
84.6
124.1
(1)/
1723.4794.6
Co-SWOTWD
33.5/25.6
$7.671399
61.3/36.0
77.0/98.
141.8 /09066
/
Co-SHQIWD-WaR Wavelet
34.9/27.7
63.4/ 14.8
62.7/92.4
86,5 / 106.7
119///203. (1)
1428.4/ 70.3
CST-ST
30.4/4411
69.1/22.1
77.51 109.3
97.0/ 126.4
179.1/254.0
1755.3/83.7
Co-MST-TWD--R
73.2/26.1
87.17
105.6/ 139.7
192.8 276.3
1863.7/94
Co-TR-TWD
36.8/21.5
19.1
66.1/91.7
82.2 102.3
121.6/ 238.6
926.8 68.8
Co-TR-TWD-Wavcelt
24/2(()/25.6
59.3/ 19.2
192.7 (*)
84.6 105.7
139.6 274.8
1031.9 72,3
Co-HHC-TWD
22.5/22.7
11.5/136
58.7/
74.1 110.1
1
1148.8 79.0
Co-HHCCTWD--wwele
28.7/36.1
$3A/196
64.2
126.8
276.4
1253.6
Co-gHHC-TW
27.9710.6
650/167
72.8 72.871122
88.5 715.5
162.7 240.7
1612.8
Co-gHHC-TWD-WWL Wieeeeeelet
32.6/16.1
78.2/
84.1 1225.
96.0 139.7
182.6 237.4
1823.5/ /
Co-UlrsRe-TW
22.5/22.1
51.9/134
58.4/813
73.1 / 19288
133.7 202.0
1294.6 78.
Co-UIraPs -WD-Waveett
26A/29.3
58.2/ 19.6
64.3/89.2
104.8
159.4/ /41.0
1432.5 89.0
QUE
22.8/19.8
57.87 10.4
71.6/67.6
88.8 72.0
93.1 173.3
906, 58.5
QUE
23.1/17.6.0
56.9 (6)9 13.7
74,2 69.4
90.6 78.5
94.5 162.7
802 4 (T) 69.0
Tre- WSV
24.6//2.9
18.2
65,4/99.2
84.2 106
/
16377.4 /
Weeeeee
27.2//2911
63A/22.5
67.2/
1204
145.7
1739.2
Alp Alg. 2
1..9/7.0
25.4337
40.4/246
$1.07
64.4/13088
$11.37
10.1 (1)/4((())
22.4 (1)/3.44
37)2 (9,4 (t)
46.6 (1) 26.6 (?)
50.9 (7))97777 99.7 (t)
489.4 (?) 45.6 (1)
indicating lower tree-likeness or increased curvature. For the hierarchical graph datasets considered in our experiments, following prior work on HGCNs [21, 45], the measured o-hyperbolicity values are as follows: DISEASE (6 = 0), AIRPORT (8 = 1), PUBMED (6 = 3.5), and CORA (Ó = 11). Among them, only DISEASE forms tree, allowing the Haar wavelet to be applied directly. For the other datasets, which are not trees, we can either convert the graph into a tree before applying Haar wavelets or construct Haar-like wavelets using the multiscale tree approximation framework [43].
For the naive transformation, we consider two strategies to approximate the graph with a tree structure: (i) we construct an MST by connecting all nodes using edges with the smallest cumulative weights; (ii) we apply agglomerative clustering to iteratively merge nodes and form a hierarchical tree, as in standard hierarchical clustering. These transformations alter the original graph topology, but the resulting tree is intended to closely approximate the structure of the original graph. Once the tree constructed, it can be used to define a Haar wavelet basis and build the corresponding wavelet filters, as described in Sec. 4.2, which are then applied to the node attributes (i.c., features). To construct Haar-like wavelets using the multiscale tree approximation framework [43], we first build hierarchical partition tree over the data using suitable distance metric (e.g., shortest path distance). Each wavelet basis function then defined over a node in the tree by assigning opposite-sipe signed, normalized values to its children. We denote wavelet filtering applied to the node attributes (i.e., features) using these approximations, without our iterative learning scheme, by "-MST" (for MST- based tree), "-HC" (for hierarchical clustering), and "-Wavelet" (for multiscale tree approximation). We incorporate these variants into HGCNs [21, 45) and evaluate their performance on link prediction and node classification tasks.
50

---

Table 8: Classification accuracy document and single-cell datasets using various TWD-based baselines with and without wavelet filtering, applied within the iterative learning scheme. The Wavelet" suffix denotes the inclusion of Haar wavelet filtering. Arrows (f) denote cases where wavelet filtering leads to an improvement compared to the unfiltered baseline. While some baselines show marginal improvement, the gains are inconsistent across datasets, highlighting the importance of using meaningful hierarchical structures for wavelet filtering to enhance representation learning. The best performance is marked in bold, and the second- underlined.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Co-Quadtree
96.21:0.4
69.620.3
95,910.2
89,4-0.2
81.7=10
80.720.3
Co So-QQuidree-- Quaveree-Wavelet Wavelet
95.6:10.6
69.6:00.2
96.040.1
89.4=110
81.4=1.1
80.6+0.4 (1)
Co-Flowtree
95.7100.9
71.510.7
95.61:0.5
91.4=0.4
84.3=0.7
83.02112
So-FFlottee-Waam Wavelet
95,5++0.8
71.0+ 0.8
95.0++0.7
90,8=005
82.7=0.88
82.2+11.0
Co- WCTWD
93.2++1.2
70.2+2.11
94.7++2.6
87.4+1.0
82.5+2.9
79.4++2.1
Co-WCTWD-Wawde
93.2411.1
70.6::19 (1)
94,5+:2.2
88.3= (1)
82.2=3.0
79.0:2.99
Co-WQTWD
95.7++1.8
70.7+-2.2
95.5++1.3
88.2+2.1
82.3+3.1
80.5+2288
Co-WQTWD-Wayee
95.0:11.6
70.2:1..8
95.24 1.8
88.24119
81.9=3.0
80.7+2.7 (1)
Co-UltraTree
1..3
70.1+:2.8
93.6122.0
86.5
85.8
84.6
95.0+1..5
70.4 (1)
13
86.3=2.7
8..7=1.3
83.9=2.0
Co-TSWD-1
88.2+11.4
70.41:1.2
94.7+00.9
16.1
80.2=1144
73.24:1.0
Co-ISWD-1-Wiin
88.0+
70.1+:1.2
9.66++1.0
85.7=0.7
80.1=112
72.9
Co-TSWD-5
88.7++1.7
71.0+1.55
96.7++0.8
91.5+004
82.0=0.9
75.4+00.7
Co-ISWD-5- Wavelet
88.5+1144
70.94:1.
96.5+00.7
91.1==0.6
82.3+:10 (t)
75.640.9 (t)
Co-TSWD-II
89.2+1111
71,4+:188
95.5++0.2
91.8+0.7
83.8=0.5
77.2+0.9
Co-TSWD-JO DESSWD-JO-T J0-Wavelet
89.0=10.9
71.21:4
95.4:10.5
91.3=0.8
83.7=0.7
76.92110
Co-SWCTWD
93.5+2..4
70.5+11.0
94.4+11.3
90.7=11.5
82.7=117
79.0=0099
Co-SWCTWD- -wwwee
93.21:1.6
71.3:0.9 (1)
94.2+:
90.3+11.7
82.3=2.0
79.020.5
Co-SWOTWD
96.2:11.2
72.4+2.1
96.0
90.6112.3
82.41:14
81331.1
Co-SWQTWD W-veeeee
95.9+10.9
71.8+:1.6
90.2+:2.0
82.0
80.9
Co-MST-TWD
88.71:2.4
68.413.3
91.312299
87.1
80.1=22.8
76.5+:
S..MSS:TWD-Waaelee
88.6+:1.3
68.74.2.
90.8++2.7
86.9+1.3
80.2+22.6 (1)
76.1+22.1
Co-TR-TWD
89.5+-1.2
70.9++1.7
93.4+-2.2
89.5+1.4
80.7=0.8
78.5++0.9
Co-TR-TWD-Waveler
90.6+0.7 (1)
70.9111.3
93,2++1.9
1.3
80.5=0099
78.510.7
Co-HHC-TWD
86.1+:2.1
70.1+:1.3
93.6+:1.5
88.5=0.5
83.2=1144
77.6
Co-HHC-TWD-Wayd
85.911.7
70.41:1. (7)
93.41:..1
89.2-0.9 (1)
82.7=100
77.5+0.9
Co-gHHC-TWD
84.0:12.0
70.41:1.
90.7111.7
87.211.9
1..9
84.2:11.22
Co-HHHHC-TWD-Wawck.
83.5
70.2+1.9
90.6:12.1
87.0=1.4
82.7=2.5
78.6211.6
Co-UltaFit-TWD
86.8100.
70.911.11
91,91110
89.9=2.0
83.7=2.9
79.1
Ww..ltaaHi-W---We... Wavelet
86.74-1.2
70.3+11.4
92. 12.1+00.00 (1)
88.7+1..8
83.4=2.7
79.0 1.7
QUE
84.7++0.5
72.4+-0.6
91.9++055
91.6=0.99
83.6#=144
82.5+t 1.9
Qaveleeelet
86.4++0.7 (f)
72.04-0.8
92.3++0.8 (t)
90.8+1100
85.1++1.0 (t)
82.4+ 1.6
WSV
85,9++1.0
71.4+:1.
92.6++0.7
89.0=1.5
81.6#2.4
77.5++1.7
Free-WSV
86.3::1.5
71.2+:19
92.41:10
88.7219
82.0=2.9
76.4+244
Tree- WSV-Wavelet
86.2+:1.1
71.5+ (1)
92.0+-088
1. (t)
82.0=2.3
78.1 (t)
Alg.
96.74-0.3
74.1-:0.5
97.3+0022
94.0=0.4
90.1=044
86.720.5
Alg. 2
97.3110.5 (1)
76.710077
97.610.1
94.2202
94.010.6
93.3
Tab. presents the performance integrating the wavelet filtering variants (*-MST", "-HC", and " Wavelet") into HGCNs for link prediction and node classification, compared to our full method that has iterative learning with wavelet filtering. The results show that the MST- and HC-based variants do not improve performance, while the improvement from the Wavelet-based variant is marginal. Applying wavelet filtering alone, without the iterative learning scheme, fails achieve performance comparable to our full approach. In contrast, our full method incorporates the relation between samples and features through TWD, and uses this structure inform wavelet filtering process. It plays key role in improving hierarchical representation and downstream task performance. To the best our knowledge, this is the first work to integrate wavelet filtering with TWD and apply within HGCNs.
51

---

Table 9: Performance comparison of HGCNs integrated with wavelet filtering variants ("-MST", "-HC", and "-Wavelet") versus our full method that combines iterative learning with wavelet filtering. Results reported for link prediction and node classification tasks across hierarchical graph datasets. ROC AUC for LP, Fl score for the DISEASE dataset, and accuracy for the AIRPORT, PUBMED, and CORA datasets for NC tasks. Arrows (1) denote cases with improvement. The highest performance is marked in bold, and the second-highest is underlined.
DISEASE
AIRPORT
PUBMED
NC
NC
NC
NC
HGONs [21]
90.8+00.3
74.5=099
96.4=0.1
90.6+:022
96.340.00
80.340.3
92.910.1
79.9+0.2
HOUNG-MST
92.3=0.2
88.4#0.4
91.7++02
779+00
90.4+0.3
77.6+066
HOCNs-HC
90.9-0.2
87.9+0.3
92.6+00.3
76.6+0.4
91.3+0.5
76.4+0..7
HGCNs
91.2=06660
76.7=110
96.0=002
90.62:0.3
96,410.1 (t)
80.14004
93.0=0.2111
80.420.311 (1)
12H UCH--GN
97.0.0.33
1..6.177
89.340.55
79,910.5
82.810.4
H2H-GCN-MST
90.5+00.3
88.6+0.4
91.4+0077
78.6++0.6
92.4403
81.3+0.5
H2H-GCN-HC
92.1=002
88.9204
92.6+0.4
79,8+02
92.5+0.5
81.7#00
H2H GCN- -CN-Waeeee
97.2-0.57 (*)
$7.918
90.1104 (1)
82.410.3
HOCN-AIg.
879-00. (1)
93.7+0.2 (t)
89.9+044
94.1++0.7
81.7+0.2 (f)
93.1+0.11 (t)
3
HGCN-AIg. 2
())4004)
89.4
97.2+0.1 (1)
92.140.3
97.2402
83,620.4
96.920.3
83.920
Table 10: Comparison of classification accuracy across word-document and scRNA-seq datasets for three variants: Alg. 1, Alg. W-Wavelet (post-hoc wavelet filtering after convergence), and Alg. (wavelet filtering applied iteratively during learning). Arrows (t) indicate improvements. We see that while post-hoc filtering offers minor improvements, iterative integration yields consistently better performance.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Alg.
96.7++0.3
74.1+0.5
97.3+0.2
94.0+004
90.1+004
86.7-0.5
Alg. 1-Wavelet
96.8+0.3 (f)
74.1+0.4
97.3+0.3
94.0+0.3
91.1+0.4 (t)
88.2+0.66 (f)
Alg. 2
97.3+0.5
76.7+0.77
97.6+0.1
94,2+002
94.0=20.6
93.3=0.7
F.6
Ablation study: effeet of post-hoc wavelet filtering vs. iterative integration
The main difference between Alg. 1 and Alg. 2 lies in the inclusion of filtering step each iteration of Alg. 2, designed to suppress noise and other undesired components of the data. This raises the question of whether applying wavelet filtering once Alg. has converged could also be beneficial.
To assess whether applying wavelet filtering as a post-processing step after the convergence of Alg. beneficial, we evaluate variant denoted Alg. 1-Wavelet. Tab. 10 demonstrates that this post-hoc filtering yields marginal improvements over the unfiltered version (Alg. 1). However, the performance remains consistently lower than that of the fully integrated approach in Alg. where wavelet filtering is applied iteratively during the learning process. This indicates that while post-pu post-processing offers slight gains, applying the filtering step within the iterative scheme more effective for hicrarchical representation learning.
F.7
Sparse approximation analysis via haar coefficient across iterations
Fig. 8 shows the evolution of the L1 norm of the Haar coefficients over iterations for all methods. Our methods exhibit consistent and monotonic decrease in the L1 norm as the iterations progress, reaching convergence. It reflects the sparse representations representation learning of both samples and features. Notably, sparsity is not explicitly enforced by our algorithm but emerges naturally from the learning process. In contrast, the competing methods are less stable over iterations. Some show slight reductions the L, norm, but the decrease not consistent and, in many cases, convergence is not reached. Even among those that converge, the final 7., norm values are higher than those attained by our methods, indicating less effective sparse approximation.
a result of the joint hierarchical
52

---

CBMC Sample Tree
Zelsel Sample Thee
1800
Alg.1 (ours)
1600
Alg.2 (ours)
1.0
160
co- Do-Quadtree
1200
co-Flowtree
1009
129
co-TSWD-1
100
co-TSWD-5
Co-TSWD-10
600
co TR TW
0
of iteration
10
iteration Zeisel Feature Tree
10
15
Co-HHC-TW
CBMC Feature Tree
QUE
1.00
Co-WCTWD
Co-WQTWD
co-SWCTWD
Co-SWQTWD
1.7
co-Ultra Tree
150
125
MSTTW
50
100
Co-gHHC-TW
0
10
10
co-UltraFit-TW
iteration
iteration CLASSIC Sample Tree
BBCSPORT Sample Tree
TWITTER Sample Tree
AMAZON Sample Tree
100
60
30
||0
||0.11
80
70-
20
60
30
40
50
#
iteration
BBCSPORT Feature Tree
TWITTER Feature Tree
CLASSIC Feature Tree
AMAZON Feature Tree
40
40
120
120
10,.11
||a./5
25
80
80
15
60
iteration
iteration
iteration
Figure 8: Evolution of the L1 norm of the Haar coefficients over iterations for all methods. Our methods exhibit a consistent decrease in the L, norm, indicating improved sparse approximation through joint hierarchical learning. Competing methods show less stable behavior, with many failing to reach convergence within 25 iterations.
F.8 Classification accuracy over iterations
We report the document and single-cell classification accuracy at each iteration for our method and the competing baselines. Fig. 9 presents the results both scRNA-seg and word-document datasets. Our method shows consistent improvement in classification accuracy over iterations, outperforning the initial performance after the first half iteration of Alg. 1, which corresponds using TWD without iterative refinement (see Tab. 6). This demonstrates the effectiveness our iterative learning scheme both data domains.
Moreover, our method consistently outperforms all competing baselines across iterations. While some baseline methods also benefit from iterative refinement, showing higher accuracy than their original TWD variants reported in Tab. 6, these gains are not uniformly observed. In addition, methods such as Co-Quadtree, Co-TSWD, Co-TR-TWD. and Co-HHC-TW: do not exhibit consistent improvements across all tasks, likely due to random sampling and initialization their tree construction procedures.
Importantly, we observe that our method converges within few iterations. While QUE and WSV also reach convergence, their final accuracy remains lower than ours. This highlights a key distinction: method explicitly learns hierarchical representation learning by jointly refining the hicrarchical
53

---

Alg. 1
60
6D Quadtree
co-Floutree
co TSWD
co TSWD 10
Co-TR.TWD
co HHCTWD
Reration
Reradoo
00CSPORT
TWITTER
CLASSIC
AMAZON
co WCTWD
Clacs &ccurayy 191
co SWCTWD
92
SWQTWD
90
co UltraTres
TWD
22
co 9-HHCTWD
UltraFi TWD
20
WSV
20
25
herstion
Peratton
Figure 9: Classification accuracy over iterations for scRNA-see seq and word-document datasets. Our method consistently improves accuracy across iterations and converges rapidly, outperforming both its non-iterative variant and all competing baselines. While some baselines benefit from iterative refinement, they show less consistent improvement and lower final accuracy.
structures of samples and features. It is further enhanced by the wavelet filtering, leading improved classification performance.
F.
Runtime analysis
10a
Ag
10a
co-TMM
co-T5MD-30
Run Time
co-TKTT
102
50 GND
CD-HHCTT
ce WMO
10a
ro-SWOTW
Co-MSTIW
grH- TV
Figure 10: Runtime comparison of our methods and competing baselines.
Fig. 10 presents the comparison of runtime between our methods and the competing baselines. Co- WMD and Co-GMD correspond to the WSV variants initialized with WMD and GMD, respectively. While our methods show slightly higher runtime than Co-TSWD-L and Co-TR-TWE TWD, they achieve significantly better sparse approximation and classification accuracy across all baselines (see Tab. I and Tab. 2). This trade-off demonstrates that our approach remains computationally efficient while providing superior performance.
F.10 Effect of Sinkhorn regularization in the iterative learning scheme
Our method adopts snowflake-baat based regularization function [57], which promotes smoothness across both rows and columns of the data matrix. T his regularizer increases monotonically over [0, oo) and penalizes small differences more strongly, encouraging local consistency. While effective, other regularization strategies may also considered. One notable alternative is entropic regularization [119]. which solves the OT problem using matrix scaling.
To explore this, we implement entropic regularization using the Sinkhorn solver [155] within our iterative learning framework. We denote these variants as Alg. 1-Sinkhorn and Alg. 2-Sinkhorn. Tab. reports the classification results on scRNA-seq and word-document datasets using the Sinkhorn regularization. In all cases, incorporating the Sinkhorn penalty leads noticcably lower
54

---

Table 11: Classification aceuracy scRNA-seq and word-document datasets using Sinkhorn- regularized variants of our methods. "-Sinkhomn" denotes the use of entropic regularization in place of the snowflake penalty. Across all datasets, the Sinkhorm-based variants perform worse than the original methods, highlighting the effectiveness of the snow flake regularization.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Alg. Alg. I-Sinkhorn
96.7+L0.3
74,1+0.5
97.3-0.22
94.0=0.4
90.14:0.4
86.7.10.5
90,310.5
70,910.4
94.010.5
91.2=0.5
86.3+0.5
82,4=0.6
Alg.
97.3:10.5
76.740.77
97.620.1
94.2=0.2
94.0110.6
93.320.7
Alg. 2-Sinkhorn
94.2+0.4
73.8+0.5
95,4+0.3
91.5-0.1
89.7+0.55
85.9+0.6
Table 12: Clustering performance on gene expression datasets (Breast Cancer and Leukemia) evalu- ated using Clustering Accuracy (CA), Normalized Mutual Information (NMI, and Adjusted Rand Index (ARI). Our methods (Alg. I and Alg. 2) outperform all baselines. Alg. 2. which incorporates wavelet filtering, achieves the highest overall performance. Results for BCOT, CCOT, CCOT-GW. and COOT are reported from the work of Fettal et al. [156]. OOM indicates methods that ran out of memory.
Dataset
Breast Cancer
Leukemia
Evaluation Metric
CA
NMI
ARI
CA
NMI
ARI
QUE
71.911.2
45.3
39.810.4
72.3
56.6106
42.910.4
WSV
72.4
33.8122
39.0108
64.31441
61.7123
44.8
Tree-WSV
73.1
36.7+1.7
41.5+09
66.2+3.0
60.4+1.8
42.643.1
BCOT
76,9
37.2++0.
26.7+000
71.2+54
59.6+6.9
39.9+6.
BCOTX
84.6=0.0
48.3+000
46.0+0.0
80.9+3.8
70.9
55,3+3.
CCOT
0OM
OOM
OOM
40.6400
0.0+n.0
0.0
CCOT-GW"
OOM
OOM
OOM
0OM
OOM
OOM
COOT
63.1
5.4
-1,2129
36.2427
14.013.6
5.4
COOT
61.5100
5.4
2.2
32.5133
8.7.227
-0.5.2.1
Alg. 1
85.6+0.7
49 19.2+0.5
50.4+0.9
75.8+0.9
64.2+110
54.2+1.
Alg. 2
87.2+0.4
51.3+0..
57.6+0.7
81.2+1.1
75.4+08
59.7+12
performance compared to the original versions of our algorithms. The classification accuracy degrades, indicating that the entropic regularization may overly smooth or distort the learned hierarchical structures. This suggests that the snowflake regularization is more suitable in our setting, where the goal is to represent fine-grained hierarchical relationships across both samples and features.
F.11 Co-Clustering performance on gene expression data
We further evaluate the effectiveness of our methods in clustering tasks using two gene expression datasets [157]. The first breast cancer (BC) dataset containing 42,945 gene expression values across 26 samples divided into two classes. The second is leukemia (LEU) dataset with 22,283 genes measured across patients divided into five classes. We apply Alg. and Alg. to the data and perform k-means clustering [158] on the resulting sample TWD matrices to obtain clustering assignments. We evaluate the OT-based co-clustering baselines including: QUE [52], WSV [59], Tree- WSV [60]. biclustering with optimal transport (BCOT and BCOT,) [156], Co-clustering through OT (CCOT and CCOT-GW) [159], and Co-optimal transport (COOT and COOT,) [160]. We follow the evaluation setup [156], reporting Clustering Accuracy (CA), Normalized Mutual Information (NMI) [161], and Adjusted Rand Index (ARI) [162]. The results of BCOT, CCOT, CCOT-GW, and COOT reported from the work of Fettal et al. [156].
Tab. shows that our methods outperform all competing baselines. Alg. already achieves strong performance, while Alg. 2 further improves upon it with consistently higher CA, NMI, and ARI scores. particular, the improvements are pronounced on the Leukemia dataset, where our method significantly exceeds the basclines. Some methods, such CCOT and CCOT-GW, run out of
55

---

Algorithm 4 Incorporating Fixed Hierarchical Distance for One Mode in Iterative Learning
Input: Data matrix X
Mg 1T Output: Trees T(W," (1) and T(W?). and TWDs w and W
R mmxm and M, R=x", "i and Yex thresholds de and 0r,
t O, x
X. Z (0) XT, W M.
Initialization
repeat
X(I-1) v(X (0; T(W2))
Tree haar wavelet filtering with updated tree Tree haar wavelet filtering with fixed tree
Z VI T(M,))
x!-1)
(((t)) (x
X!-1)
/ (z
w
t 2 (Xx (x(++));
(TW!!)
Iterative update with updated tree Iterative update with fixed tree
W!
$ (X (1+1) -T(M))
1+1+1
until convergence
Table 13:
Link prediction and node classification performance on hierarchical graph datasets using different variants of the proposed HGCN framework. HGCN-Alg. uses only the reference update without filtering. HGCN-Alg. updates both filtering and referencing functions adaptively. HGCN- Alg. 4 applies filtering with fixed sample distance matrix M, dи- Results show that while fixing M, improves over the reference-only variant, superior performance is achieved when both filtering and referencing steps are adaptively updated throughout the iterations.
DISEASE
AIRPORT
PUBMED
CORA
LP
NC
LP
NC
LP
NC
LP
NC
HGCN-Alg.
93.211..6
87.910.7
93.7=0.2
89.910.4
94.110.7
81.7=0.2
93.110.1
82.910.3
HGCN-AIg. 4
95.1+0.55
88.2+0.9
95.9-0.1
90.7+0.3
95.710.4
82.0-0.3
95.3-0.1
83.2+10.4
HGCN-AIg, 2
98.4:00.4
89.4+0.33
97.2=0.1
92.11103
97.240.2
83.6=0.4
83.9:10.2 96.9-0.3
memory (OOM), and others like COOT and COOT produce poor clustering results. These findings demonstrate the robustness and effecttivneess of our method, especially when enhanced by wavelet filtering at each iteration.
F.12 Incorporating fixed hierarchical distances in the iterative learning scheme for HGCNs
In Sec. 5, we demonstrated how prior knowledge of a hierarchical structure be incorporated by initializing the sample distance matrix M, with the shortest- shortest-path distances dæ induced by given hierarchical graph H. This initialization graph data using the HGCN framework. Beyond initialization, another way to incorporate dæ into our iterative learning scheme is to fix the sample pairwise distance matrix M, = du during the filtering step (-) and the joint reference computation (-.).
of the modes can
evaluated in the context of hierarchical
However, it is important to distinguish between the two variants of our learning scheme. The iterative algorithm in Alg. only involves the reference computation (--). Therefore, if M, is fixed within (-)), no update occurs after the first iteration, and the procedure halts. In contrast, Alg. 2 includes an additional filtering step (-) that continues to evolve even when the used distance in ((-) remains fixed. As result, fixing M, is only meaningful within Alg. 2, where the iterative updates can still proceed. We summarize this approach in Alg. 4 as an alternative variant that allows the incorporation of external hierarchical information. This approach is also evaluated on hierarchical graph data using the HGCNs framework.
Tab. 13 reports the link prediction and node classification results when the sample distance matrix M, is fixed the iterative learning scheme (denoted as HGCN-Alg. 4). Compared to HGCN-Alg. 1, fixing M, and applying filtering improves performance. This shows the benefit of incorporating the
56

---

Table 14: Link prediction and node classification results for GCN, GRAPHSAGE, and HGCNs, with and without initialization using our learned hierarchical representations ("-Alg. 2").
DISEASE
PURMED
LP
NC
L.P
NC
I.P
NC
LP
NC
GCN
64.710.5
69.7104
89.3.1044
8141066
91.120.5
78.120.2
90.4.0.2
813100.3
OCN-Alg.
68.110.4 (*)
(1,2=0.(())
90.610.2
(3,0109(1)
90.520.3
78.710.1 (h)
91.2104
81910.2
GRAPHSAGE
63.9100.3
69.12006
90.110.55
$2.110.5
86.21110
7741222
85.5.006
77.912.1
GRAPHSAGE-ALg.
67.8
70.5=00.3
90.2:0.4
83.740.5 (f)
90.3= (f)
78.81177
89.640.8 (f)
$0,5#1.6(h)
HGCNs
90.840.3
74.5409
964420.1
90.40022
96.3=00
80.3-0.33
929--011
79.940.2
97.010.3
88.61117
89.3.0.5
96:0=:0.0
82.8.0.44
HGCN-AIg
(*)4+04())
89.4-0.3
97.2+0.1 (1)
92,1+003((1)
91.2-0.2 (1)
836+0441) (1)
96.90.33 (t)
83.9+02 (f)
Table 15: Comparison of representation quality for Independent MST, Co-MST-TWD, and our full method (Alg. 2), measured by the L1 norm of Haar coefficients (lower is better). Results are reported as samples features for each dataset. Lower values indicate that the learned hierarchical representations more effectively capture the hierarchical information across the two data modes.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Independent MST
37.9/42) 42.8
82.4/311
83.6/ 131.4
109.6 / 144.1
194.3 /
1928.7 / 197.6
Co-MST-TWD Alg.
34.1
69.1/22.1
77.5/ 109.3
97.0/ 126.4
179.1 /25400
1755.3 / 83.7
10.1/4.8
22.4/3.4
37./ 19.4
46.6 / 26.6
50.9 /9.9
489.4 / 45.6
hierarchical structure during filtering. However, HGCN-Alg. 2, which adaptively updates both filtering and referencing functions via the evolving ( along the iteration, consistently outperforms the fxed-distance variant. This highlights the importance of mutual adaptation between the sample and feature structures during the iterative learning process, where the joint updates allow the hierarchical relations be refined jointly across iterations.
F.13 Evaluating the preprocessing step with additional backbones
To evaluate the effectiveness of our method preprocessing step across different neural network backbones, we extend our experiments to include GCN and GraphS in addition to the hyperbolic models considered in Tab. 3. Tab. 14 reports link prediction and node classification results, where models initialized with our Icarned hierarchical representations (denoted " "-Alg. 2") consistently outperform their baselines. Improvements are indicated with an upward arrow (t). These findings demonstrate the broad applicability of our approach across diverse neural architectures. Among all models, HGCN-Alg. 2 achieves the highest overall performance, which we attribute to the advantages of hyperbolic representations in modeling hierarchical structures, distinguishing HGCNs from their Euclidean counterparts.
E114 Discussion and comparison with other tree distance metrics
In our framework, TWD plays a central role due to its ability compare probability distributions supported on hierarchical structures. This is essential for our alternating, joint-learning process, where the hierarchical structure from one mode (e.g., samples) used to inform and refine the hierarchy in the other mode (e.g., features), and vice versa. That is, the tree structure provides the foundation for expressing hicrarchy, while the Wasserstein metric provides the mechanism for distributional comparison and transfer from mode to the other. This cross-mode interaction is central to method and relies on the Wasserstein nature of TWD.
TWD differs from other tree distance metrics such as minimum spanning tree (MST) [163], neighbor joining, or low-stretch trees, which typically operate on tree structures or discrete node sets but do not directly support distribution- based comparison needed for cross-mode refinement, To clarify this difference, we consider two alternative setups using MST: (1) using independently constructed MSTs for samples and features without cross-mode interaction, and (2) using an MST as the tree structure in the TWD framework for our iterative learning process (Co-MST-TWD), We compare these alternatives to full method (Alg. 2) by measuring representation quality using the L
57

---

Initial
Iter
User
Video
User
Vidco
Iter 5
Iter
User
Video
User
Video
Cookap ucn
kep-irrir
brccing banly Caval-icnen Corrraceoos
Ialry Bigpsty
hussy
Camoty
Drara
Acki
Figure 11: Stepwise hierarchical refinement on the toy user-video recommendation in App. F.I.
norm of Haar coefficients. Results are shown Tab. 15 (format: samples features). see that both alternatives yield substantially higher L1 norms. These results suggest that the learned joint hierarchical representation from our iterative process better reflects the hierarchical information of the two data modes. We attribute this advantage to the specific design of proposed method, where the TWD supports iterative cross-mode refinement, which is not feasible with independent tree metries.
Replacing TWD with other graph-based distances is conceptually possible, but there are important considerations:
Most graph-basee based distances are designed to compare node positions or global connectivity in general graphs, measure how distributions supported on tree reveal explicit nested groupings.
Gromov- Wasserstein distance [164] could, in principle, be integrated into an iterative scheme compares distributions across different metric spaces graphs. However, such approaches are much more computationally intensive.
Another possible direction is construct graph from high- sigh-dimensional data and transform this general graph into unrolling trees [165]. Further investigation would be needed to determine if and how such transformations can be integrated within iterative learning process.
E.15 Stepwise hierarchical visualizations in recommendation toy problem
Fig. 11 presents the stepwise visualizations based the toy video recommendation example in App. F.1 for hierarchical refinement. The iterative refinement process starts with an initial tree structure, which is unstructured since the preliminary grouping based on initial data relationships (M, or M,). At this point, major clusters are not well formed, and samples or features from different categories are mixed. After the first iteration, the algorithm uses the current (e.g., users) to compute TWD for the other mode (e.g.. videos). This leads to newly induced tree that starts to separate the major categories, and major branches start to reflect groupings, but subcategory assignments are not fully sorted. As the algorithm proceeds through subsequent iterations, alternately refines the trees for users and videos. The updated tree for one mode informs a more accurate computation of TWD in the other mode, enabling finer divisions within the existing branches. Subcategories like action, drama, and sci-fi (for fiction videos) become more distinct and better grouped. After several iterations, the subcategories become more cohesive and consistently organized, and the tree structures stabilize. It illustrates how the trees evolve with each iteration and points out key groupings formed in the learned trees.
from one mode
G Additional remarks
provide additional remarks on our methods, including the motivation for using wavelet filtering, and the comparison with WSV and Tree-WSV.
58

---

G.1 Motivation for tree-based wavelet filtering
To refine the hierarchical representations learned through our iterative learning scheme, we incorporate filtering step (-) at each iteration. This filtering is based on Haar wavelets [38-41] induced by the learned trees. Each sample X E 1n can be viewed as a signal defined over the feature tree T(Mc), whose leaves correspond to features [42]. Given feature tree T T(M,), we construct an orthonormal Haar basis B. E Rm xm as described in Sec. 3 and App. A. Each sample is expanded in this basis, yielding coefficients O; = (X,,, B.)) E R where (j) (X 3y) and (3, is the j-th wavelet in B,. To define a wavelet filter, we select subset of the Haar basis vectors as follows. For each coefficient index j. compute the aggregate I. norm Σi=1 [xx( across samples and sort the indices in descending order. We then select the subset Ω such that the cumulative contribution 19 = Eee اi=1 (q) exceeds a threshold 9. > 0. similar filtering procedure can be applied to features using the Haar wavelet construct by the sample tree.
is important to note that this L -bssd filtering retains the components that contribute most to the hierarchical decomposition, as represented by the tree [69]. We postulate that the residual, i.e., the remaining coefficients outside the set 22 the variation that is not aligned with the tree structure. These components are treated as noise or nuisance terms that are irrelevant for the hierarchical representation learning. By discarding this residual, the filter emphasizes signal content that is structurally consistent with the hierarchical structure, leading to representations that are more stable and meaningful across iterations [166, 167]. Our filtering step adapts to the learned trees at each iteration. That is, the Haar wavelets are constructed from the trees T T()()) and 'T (W()) via TWD, making the filtering aligned with the evolving hierarchical representations. Empirically, we show that applying the wavelet filtering in the iterative learning scheme improves performance in various tasks on hierarchical data, including link prediction and node classification on hierarchical graph data, as well as sparse approximation and metric learning on scRNA-seg and document datasets.
G.2 Advantages of Haar wavelet filtering over Laplacian-basee filtering for hierarchical data
In Sec. 4.2, we introduced a filtering step based on the Haar wavelet, where the wavelet basis is induced by the tree structure [43] learned during our iterative scheme. This filtering step is designed reflect the intrinsic hicrarchical structure of the data, enhancing meaningful hicrarchical representations while suppressing noise and other nuisance components. An alternative filtering approach could involve using the Laplacian matrix derived from the tree graph. This would allow spectral filtering based on the eigendecomposition of the tree [42, 168], offering another way to apply filtering that incorporates the hierarchical information from the tree [169].
Consider a feature tree T(M.) with the affinity matrix A. RRmx. We can construct the Lapla- cian matrix for the feature tree by L, = D. A,: where D. is a diagonal degree matrix with D.(j.j) = A(i The eigenvectors of the Laplacian matrix are then treated as an orthonormal basis. often referred to as the graph Fourier basis [170-172]. Subsequently, each sample can be expanded in this eigenbasis, and we denote ā; (X,,. U)) E R as the vector of the expansion coefficients of the i-th sample, where U, is the matrix consisting of eigenvectors of L.. This expan- sion allows the interpretation of the Laplacian expansion coefficient āx as the contribution of the j-th frequency component to the signal. By sorting these coefficients according to the corresponding Laplacian eigenvalues, we can construct low-pass filter (containing the first few expansions) that suppresses
components and noise while promoting signal smoothness [173]. A similar decomposition can also be applied to the features using the Laplacian matrix of the sample tree.
We evaluate this alternative filtering technique based on the Laplacian matrix of the tree graph. Specifically, we modify the iterative scheme in Alg. 2 by replacing the Haar wavelet filtering step with the low-pass Laplacian-bart based filtering, and refer to this variant as "Alg. 2-Laplacian." " We test this variant on both document and single-cell classification tasks, as well as on hierarchical graph data for link prediction and node classification. Tab. 16 reports the classification accuracy for the document and single-cell tasks. We see that compared to the wavelet- based approach, the Laplacian variant consistently underperforms, showing a noticeable drop. This suggests that while Laplacian filters offer spectral interpretation, they are less effective at representing hierarchical information in our setting. Similarly, Tab. 17 presents link prediction and node classification performance for integrating the Laplacian filtering into our iterative learning scheme and applied to HGCNs. Here too, we observe that AGgN-Alg, 2-Laplacian performs worse than HGCN-Alg. 2. These results show the
59

---

Table 16: Classification accuracy document and single-cell datasets. We compare the performance of the iterative method without filtering step (Alg. 1), our full method with Haar wavelet filtering (Alg. 2), and the Laplacian filtering variant (Alg. 2-Laplacian). Haar wavelet filtering consistently achieves the better performance.
BBCSPORT
TWITTER
CLASSIC
AMAZON
ZEISEL
CBMC
Alg.
96.7-0.33
74.140.5
97.3=1022
94.0::044
90.1
86.7+0..5
Alg. 2-Lappacian 95.2-0.7
72.6+0.9
94,4+1.0
91.2++0.4
83.7+009
81.6+11.2
Alg.
97.3-0.5
76.7+0.77
97.6+0.1
94.2++0.2
94.0++0.6
93.3+0.77
Table 17: Performance of different filtering variants in our iterative learning scheme applied to hierarchical graph data within the HGCN framework for link prediction and node classification. We compare HGCNs using the hierarchical representations from Alg. I, the variant with Laplacian filtering (HGCN-Alg. 2-Laplacian), and the full method with Haar wavelet filtering (HGCN-Alg. 2). The Haar wavelet variant consistently achieves superior performance.
DISEASE
AIRPORT
PUBMED
CORA
LP
NC
LP
NC
LP
NC
LP
NC
HGCN-AIg.
93.2++0.6
87.9+0.7
93.7+0.2
89.9+0.4
94.1+00.7
81.7+0.2
93.1+0.1
82,9+0.3
HGCN-AIg. 2-Laplacian
89,4++0.3
79.2+1111
93.1
88.4100.7
93,540.4
80.9+0.2
92.7=0.2
81.940.3
HGCN-AIg.
98.4++0.4
894++..3
97.2+0.1
92.1+0.3
97.2+0.2
83.6=0.4
96.9=0.3
83.9+0.2
benefit Haar wavelet filtering, which more directly reflects the hierarchical structure and is more effective for learning meaningful representations.
We attribute this performance gap to the differences between the two filtering methods. Laplacian- based filtering is not designed to represent hierarchical structure explicitly [42, 72]. Its eigenvectors global and do not distinguish between different levels of abstraction in the data. In addition, from a computational perspective, computing the eigendecomposition of the graph Laplacian can be expensive: full decomposition scales cubically with the number of nodes, while partial decomposition (e.g.. top k eigenvectors) still incurs significant cost, often quadratic or worse in (depending on the method and graph sparsity). Furthermore, it has been shown that the expansion coefficients in the Haar wavelet basis decay at an exponential rate (Prop. A. in contrast to the typically polynomial rate of decay observed with Laplacian cigenfunction coefficients [43]. This rapid decay reflects the ability of Haar wavelets to provide compact representations of piecewise-smooth signals, using the learned tree structure to concentrate energy in fewer coefficients. Finally, Haar wavelets on trees yield interpretable components [38]: each wavelet coefficient corresponds to specific scale and position within the tree, often representing differences between sibling nodes or between a node and its parent. This locality and hierarchical structure are absent in Laplacian eigenbases. Therefore, for data where a hierarchical relationship between samples or features important, Haar wavelets provide a more structurally aligned approach compared graph signal filtering based on the Laplacian.
We note that more advanced filtering strategies, such as those learned through neural networks [96]. could also be explored. We leave this extension for future work.
G.3 Comparison with Tree-WSV and the role of hierarchical representation learning
While preparing this manuscript, a new method called Tree-WSV [60] was introduced. It uses Tree- Wasserstein Distance (TWD) to place samples and features leaves trees and was introduced to reduce the complexity that follows ideas from existing TWD research [35, 34j. Specifically, this approach speeds up WSV [59] for unsupervised ground metric learning.
The main idea of WSV is to simultaneously compute the Wasserstein distance among samples (resp. features) by using the Wasserstein distance among features (resp. samples) as the ground metric. This design conceptually aligns with co-manifold learning methods, which use relationships among samples to guide the relationships among features [52-55, 7, 56, 57, 63]. Tree-WSV extends WSV by using TWD as a low-rank approximation of the Wasserstein distance. does this by computing TWD among samples (resp. features) while using TWD among features (resp. samples) as the
60

---

ground metric. Tree-WSV reduces the complexity of WSV from O(n2m² n + m logm). m)) to O(n"+ +m+++m+-- +m²++m) mm)), making more efficient practice. In this context, the constructed trees primarily serve as computational tools for approximating Wasserstein distances.
In contrast, work takes geometric approach where the trees induced by TWD used to represent hierarchical structures among samples and features. Based on this geometric perspective, we propose Alg. and Alg. 2, which jointly learn hierarchical representations for both samples and features. Rather than using TWD merely as a computationally efficient approximation, we treat it as an input distance to construct diffusion operators that encode manifold information [68]. These operators are then used to construct trees based on hyperbolic embeddings and diffusion densities [37], which recursively guide the computation of TWD for the other data mode.
This geometric view relates more closely to co-manifold learning in spirit, which use the geometry of one data mode inform the geometry of the other, and was explored in the WSV Tree- WSV frameworks. Our primary goal is not unsupervised ground metric learning (though this arises naturally) but rather joint samples and features hierarchical representations. We also use multiple resolutions to better represent the hierarchical structure of the data, similar to ideas found in tree-based Earth mover's distance [52, 57J.
Empirically, our methods outperform WSV and Tree-WSV across tasks involving hierarchical data, including sparse approximation, document classification, and single-cell analysis, as detailed in Sec. 6. These improvements highlight the benefit of explicitly modeling and learning hierarchical representations rather than using trees solely for efficient computation. Furthermore, demonstrate that the learned hierarchical representations can serve as effective inputs for downstream applica- tions, such hyperbolic graph neural networks (HGCNs), where incorporating our wavelet-based hierarchical representations leads to improved link prediction and node classification performance. This integration further illustrates the utility and effectiveness of our approach beyond unsupervised ground metric learning.
G.4 On the interpretability of TWDs
Note that strength of TWD lies in its transparent and decomposable structure. Given two distributions (e.g.. feature or sample vectors), the TWD quantifies the minimal "mass transport" needed to align them over the tree, with transport cost reflecting branch length. This formulation enables users to understand where and at what scale differences occur, For instance, if most of the transport cost is concentrated at high-level branches, the two distributions differ in broad terms; if the cost is concentrated near the leaves, they are largely similar but differ in fine details. As a result, TWD not only quantifies similarity but also supports qualitative interpretation: users can trace which parts of the tree contribute most to the distance between any two data points. In classification tasks, this interpretability reveals which regions of the hierarchy drive class separation and can help identify key transitions in the data.
61


Reviews:
[{"review_id": "6QALKqbPfW", "reviewer": "Reviewer_pif9", "summary": "This paper proposes a novel representation learning method for data with hierarchical structure in both sample dimension and feature dimension. The method proceeds by refining Tree-Wasserstein Distance for both dimension in an alternative manner. The authors also propose to use Haar filter to reduce noises and even show that the proposed method can work well on graph data. Overall, this is a novel and interesting work with solid theoretical insights and extensive experiment results.", "questions": "Please see weaknesses", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n- The paper is targeting an important and common representation learning scenario with a novel TWD-based method; the alternating refinement step with proved convergence guarantee is reasonable\n- The experiments are extensive and solid - I especially like the one that shows proposed method can work well as an initialization step for hyperbolic graph data, which is an interesting application scenario.\n\n**Weaknesses:**\n - I think some terms in this paper could be presented with more clarity. For example, it was not clear to me what the \"two data modes\" of \"samples and features\" mean until I read through to the second paragraph. It could be mentioned earlier that they correspond to rows and columns of an m by n data matrix. \n-  The authors may consider further integrate the Haar wavelet filtering into Figure 1's illuatration", "was_revised": true, "final_justification": "This is a novel paper and my concerns have been resolved as acknowledged. I will maintain the recommended score.", "rebuttal": "We thank the reviewer for the thoughtful and encouraging review, and for recognising the relevance of addressing an important and common representation learning scenario, as well as its application to hyperbolic graph data. In terms of weaknesses: \n\n>**Clarifying the terminology regarding data modes**\n\nThank you for this valuable suggestion. We agree that introducing \"samples\" and \"features\" as corresponding to the rows and columns of an $n \\times m$ data matrix at the very beginning would improve the clarity and accessibility of our paper. In the revision, we will update the first paragraph of the introduction to clearly state early on that the \"two data modes\" refer to the rows (\"samples\") and columns (\"features\") of the data matrix. We hope this clarification will help readers immediately understand the problem setting. If there are any further suggestions for clarifying the terminology, we would be happy to incorporate them. Thank you! \n\n\n>**Including Haar wavelet filtering in the illustration in Figure 1**\n\nWe thank the reviewer for the helpful suggestion. We agree that including the Haar wavelet filtering step into Figure 1 would provide a more complete visual summary of the overall method. In the revised version, we will add a filtering block below each learned tree to represent the application of wavelet filtering to the data. This filtered data is then used in the computation of the TWDs for the subsequent iteration. We hope this addition will enhance how the filtering is induced by the learned tree, how it modifies the data between iterations, and how it contributes to refining the joint hierarchical representation. Thank you again for this helpful recommendation.\n\n---\n\nWe thank the reviewer again for their constructive comments and hope to discuss any ongoing suggestions further."}, {"review_id": "mme7QRFQhJ", "reviewer": "Reviewer_7yBG", "summary": "This paper focuses on learning hierachical representations via Tree Wasserstein Distance (TWD). The paper proposes a novel tree-based representation algorithm that iteratively aligns the feature and sample hierachy by optimizing the TWD distance with coordinate descent method. A variant based on Haar wavelets is introduced to reduce the noisy components in the constructed feature tree. Experiments on word-document and sRNA-seq datasets verify the quality and effectiveness of the proposed representation algorithms.", "questions": "1. Is it possible to replace TWD with other graph-based distance while preserving the hierachical features in the outputs?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n+ A novel framework that generates informative document and word-level hierachical representations.\n+ The convergence of the methods is grounded by theoretical analyses and demonstrations.\n+ Comprehensive experiments validate the effectiveness of the algorithms.\n\nWeaknesses:\n- As the core of the optimization algorithms, the tree construction approach and the TWD metric are both adopted from prior research works.\n- Due to the discerete nature of TWD, it appears challenging to integrate the proposed optimization objective into other neural network-based models. On the other hand, the effectiveness of the algorithm as a preprocessing step should be evaluated among a wider range of backbone models, such as GCN and GraphSAGE.\n- While only the tree-based Wasserstein distance is utilized as the similarity metric, its advantage over other tree distance metrics is not adequately discussed or explored via experiments.", "was_revised": true, "final_justification": "The rebuttal message addresses the majority of my previous concerns and clarrify the main contribution of the paper, therefore I would raise the score.", "rebuttal": "We thank the reviewer for their insightful and in-depth evaluation of our paper, and for recognising the theoretical grounding of our method. In terms of the weaknesses and questions: \n\n>**Use of established methods and our novel contributions**\n\nWe thank the reviewer for commenting on this point. We agree that our method indeed builds on existing components: the TWD as introduced in prior work by Indyk and Thaper (2003), and tree construction methods based on diffusion geometry and hyperbolic embedding as proposed by Lin et al. (2025). However, our main contribution lies in the formulation of the joint hierarchical representation learning problem for both samples and features, and in the development of a novel iterative algorithm that alternates between these two modes, so that each hierarchy helps improve the other and jointly refines both. The proposed algorithm is new, as are the theoretical convergence guarantees we establish and the utility of our methods as a preprocessing step for HGCNs. This alternating optimisation, along with the Haar wavelet filtering over the learned trees, constitutes a new approach that extends beyond prior uses of TWD or hierarchical representations. In addition, the iterative process is broadly applicable to various TWD-based methods and lays a foundation for future extensions to more general geometric representations beyond trees, while retaining the proposed iterative procedure based on the Wasserstein distance learned between samples and features. We will highlight more explicitly what is novel in our formulation, optimisation scheme, and applications in the revision. Thank you for bringing this up, and we will better clarify the novel contribution of this work in the revision.\n\n>**On the integration of TWD with neural architectures**\n\nWe thank the reviewer for this insightful feedback. These are valuable suggestions that could extend the scope of our work.\n\nIndeed, due to the discrete nature of TWD, direct integration into end-to-end neural network training is non-trivial. In this work, we focus on using our method as a preprocessing step that produces hierarchical representations that can then be used as input to initialise hyperbolic graph convolutional networks, where we observe performance improvements. While this is the current scope of our work, the reviewer's comment actually motivates a promising direction: exploring a differentiable variant of TWD, such as the soft TWD proposed by Takezawa et al. (2021), could enable integration of our iterative process into neural architectures. In Takezawa's work, the TWD is differentiable with respect to the probability of the parent–child relationships of a tree.  In our setting, these relationships are derived from the diffusion operator introduced by Coifman and Lafon (2006), for which gradients can be computed. While a full integration would require further investigation, this connection suggests that a differentiable extension of our iterative framework may be feasible and could further extend the applicability of our method within neural architectures. We will add this direction for future work in the revised manuscript. Thank you again for motivating us to consider this line of research.\n\n>**Evaluating the preprocessing step with additional backbones**\n\n In terms of evaluating the effectiveness of our method as a preprocessing step across different neural network backbones, we appreciate the suggestion to extend our experiments beyond the set of models originally included. In response, we expanded our experimental analysis to include GCN and GraphSAGE, in addition to the hyperbolic models originally considered. We report the link prediction and node classification performance in the table below. We find that initialising these models with our learned hierarchical representations  (denoted \"-Alg.2\") consistently improves performance in node classification and link prediction tasks. Improvements are marked with an upward arrow ($\\uparrow$). These results support the general applicability and effectiveness of our method as a preprocessing step across different neural architectures. Notably, we see that HGCN-Alg.2, which incorporates our method as an initialisation within HGCNs, still achieves the best overall results. We attribute this to the unique advantages of hyperbolic representations for modelling hierarchical structures, which distinguishes HGCNs from other Euclidean-based models. We will include these results in the revised manuscript. Thank you again for this valuable suggestion.\n \n | | DISEASE LP| DISEASE NC| AIRPORT LP| AIRPORT NC| PUBMED LP| PUBMED NC| CORA LP| CORA NC|\n|------------------|-------------------|-------------------|-------------------|-------------------|-------------------|-------------------|-------------------|-------------------|\n| GCN| 64.7±0.5| 69.7±0.4| 89.3±0.4| 81.4±0.6| 91.1±0.5| 78.1±0.2| 90.4±0.2| 81.3±0.3|\n| GCN-Alg.2| 68.1±0.4 (↑)| 71.2±0.5 (↑)| 90.6±0.2 (↑)| 83.0±0.9 (↑)| 90.5±0.3| 78.7±0.1 (↑)| 91.2±0.4 (↑)| 81.9±0.2(↑)|\n| GraphSAGE| 65.9±0.3| 69.1±0.6| 90.4±0.5| 82.1±0.5| 86.2±1.0| 77.4±2.2| 85.5±0.6| 77.9±2.4|\n| GraphSAGE-Alg.2  | 67.8±0.1 (↑)| 70.5±0.3 (↑)| 90.2±0.4| 83.7±0.5 (↑)| 90.3±1.1 (↑)| 78.8±1.7 (↑)| 89.6±0.8 (↑)| 80.5±1.6 (↑)|\n| HGCNs| 90.8±0.3| 74.5±0.9| 96.4±0.1| 90.6±0.2| 96.3±0.0| 80.3±0.3| 92.9±0.1| 79.9±0.2|\n| H2H-GCN| 97.0±0.3| 88.6±1.7| 96.4±0.1| 89.3±0.5| 96.9±0.0| 79.9±0.5| 95.0±0.0| 82.8±0.4|\n| HGCN-Alg.2| **98.4±0.4 (↑)**| **89.4±0.3 (↑)**| **97.2±0.1 (↑)**| **92.1±0.3 (↑)**| **97.2±0.2 (↑)**| **83.6±0.4 (↑)**| **96.9±0.3 (↑)**| **83.9±0.2 (↑)**|\n\n>**Discussion and comparison with other tree distance metrics, and replacing TWD with other graph-based distances while preserving hierarchical features**\n\nWe thank the reviewer for the insightful comment and question.  In our framework, TWD plays a central role due to its ability to compare probability distributions supported on hierarchical structures. This is essential for our alternating, joint-learning process, where the hierarchical structure from one mode (e.g., samples) is used to inform and refine the hierarchy in the other mode (e.g., features), and vice versa. That is, the tree structure provides the foundation for expressing hierarchy, while the Wasserstein metric provides the mechanism for distributional comparison and transfer from one mode to the other. This cross-mode interaction is central to our method and relies on the Wasserstein nature of TWD.\n\nWe agree that the paper would benefit from a more explicit discussion of how TWD differs from other tree distance metrics, e.g., minimum spanning trees, neighbour joining or low-stretch trees. These tree distance metrics typically operate on tree structures or discrete node sets, but do not directly support distribution-based comparison that allows for cross-mode refinement. To clarify this difference, we consider two alternative setups using MST: (1) using independently constructed MSTs for samples and features without cross-mode interaction, and (2) using an MST as the tree structure in the TWD framework for our iterative learning process (Co-MST-TWD). We compare these alternatives to our full method (Alg. 2) by measuring representation quality using the $L_1$ norm of Haar coefficients. Results are shown in the table below (format: sample / feature). We see that both alternatives yield substantially higher $L_1$ norms. These results suggest that the learned joint hierarchical representation from our iterative process better reflects the hierarchical information of the two data modes. We attribute this advantage to the specific design of our proposed method. In the revised manuscript, we will emphasise more how TWD supports iterative cross-mode refinement, which is not feasible with independent tree metrics, and report the empirical comparison to support our use of TWD. We thank the reviewer again for encouraging us to address this important point more thoroughly.\n\n||BBCSPORT|TWITTER|CLASSIC|AMAZON|ZEISEL|CBMC|\n|------------------|---------------|---------------|---------------|---------------|----------------|----------------|\n| Independent MST| 37.9 / 42.8| 82.4 / 31.1| 83.6 / 131.4| 109.6 / 144.1 | 194.3 / 281.5| 1928.7 / 97.6  |\n| Co-MST-TWD| 30.4 / 34.1| 69.1 / 22.1| 77.5 / 109.3  | 97.0 / 126.4  | 179.1 / 254.0  | 1755.3 / 83.7  |\n| Alg.2| **10.1 / 4.8**| **22.4 / 3.4**| **37.2 / 19.4**| **46.6 / 26.6**| **50.9 / 93.7**| **489.4 / 45.6**|\n\nRegarding graph-based distances, from a conceptual standpoint, using other graph-based distances is possible, but there are important considerations:\n\n- Most graph-based distances are designed to compare node positions or global connectivity in general graphs, not to measure how distributions are supported on a tree or reveal explicit nested groupings\n- Gromov-Wasserstein distance could, in principle, be integrated into an iterative scheme as it compares distributions across different metric spaces or graphs. However, such approaches are much more computationally intensive\n- Another possible direction is to construct a graph from high-dimensional data and transform this general graph into unrolling trees, as explored in recent works such as Chuang \\& Jegelka (2022), Rauchwerger et al. (2025), and Vasileiou et al. (2025). Further investigation would be needed to determine if and how such transformations can be integrated within our iterative learning process\n\nWe will add this discussion as a potential extension for future work to the revision. We appreciate the reviewer’s suggestion to explore broader classes of hierarchical metrics and would be happy to consider specific alternatives if the reviewer has recommendations.\n\n---\nWe thank the reviewer again for their valuable comments and questions. We look forward to discussing any ongoing concerns and suggestions further."}, {"review_id": "l27l2zwmeS", "reviewer": "Reviewer_DBYz", "summary": "This paper presents an unsupervised method that jointly learns hierarchical representations of samples and features using the Tree Wasserstein Distance (TWD). The approach is shown to perform well across a variety of datasets, demonstrating its general applicability. Additionally, the authors provide a theoretical analysis proving the convergence of the proposed method.", "questions": "While the method is generally well-motivated, the paper would benefit from a more detailed explanation of its applicability to specific data domains. In particular:\n\nWhat is the intuition or motivation for applying this method to document and scRNA-seq datasets and how many features and samples for these datasets?\n\nWhat downstream analyses or tasks could benefit from the learned hierarchical representations? For example, can these representations improve clustering, trajectory inference, or biomarker discovery in the scRNA-seq context?\n\nProviding more domain-specific insights and potential applications would enhance the clarity and impact of the work.", "limitations": "While the proposed method is technically sound, the paper would benefit from a clearer discussion of the practical implications of the learned hierarchical structures. Specifically, it would be helpful to elaborate on what these hierarchies represent in various real-world datasets, and why they are expected to be meaningful or useful. Additionally, outlining potential downstream applications—such as improved interpretability, clustering, or other analytical tasks—could further strengthen the motivation and demonstrate the broader utility of the approach.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\n\nThe paper addresses an important and interesting problem—learning latent hierarchical structures in both data samples and features. The methodology is clearly presented, and the authors provide a sound theoretical proof of convergence. The effectiveness of the proposed algorithm is demonstrated through a diverse set of experiments, including classification and prediction tasks across multiple data types such as word-document corpora, single-cell RNA-sequencing data, and graph datasets. The method is also thoroughly compared against a range of existing approaches.\n\nWeakness:\nThe paper would benefit from a clearer explanation of the motivation for applying this method to specific datasets, such as documents data and scRNA-seq data. In particular, it is unclear what hierarchical structure in samples and features means in the context of scRNA-seq, and why this structure is expected to be meaningful or beneficial. Providing more domain-specific insights or interpretations would help justify the applicability and impact of the method in these settings.", "was_revised": true, "final_justification": "I will keep my score as accept", "rebuttal": "We thank the reviewer for their insightful review and for recognising the soundness of our theoretical results and the significance of the problem we address. In terms of the weaknesses and questions:\n\n>**Motivation and hierarchical structure in domain-specific applications**\n\nWe thank the reviewer for this helpful suggestion. We agree that the paper would benefit from a more explicit discussion of why our method is well-suited to document and scRNA-seq datasets. In document data, the features (words) exhibit hierarchical organisation through semantic groupings (e.g., \"animal\" → \"mammal\" → \"dog\"), while the samples (documents) are structured by topics and subtopics (e.g., science → biology → genomics). Representing these hierarchical relationships can improve representation learning in unsupervised settings. We evaluate our method on four word-document benchmarks used in the word mover’s distance and TWD literature (dataset sizes are reported in Appendix E): (i) the BBCSPORT dataset, consisting of 13,243 bags of words (BOW) and 517 articles, (ii) the TWITTER dataset, comprising 6,344 BOW and 2,176 tweets, (iii) the CLASSIC dataset, including 24,277 BOW and 4,965 academic papers, and (iv) the AMAZON dataset, containing 42,063 BOW and 5,600 reviews.\n    \nIn scRNA-seq data, the features are genes, which are functionally related through biological pathways and gene ontologies, forming hierarchical groupings (e.g., immune response genes → cytokine genes → specific interleukins). The samples, i.e., individual cells, follow developmental trajectories or cell-type taxonomies (e.g., hematopoietic stem cell → progenitor cell → mature blood cell types), which are modelled as hierarchical structures in computational biology. For example, in the ZEISEL dataset, which includes 4,000 marker genes across 3,005 single cells from the mouse cortex and hippocampus, hierarchical gene programs underlie neuronal specialisation and brain region identity. In the CBMC dataset, which contains 500 gene markers and 8,617 cells from cord blood mononuclear cells, gene modules related to immune function form structured groupings, while the cells follow a developmental hierarchy from cell types. Representing these hierarchical structures can provide insights into lineage progression, cell type classification, and gene module activation across conditions.\n\nOur method jointly learns hierarchies across both modes, samples and features, making it suitable for these domains. We will incorporate this domain-specific motivation into the revised introduction. We thank the reviewer again for encouraging this clarification.\n\n>**Downstream analyses benefiting from the learned hierarchical representations**\n\nWe thank the reviewer for this helpful suggestion. We agree that explicitly outlining downstream applications would further strengthen the motivation and demonstrate the utility of our method. The hierarchical representations learned by our framework support a variety of downstream tasks. In Section 6, we demonstrated improvements in document and cell classification, as well as link prediction and node classification in hierarchical graph data. Additionally, in Appendix F, we reported results on clustering performance on gene expression data, showing that our method provides meaningful groupings in biological domains.\n\nFrom an application perspective, we believe our method is relevant for biological and semantic data. In scRNA-seq analysis, the sample hierarchy represents biological structure such as cell types or subpopulations, often supporting improved clustering compared to flat representations. In document data, the feature trees over words and the sample trees over documents capture semantic and topical organisation, which is valuable for tasks like topic modelling and document classification. In addition, the feature hierarchy learned from gene expression data groups genes by co-expression and functional relevance. This supports tasks like gene module discovery and potentially biomarker identification at multiple levels of biological resolution. In biomedical applications, we learned that such hierarchical representation could also facilitate visualisation and hypothesis generation.\n\nOur method is also directly applicable to hierarchical graph datasets, where it can be used as a preprocessing step to initialise HGCNs. We demonstrate this on several benchmark datasets: the DISEASE dataset, which simulates the SIR disease spreading process with infection status and susceptibility features; the AIRPORT dataset, representing flight networks between airports; the CORA citation network, where nodes represent research papers classified into machine learning categories; and the PUBMED citation network, which includes three classes of medicine publications. In these settings, initialising HGCNs with our learned hierarchical representation leads to consistent performance gains for graph learning tasks, suggesting the utility of our method for hyperbolic graph learning.\n\nBeyond these domains, our approach could be beneficial in several other areas. For example, in social network analysis, users (samples) and attributes or relationships (features) often display group or role hierarchies; our method supports community detection. In chemoinformatics and drug discovery, both molecules and their substructural or chemical property features can reflect hierarchical organisations, assisting in clustering and identification. In education analytics, student data and knowledge components, such as skills or assessment questions, could exhibit hierarchical structure as well, enabling adaptive curriculum design and hierarchical clustering of learners or resources.\n\nIn the context of trajectory inference in scRNA-seq, our method produces binary trees. While internal nodes do not always correspond to domain-driven groupings, we observe that the primary branches near the root can sometimes align with prominent structures in the data. In the case of scRNA-seq, this suggests a possible connection to major developmental stages or lineage bifurcations, though a more systematic evaluation would be needed to confirm such alignment. Investigating the biological interpretability of these coarse-scale branches remains an interesting direction for future work. Although our method is not explicitly designed for trajectory inference, the hierarchical structure it produces could potentially serve as a meaningful scaffold for reconstructing developmental processes, without requiring additional modelling assumptions. These observations also point to a promising future direction: relaxing the binary tree constraint and allowing more flexible branching, for example, using Haar coefficients to adjust the number of internal nodes automatically. This could improve the interpretability of finer internal groupings.\n\n    \nWe will incorporate these domain-specific applications and future directions into the revised manuscript. We thank the reviewer again for encouraging us to include the broader relevance of our method.\n\n---\n\nWe thank the reviewer again for their helpful comments and hope to discuss any ongoing concerns or suggestions further."}, {"review_id": "L5PcixkDR9", "reviewer": "Reviewer_GK8p", "summary": "The paper proposes a novel approach to simultaneously modeling hierarchies in both samples and features. The joint learning process facilitates the convergence to the shared optimal point, which is theoretically proved and discussed in the paper, leading to an accurate and efficient output. The alternating algorithm, which iteratively constructs trees for one mode and leverages TWD to refine the other applies the joint learning criteria to the field of representation learning successfully.. This bidirectional refinement captures the interplay between data modes effectively, a critical insight for high-dimensional datasets like word-document matrices and single-cell RNA-sequencing data.", "questions": "1. The writing around the citations should be improved. The text should remain correctly readable in case of the removal of the citations (e.g., on line 94, , the method recently proposed in [34] computes multiscale diffusion densities...).\n\n2. The algorithms should be presented better. The current representation is concise, but the level of description is very limited. The algorithms can also benefit from an informative title, instead of just a number.\n\n3. \"We provide detailed descriptions of our method’s implementation within the paper and report the hyperparameters in App. E. Our source code will be made publicly available on GitHub upon publication.\", it would be better if the code and the instructions to reproduce the results were available for the review. It could be easily provided by using an anonymous GitHub link using services like https://anonymous.4open.science/ .\n\n4. The trees and TWDs produce hierarchical structures, but how interpretable are they to end users? This investigation could further improve the paper, especially with regards to the application in classificaiton. \n\n5. It would be interesting to see the visualization of the different levels of the trees to visually investigate how they are extracted and refined in each step or level.", "limitations": "The main limitations of the paper is the lack of the implementation and code to reproduce the results. The investigation of the structure in different hierarchies and the explaianbility of the hierarchical structure is also expected.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper’s primary strength is the novel approach to simultaneously modeling hierarchies in both samples and features.\n\nhe proposed method seems very sound. The theoretical analysis of convergence is a standout feature. The authors provide a robust proof that the iterative process stabilizes, providing confidence to the method’s reliability. There are many theoretical and empirical analyses presented in the Appendix which shed light further onto the different properties of the proposed method.\n\nEmpirically, the experimental results show that the proposed method can be successfully applied to different fields such as document and single-cell classification, link predication and node classification.\n\nOn the other hand, no implementation or code in order to reproduce and validate the results are included, which is expected from a paper of this type.", "was_revised": true, "final_justification": "The author rebuttals and other reviews have not changed my perception that this an very good to excellent contribution to the field. I have, however, raised my confidence.", "rebuttal": "We thank the reviewer for their thorough and valuable review, and for recognising the novelty of our work for simultaneously modelling hierarchies in both samples and features. In terms of the weaknesses and questions: \n\n>**Improving citation-dependent writing**\n\nWe thank the reviewer for highlighting the need to improve the readability of the manuscript around citations.  We will revise the manuscript to make sure that the text remains clear and grammatically correct even if the citation markers are removed. This revision will be applied consistently throughout the manuscript to improve readability and maintain clarity independently of the citation placement. Thank you for pointing it out. \n\n>**Improving algorithm presentation and descriptions**\n\nThank you for your helpful feedback. We agree that the current algorithm presentation is concise but could benefit from clearer structure and more descriptive detail. In the revised version, we will:\n\n- Add informative titles to each algorithm to clarify its purpose (e.g., \"Iterative Joint Hierarchical Representation Learning via Tree-Wasserstein Distance\" instead of a generic \"Algorithm 1\" and \"Haar Filtering over Alternating Tree Refinement\"  instead of \"Algorithm 2\")\n- Expand the descriptions in each line within the algorithm boxes, making key steps more explicit and easier to follow (e.g., construct an initial sample and feature trees, compute sample and pairwise TWD matrices,  Tree Haar Wavelet filtering, etc)\n\nWe hope that these changes will improve clarity without sacrificing conciseness, and make the algorithmic components more accessible and self-contained. If there are any additional specific suggestions for improving the algorithm presentation, we would be happy to incorporate them.\n\n>**Code availability for reproducibility**\n\nWe thank the reviewer for highlighting the importance of reproducibility. We fully agree and are currently finalising and cleaning our code to ensure it is clear and easy to use for the community. The full implementation will be released, and we will include a link in the final version of the paper. Additionally, we have provided detailed implementation information in Appendix E, including model settings, training procedures, and hyperparameters, allowing one to fully reproduce our results. If there are any specific details that are missing or unclear, we would be happy to clarify them.  We apologise for any inconvenience caused by releasing our code upon publication, and appreciate the reviewer’s understanding.\n\n>**On the interpretability of learned trees and TWDs for end users**\n\nWe thank the reviewer for raising this important aspect. This is indeed a valuable question. Our motivation for applying the method to document and scRNA-seq datasets because the hierarchical structures present in both domains. In document data, the features (words) exhibit hierarchical organisation through semantic groupings (e.g., \"animal\" → \"mammal\" → \"dog\"), while the samples (documents) are structured by topics and subtopics (e.g., science → biology → genomics). In scRNA-seq data, the features are genes, which are functionally related through biological pathways and gene ontologies, forming hierarchical groupings (e.g., immune response genes → cytokine genes → specific interleukins). The samples, i.e., individual cells, follow developmental trajectories or cell-type taxonomies (e.g., hematopoietic stem cell → progenitor cell → mature blood cell types), which are modelled as hierarchical structures in computational biology.\n\n\nOur method produces binary trees. While internal nodes do not always correspond to domain-driven groupings, we observe that the primary branches near the root can sometimes align with prominent structures in the data. In the case of scRNA-seq, this suggests a possible connection to major developmental stages or lineage bifurcations, though a more systematic evaluation would be needed to confirm such alignment. Investigating the biological interpretability of these branches remains an interesting direction for future work. In applications such as clustering or classification, we find that the coarse structure of the binary tree can highlight major distinctions in the data. These observations also point to a promising future direction: relaxing the binary tree constraint and allowing more flexible branching, for example, using Haar coefficients to automatically adjust the number of internal nodes. This could improve the interpretability of finer internal groupings. \n    \n    \nIn terms of the TWD, its strength lies in its transparent and decomposable structure. Given two distributions (e.g., feature or sample vectors), the TWD quantifies the minimal \"mass transport\" needed to align them over the tree, with transport cost reflecting branch length. This formulation enables users to understand where and at what scale differences occur. For instance, if most of the transport cost is concentrated at high-level branches, the two distributions differ in broad terms; if the cost is concentrated near the leaves, they are largely similar but differ in fine details. As a result, TWD not only quantifies similarity but also supports qualitative interpretation: users can trace which parts of the tree contribute most to the distance between any two data points. In classification tasks, this interpretability reveals which regions of the hierarchy drive class separation and can help identify key transitions in the data. \n    \n    \nTo further highlight this interpretability, we will include in the revised manuscript visual breakdowns of TWD computation and learned tree structures for real datasets, demonstrating how the hierarchy enables multiscale reasoning and class separation. We hope these visualisations will help clarify how particular branches, at varying scales, relate to potential biological or functional groupings and will give users insight into the decision boundaries and similarity patterns uncovered by the model.\n\n>**On visualising tree levels and their refinement**\n\nWe thank the reviewer for this valuable comment regarding the visualisation of different levels of the learned trees and how they evolve during each iteration in our method. We agree that stepwise visualisations would offer valuable insight into the process of hierarchical refinement. Thank you for this suggestion! Since, at this stage, we are unable to include new figures in the PDF, we provide here a detailed textual description of the tree refinement throughout the iterative process, based on the toy video recommendation example in Appendix F.1.\n\nThe iterative refinement process starts with an initial tree structure, which is unstructured since the preliminary grouping is based on initial data relationships ($\\mathbf{M}_c$ or $\\mathbf{M}_r$). At this point, major clusters are not well formed, and samples or features from different categories are mixed. After the first iteration, the algorithm uses the current tree from one mode (e.g., users) to compute a TWD for the other mode (e.g., videos). This leads to a newly induced tree that starts to separate the major categories, and major branches start to reflect groupings such as fiction and documentary in videos, but subcategory assignments are not fully sorted. As the algorithm proceeds through subsequent iterations, it alternately refines the trees for users and videos. The updated tree for one mode informs a more accurate computation of TWD in the other mode, enabling finer divisions within the existing branches. Subcategories like action, drama, and sci-fi (for fiction videos) become more distinct and better grouped together. After several iterations, the subcategories become more cohesive and consistently organised, and the tree structures stabilise. \n\nIn our revised manuscript, we will include these visualisations that illustrate how the trees evolve with each iteration and point out key groupings formed in the learned trees. We hope these new figures will make it possible to directly observe how broad groupings emerge in early iterations and how finer subcategory assignments are progressively refined, ultimately converging to clearly organised hierarchies.\n\n---\n\nWe thank the reviewer again for their insightful comments and questions. We look forward to discussing any ongoing concerns and suggestions further."}]

Author rebuttal:

Meta-review: This manuscript proposes a methodology for learning hierarchical representations of data along both features and samples. In addition to motivating and introducing the idea, the authors provide theoretical convergence results for their method and an empirical study showing improvement in a variety of learning tasks.

This work was received in extremely regard by the reviewers. In particular, the reviewers were impressed by the paper's problem setting -- learning hierarchical representations for both features and samples -- and equally impressed by the design and performance of the proposed methodology.

Due to the novelty of this setting, and the universally positive reception of this work, I believe this paper is appropriate for spotlight presentation or higher.
Decision: Accept (spotlight)