MS-IMAP - A Multi-Scale Graph Embedding Approach for Interpretable Manifold Learning
Authors: Shay Deutsch
Lionel Yelibi
Alex Tong Lin
Arjun Ravi Kannan
Venue: Accepted by TMLR ()

Abstract: Deriving meaningful representations from complex, high-dimensional data in unsupervised settings is crucial across diverse machine learning applications. This paper introduces a framework for multi-scale graph network embedding based on spectral graph wavelets that employs a contrastive learning approach. We theoretically show that in Paley-Wiener spaces on combinatorial graphs, the spectral graph wavelets operator provides greater flexibility and control over smoothness compared to the Laplacian operator, motivating our approach. A key advantage of the proposed embedding is its ability to establish a correspondence between the embedding and input feature spaces, enabling the derivation of feature importance. We validate the effectiveness of our graph embedding framework on multiple public datasets across various downstream tasks, including clustering and unsupervised feature importance.

Published in Transactions on Machine Learning Research (04/2026)
MS-IMAP: A Multi-Scale Graph Embedding Approach for Interpretable Manifold Learning
Shay Deutsch** *
shaydeu @ymail.com
Lionel Yclibits
ylblio001 @emjuct..cccc
Alex Tong Lint8
atlin271 @gmmil com
Arjun Ravi Kannan
arjun.kenmann@@mail.comm
Reviewed on OpenReview: https: // openreview. net/ form? id=pc6BgürCjj
Abstract
Deriving meaningful representations from complex, high-dimensional data in unsupervised scttings crucial across diverse machine learning applications. This paper introduces framework for multi-scale graph network embedding based on spectral graph wavelets that employs contrastive learning approach. We theoretically show that in Paley-Wiend spaces on combinatorial graphs, the spectral graph wavelets operator provides greater flexibility and control over smoothness compared to the Laplacian operator, motivating our approach. A key advantage of the proposed embedding is its ability to establish a correspondence between the embedding and input feature spaces, enabling the derivation of feature importance. We validate the offectiveness of our graph embodding framework on multiple public datasots across various downstream tasks, including clustering and unsupervised feature importance.
1 Introduction
Graph Embeddings and Manifold Learning (Roweis & Saul (2000); Belkin & Niyogi (2003); Tenenbaum et al. (2000); der Maaten & Hinton (2008); Coifman et al. (2005)) play pivotal role in analyzing complex data structures encountered in wide range of machine learning applications. The representations learned by these techniques uncover the intrinsic geometry of the data and support downstream tasks such as clustering and ietalizaaiionsspcciiale
scenarios where labels are unavailable or unreliable. Most manifold learning methods typically perform nonlinear dimensionality reduction to embed data into lower-dimensionall spaces, while some also explore high-dimensional embeddings (Gama al. (2019); Sun et al. (2009)). Nonetheless, most such methods lack a direct connection to the input features, preventing the evaluation of individual feature contributions and resulting in embeddings with limited interpretability. This limitation is critical in domains like finance, where key behavioral indicators must be traced back to individual features, or in biology, where interpretability is vital for uncovering mechanisms within complex systems. To address these limitations, this paper makes three primary contributions:
We introduce framework that leverages multi-scale graph representation using 3D tensor and contrastive learning. This enhances the expressiveness of the embedding by jointly optimizing over both low- and high-frequency components, enabling the capture of fine-grained data structure.
We propose manifold-based interpretation method for fcature importance, allowing each dimension of the learned embedding to be directly linked back to unique original feature.
*Corresponding author. Equual second authors. 9hhis work was conducted Discover Financial Services.

---

Published in Transactions on Machine Learning Research (04/2026)
4
MS-IMAP Regularized Embedding 42
Most Important features
K-means Clustering
Festure Importance
NS-IMAP Embrdding &
(a) Data X
(b) Encoding using SGW 4x
(c) Regularized Embedding x
Figure 1: Ilusttationn of the proposed method. (Left) Initial encoding using SGW in (b) and regularized embedding in (c). (Right) Feature embedding space representation, which has a one-to-one correspondence with the input features, can be used along with clustering for interpretability models, such as feature im- portance.
We characterize the theoretical representational capacity of the Spectral Graph Wavelet (SGW) (Hammond et al. (2011)) operator by analyzing functions from Paley-Wiener spaces (Pesenson (2008)) on combinatorial graphs. This analysis demonstrates that SGW provides greater flexibility and smoother control compared to traditional Laplacian-base operators.
Our method, MS-IMAP (Multi-Scale Interpretable Manifold Approximation and Projection), is designed to associate each embedding coordinate with a unique input feature, thereby enabling the estimation of feature importance relative to the learned embedding space. This interpretability arises from a deliberate design choice: each embodding coordinate corresponds to a single, unique input feature and is optimized through spectral graph wavelets to align with the underlying graph structure. Unlike standard manifold learning and graph embedding methods, where each embedding dimension typically blends information from all input features, MS-IMAP preserves the semantic identity of each feature by assigning it to distinct embedding dimension.
Another key distinction of our approach is how it uses the graph Laplacian. Classical methods such as Laplacian Eigenmaps (Belkin & Niyogi (2003)) and many modern techniques such as UMAP (Mclnnes et al. (2018)) locus mainly on low-frequency components of the Laplacian and therefore emphasize coarse structural information. In contrast, MS-IMAP relies on multi-scale graph representations via SGW within contrastive learning framework. This design allows the embedding to retain the original dimensionality while still preserving clear feature-wise interpretation.
Optimization is performed with stochastic gradient descent (SGD) on a 3D tensor based contrastive objec- tive, so that the learned embeddings capture both local and global graph structure. Empirically, we observe that applying standard feature importance criteria MS-IMAP embeddings yields more informative feature subsets for clustering than applying the same criteria directly to the raw input features. This indicates a stronger semantic alignment between the embedding and the underlying data structure, which is often lack- ing in nonlinear deep models conventional manifold-learingg methods. MS-IMAP achieves interpretable and semantically faithful embeddings, while remaining competitive with state-of-the-ar graph embedding techniques across a range of tasks.
2 Related Work
Extensive efforts have been dedicated to exploring manifold learning methods that aim to produce low- dimensional embeddings while preserving the intrinsic structure of high-dimensionall data. Classical manifold learning techniques include Locally Linear Embedding (LLE) (Roweis & Saul (2000)), Laplacian Eigenmaps (Belkin & Niyogi (2003)), Isomap (Tenenbaum et. al. (2000)), Diffusion Maps (Coifman al. (2005)), and
2

---

Published in Transactions on Machine Learning Research (04/2026)
t-distributed Stochastic Neighbor Embedding (t-SNE) (van der Maaten & Hinton (2008)). One of the central challenges in these methods is maintaining a balance between local and global structure, especially in the presence of noise, To address this, various denoising strategies have been developed (Hein & Maier, 2006; Deutsch & Modioni, 2017; Deutsch et al., 2018; Deutsch & Modioni, 2015a;b), aiming to improve robustness and extend the applicability of these techniques to real-world, noisy datasets.
UMAP (Melnnes et al. (2018)) has shown strong empirical performance for visualization and cluster separa- tion. It typically employs Laplacian Eigenmaps initialization, which emphasizes low- fwwereeuency graph signals while neglecting freguery patterns that may carry important information.
Recent work (Damrich & Hamprecht (2021)) has shown that negative and positive sampling in UMAP significantly influences the effective loss function. These insights have led lo studies on alternative sampling schemes (Yan et al. (2023)), bounds on generalization error in embedding methods (Suzuki et al. (2021)), and biascs introduced by high-degree nodes during graph construction and sampling (Kojaku et al. (2021)).
Manifold learning and non-linear dimensionality reduction serve purposes beyond visualization; one of their key goals is to learn representations that capture intrinsic geometric structure and approximate geodesic distances. The Icarned embedding space dimensionality can go beyond two or three dimensions. However, increasing dimensionality does not necessarily improve robustness: noise and the curse of dimensionality can introduce instability. These methods also suffer from computational limitations. Moreover, to compute em- bedding in higher dimensions many algorithms, such as Laplacian Eigenmaps and UMAP require computing large portions of the Laplacian spectrum, making them computationally intensive for large datasets. Our method addresses these issues by producing interpretable and efficient. embeddings while leveraging spectral information across multiple frequency scales.
A key limitation of many manifold learning and graph embedding methods is the absence of an explicit and interpretable mapping, which prevents evaluating the contribution of each input feature to the resulting low-dimensional embedding. This gap restricts their utility in domains where interpretability is crucial. Providing such a link can enable improved feature selection, clearer interpretation of embeddings, and greater robustness to noise (Hosseini & Hammer (2020); Frye et al. (2021); Zharmagambetor & Carreira-Perpinan (2022)).
SGW are central to our method (see Section 3.1), offering localization in both spectral and vertex domains. SGW apply a family of bandpass filters to the graph Laplacian's spectrum, allowing for multi-resolutionn analysis of graph signals and providing richer representations across frequency bands. Our approach is also related to graph embedding methods that exploit multiscale structure, such as Diffusion Wavelets (Ronald R. Coifman (2006)), Geometric Scattering for Graph Data Analysis (Gao et al. (2018)), and the Graph Scattering Transform (Gama et al. (2019)), which extract graph features for downstream tasks like classification. These methods often embed data in higher-dimenn dimeessoinni spaces to improve representation power.
2.1 Preliminaries
Consider a set of points x = {xn} , = 1,,.,N, E ERD, where each Xn E R,D, sampled from an unknown smooth manifold M C RDD To model the geometric structure of the data, we construct an undirected, weighted graph G = (V,W) over the point set x, where V is the set of nodes and W = =(wm) is the weighted adjacency matrix encoding edge weights between nodes. The weights Wam can be computed using various techniques. A common approach is to apply a Gaussian kernel to pairs of points x,, and Xm within the k nearest neighbors of xn, denoted as kNN(x,). In this work, we adopt the fuzzy graph in the input exports
space using a smooth exponential kernel (as in McInnes et al. (2018) with local scaling parameter per point, though other graph construction methods may also be employed. Let L denote the combinatorial graph Laplacian defined as L = D W, where D is the diagonal degree matrix with entries dn.. d(n). and the degree d(n)) of vertex n is the sum of the weights of all edges incident to n. We use the normalized which has real eigenvalues in
graph Laplacian, defined as: LN = D D-1/2DD 1/2 = I D the interval [0, 2). For simplicity, we refer to Ly as L for the remainder of the paper. Let Aj...., ; An and PI....,yy denote the eigenvalues and corresponding eigenvectors of L. We denote the matrix of eigenvectors Φ N and the diagonal matrix of eigenvalues as A = diag(A1,-... AN). In our setting, each input
3

---

Published in Transactions on Machine Learning Research (04/2026)
coordinate fa R. viewed as graph signal, is defined across the N nodes of the graph. The graph Fourier transform (GFT) f of the graph signal f is defined as the expansion of f in terms of the eigenvectors & of the graph Laplacian, that for frequency X we have: f():) f
The original data matrix x can thus be expressed as a collection of such signals: x = f2,.,, Sa, (a))..)) Our goal is to derive multi-scale embedding from these signals that approximates the intrinsic manifold coordinates, using spectral information derived from the graph Laplacian.
2.2 Multi-scale representations using SGW
In the past two decades, several multi-scale representations for data on irregular graphs have been devel- oped, including Spectral Graph Wavelets (SGW) (Hammond et al. (2011)) and Diffusion Wavelets (Ronald R. Coifmman (2006)). In this work, we focus on the SGW-based multi-scale graph transform. SGW offers a principled way to balance spectral and spatial resolution, as its coefficients are localized in both domains. These wavelets are constructed using kernel function operator g(L), which acts on a signal by modulating each of its spectral components (i.e., Fourier modes). This design enables a trade-off between vertex-domain (spatial) and frequency-domain (spectral) localization. The spatial localization is implicitly controlled by a single scale parameter defined in the spectral domain: greater localization in the vertex domain corresponds to a broader spectral bandwidth. To represent a signal f E TR at multiple scales S = [81, 12,,,82,-0 .2,..., .8]], the SGW transform is defined as follows. Let g(A) be a bandpass filter in the spectral domain. Let 8n E RN denole the delta function centered at verlex E G, where 8, (m) = 1 if m = n, and 6n (m) = 0 otherwise. A wavelet centered at node and scale is then given by: (s,,n) n) Dg(s>) 8n. The value of vi n) at vertex m can be written explicitly as: (bss,) (m) = E (( (n) Given a graph signal f E RN the SGW W coefficient at node n and scale defined as:
N (s, - ((())))(()))((())
(1)
Each scale corresponds to a distinct spectral band whose contribution is determined by the SGW filter's response g(s)e) over the Laplacian spectrum, which is adaptive to the graph structure. Fast computation using Chebyshev polynomials: Direct computation of SGW coefficients is computationally expensive, requiring O(N3) operations for N nodes. To address this, SGW proposed an efficient algorithm based on approximating the scaled generating kernels g(s)) using low-order Chebyshev polynomials of the Laplacian L, which are applied directly to the input signal. This significantly reduces computational cost (see Appendix for further details).
3 Our Proposed framework: Multi-Scale IMAP
In this section, we introduce Multi-Scale IMAP (MS-LMAP), a framework for interpretable manifold-basee embedding that leverages multi-scale graph representation. The method defines differentiable mapping h : + t/x supported on discrete graph G = (V, W), where vix denotes multi-scale graph transform of x. MS-IMAP consists of two main steps.
Step 1 (multi-scale representation). We construct a multi-scale graph representation by applying the Spectral Graph Wavelet (SGW) transform to feature signals across multiple scales and graph frequencies. Coneretely, we form 3D tensor of SGW coefficients with dimensions K N x D (number of scales x nodes x features), as detailed in Algorithm 1. This tensor encodes the multi-scale transform in a way that allows to align the optimized embedding with all scales simultaneously and to enforce differentiable structure on the transformed features.
The 3D tensor representation is beneficial for two main reasons. First, it enables joint optimization across scales: using the SGW coefficients at all spectral bands within the contrastive loss, we smooth low- and high-frequency information simultancously over a fixed graph structure. This reduces the need for ad hoc thresholding of high-frequency bands and ensures that all scales are optimized consistently for each pair of nodes. Although neighboring scales may have correlated wavelets, each scale corresponds to distinct
4

---

Published in Transactions on Machine Learning Research (04/2026)
Algorithm 1: Encoding multi-scale structure in the optimized embeddings
Input: Set of points J Initial Node Embeddings 4xx
Output:
Step 1: Construct G (V, W) from S
Step 2: Construct the Laplacian L.
Step 3: Compute Amax (L).
Step 4: Initial Embedding construction:
For d d=1,., associated with the feature signals {fa) do:
Compute f (s).:) :) at scales 8k: k I...K using Chebyshev approximation.
Step 5: Concatenate ( for d = 1, .... D in matrix form (*x,,) for each fixed scale 8k1 representing all scales using 3D tensor e/' E RK X XRN.
bandwidth and captures graph structure at different neighborhood radius, so combining them yields representation that balances local and global structure. Second, the tensor is organized feature-wise, each embedding coordinate is aligned with single input feature. This one-to-one correspondence is what underpins the interprotability results Section 4, where we use the embedding to derive feature-importanco rankings.
Step 2 (embedding optimization). We then integrate the multi-scale representation into an SGD- optimization, using 3D-tensor contrastive objective that leverages both low- and high-frequency informa- tion. This optimization refines the multi-scale coefficients into an embedding that is jointly regularized by the graph structure and remains aligned with the original features, leading to fine-grained manifold regularization and improved robustness in downstream tasks.
3.1 Step 1: Feature Representation to encode multi-scale structure
To encode multi-scale representations used for subsequent optimization we focus on a 3D tensor based method. Note that each dimension in the embedding space is constructed using a single feature in the original feature set, which is an essential characteristic that can leveraged for interpreting graph embeddings. 3D-Tensor Encoding: The 3D-Tensor Encoding method concatenates the multi-scale representation for all features at fixed scale, generating a matrix representation ;s) RD x R for each scale & encoding SGW coefficients across features and nodes. After concatenating all fa at fixed scale 8k, we have
( S (Sk: (ll: (Sk. :))x::))l-)) )11-
where we designate the concatenation using ll. with c(1b,, (sk, V,,, (sk. = Vy (sk. :I Wy, (sk, :::. The optimization can be performed separately for each scale Sk. or jointly for all scales Sy using the 3D tensor vix E R TRD RV.. We detail the proposed feature representation encoding methods the pseudo code algorithm 1.
3.2 Step 2: Optimization design using multi-scale network structure.
We propose use optimization basod on stochastic gradient descent (SGD) that begins with the initial embedding within the spectral graph domain. While the initial spectral graph wavelet (SGW) representation already encodes multi-scale information, it does not yet explicitly incorporate the graph-based structural loss used by MS-IMAP. In this step, we refine the initial embedding by optimizing a multi-scale structural objective across the graph.
3D-tensor based optimization for graph embedding. Encoding the representation in a 3D tensor (1) it naturally captures the multi-scale structure induced by SGW; and (2) it makes it possible to define a scale- wise structural loss that can be optimized jointly across all nodes and features.
offers two key advantages:
5

---

Published in Transactions on Machine Learning Research (04/2026)
We write the SGW tensor as
E E (1....,KK, K}, de (I...., D). n (V).
where m index nodes, d indexes features, and &o indexes SGW scales. Equivalently, for each scale Sk, have slice ( E RDX,, whose d-th row, (x d,:) RN, contains the coefficients of feature d across nodes, and whose n-th column, (sk,,, E R", contains the coefficients of all features at node n. An SGD update rule at iteration t, applied to the 3D tensor, is given by
(+tt))) )kdd = (W )add f(1/x)k,d.
(2)
where a is the learning rate and L is the structural loss defined below, and )h xxx for t t= 0. We employ cross-entropy loss that compares the empirical graph weights Wnm to a similarity term derived from the SGW representation:
K
C(Ux W) =
Wnm log
W Ware
+ (1 Wmm) log
Wnm
(3)
Upm (5k)
Unm (5k)
n,m==
k=
where
V
r0m
(8k)
(s,, n) 1 (sk, m)ll
(4)
Thus, eiin (5k) encodes structural similarity between nodes n and m at scale Sk, based on the SGW coefficients of all features.
We apply the optimization based on SGD with respect each scale sk, where the gradient of the loss with respect to the SGW representation node can be written in the form
ac (l((x(5kii)))) ()
Σ Wnen uye (8k) F'rm (8k)
(
(5)
O1 (1kx(88k;,1)))
Umm (Sk) (nm((s))
rnm (sk) 1|
where
Tom (5x) ('s( (s)...m)
(6)
This expression makes explicit that the gradient at node n depends on SGW differences Fum (8x) across neighbors m and across all features.
Feature updates. Our implementation employs an optimization process that operates jointly across all scales 8& and all nodes, and then aggregates the optimized embeddings with respect to each scale. We define the final embedding tensor as
1x
Σ (st,),);
so that the dimensionality of the final embedding is dx E DDN maintaining a one-to-one correspondence between the D embedding coordinates and the D input features. For each feature index d, the d-th embed- ding coordinate is aligned with input feature and is updated in aligure-aligned manner via SGD: the gradient is taken with respect to the entries (* d, n),, while the structural similarity term wxix (sx) and the graph construction depend jointly on all features through (sx(sk,;,))) and ()x(Sk, m). Thus, MS-IMAP enforces feature-wise aligument of embedding coordinates with input features, while still incorporating joint structural coupling through the graph and the structural term. This preserves the one-to-one correspondence that is crucial for interpretability.
Remark 1. While we adopt the SGW formulation by Hammond (2011), our framework compatible with any multi-scale spectral graph transform that yields tensor of coefficients indexed by scales, features,

---

Published in Transactions on Machine Learning Research (04/2026)
and nodes. Alternative filter banks or polynomial approximations of the graph Laplacian can be substituted without changing the optimization structure.
Remark 2. While the manifold hypothesis suggests that the data lie near lover-dimensional manifold, this docs not require the embodding to be explicitly low-dimensional. In MS-IMAP, we retain the origi- nal feature dimensionality to preserve leature-wise interpretability, while still leveraging multi-scale graph structure through the SGW-based optimization. This allows us to avoid repeated large-scale eigenspectrum computations, in classical methods such as Laplacian Eigenmaps or UMAP, and focus computational effort on the multi-scale structural refinement of the embedding.
4 Feature Importance using Explicit Correspondence with Embedding Features
Feature importance is one of the simplest yet most fundamental tools for model interpretability: it associates imput features with model outputs or representations, highlighting which features most strongly influence the learned structure. In MS-IMAP, each embedding coordinate has a one-to-one correspondence with an original input ura This explicit alignment. allows us to apply standard feature-imporlance criteria directly in the embedding space and map the selected dimensions back to specific features.
We show that feature subsets selected from the MS-IMAP embedding lead to superior clustering performance compared to subsets selected from the input space using the same methods. Unlike competing approaches, MS-IMAP retains semantic alignment and supports a direct correspondence between embedding dimensions and original features, making it an interpretable alternative to black-box embeddings. Empirically, apply- ing fepturc-in techniques to the MS-IMAP embedding yields more informative feature subsets for clustering than applying the same techniques to the raw input. This provides evidence that our method pre- serves semantic consistency, a property that traditional nonlinear embeddings typically do not provide. We explore two approaches for estimating feature importance that leverage this explicit correspondence, both in unsupervised self-supervised settings where feature selection is performed without access to ground-truth labels.
Laplacian Score (LS) on embedding features. The first method uses the Laplacian Score (LS) (He et al. (2005)), which assesses the importance of each feature based on how smoothly it varies over a graph constructed from the data. In our setting, the variable corresponding to the d-th embedding coordinate, denoted fa E R.N. is inferred from the coordinate Thyy in the MS-IMAP embedding. Details on how we compute LS with respect the embodding fcatures vxx are provided in Appendix A.5.
Mutuallinfformation-B based feature importance. We also consider self-supervised approach based Mutual Information (MI). First, we apply k-means to the MS-IMAP embedding to obtain clusters, which serve as pseudo-labels. We then compute MI between each embedding feature and these pseudo-labels, analogous to standard MI-based feature selection with true labels. This quantifies how much each feature contributes to cluster assignments and highlights those most informative for distinguishing clusters. The MI-based procedure is summarized in Algorithm 2 (see Appendix A.5).
Comparison: LS vs. MI-based feature importance. Laplacian Score and MI-based feature impor- tance are complementary. LS is fully unsupervised, geometry-driven criterion that favors features which vary smoothly over the graph and align with manifold structure; however, it is not task-aware and can overvalue smoothly varying but irrelevant features. In contrast, MI explicitly measures the statistical depen- dency between each feature and cluster assignments (obtained from k-means on either the raw features or the MS-IMAP embeddings), making it more directly aligned with clustering performance. In our experiments, both criteria become more informative when applied in the MS-IMAP embedding space than in the raw feature space.
5 Experimental Results
We evaluated our approach on both synthetic and real datasets in two application domains. First, we demonstrate the of MS-IMAP embeddings for feature importance, novel application that highlights
7

---

Published in Transactions on Machine Learning Research (04/2026)
Number Produris Kest
vs Namber Rest
Dead Eits
Wurber of Prasuures Kept
Nunte of Peatures Koct
(2.a) Cora dataset
(2.b) Zillions dataset
Figure 2: Clustering results using Adjusted Rand Index (ARI) on Cora and Zillions datasets with subset of features based on feature importance, comparing the MS-IMAP embedding space and the input features.
the interpretability of our method. Next, we assess the performance of MS-IMAP in clustering tasks using the full embedding space, comparing it against a range of manifold learning- based embedding techniques. Datasets We study the performance of MS-IMAP compared to other methods for real datasets. We chose mix of datasets from varied fields: the Census dataset (Dua & Graff (2017)) is financial dataset containing information about individuals extracted from the 1994 US Census; the Zilionis dataset (Zilionis el al. (2019)) is a biological dataset containing single-cell sequencing data from different types of cells; the Animals with Attributes (AWA) (Xian et al. (2017)) dataset is an image dataset containing images of animals; The Human Activity Recoguition (HAR) (Anguita et al. (2013)) dataset which consists of sensor data collected from subjects performing six different activities:; and the Cora dataset (Sen et al. (2008); McCallum et al. (2000)) is a graph dataset of a citation network consisting of scientific publications.
5.1 Assessing Feature Importance via Embedding Correspondence
We demonstrate the application of MS-IMAP for feature importance selection by evaluating clustering per- formance on subsets of the most significant features. Feature selection performed with respect to both the input space and the MS-IMAP embedding space, using two alternative methods: Laplacian Score and Mutual Information-based selection. We compare clustering results obtained by applying k-means to subsets of features selected from the embedding, against those selected directly from the original feature space using the same methods. Unlike most manifold learning and graph embedding techniques, MS-IMAP produces embeddings in which each dimension corresponds directly to an input feature. This one-to-one mapping enables direct assessment of individual feature contributions, preserving interpretability.
Evaluation metrics: Clustering performance is measured using Adjusted Rand Index (ARI) and Adjusted Mutual Information (AMI). In all experiments, k-means is applied in the embedding space. On the Cora and Zillions datasets, Figures 2.a, 2.b, 6.a, and 6.b show that, for both feature selection methods and for both evaluation metrics, clustering accuracy is consistently higher when using the MS-IMAP embedding compared to the original features across varying feature subset sizes. This improvement highlights the obbustness of the MS-IMAP embedding, while maintaining interpretability through its feature-wise alignment. Additional experiments on the AWA and HAR datasets, presented in the Appendix, further confirm these gains. In addition to the ARI and AMI evaluations, we further examine the internal clustering structure induced by the selected feature subsets. As detailed in Appendix A.3, we report the silhouette score computed on these subsets to provide a label-free validation of cluster coherence. This complementary analysis confirms that the subsets selected via MS-IMAP yield more stable and well-separated clusters compared to those from the raw feature space.
5.2 Clustering Output Embedding Space
We evaluate our method by testing its clustering performance using the full output embodding and compar- ing it to several representative methods, including UMAP, I-SNE, Diffusion Maps, TSOMAP, PHATE (Moon
8

---

Published in Transactions on Machine Learning Research (04/2026)
et al. (2019)), PaCMAP (Wang et al. (2021)), TriMap (Amid & Warmuth (2019)) and HeatGeo (Huguet et al. (2023)). While most of these methods perform dimensionality reduction, we emphasize that dimen- sionality reduction is a methodological design choice, not an inherent advantage. Many manifold learning techniques default to low-dimensional embeddings (c.g., 2D or 3D) due to computational constraints such as eigendecompositionn, yet. this often comes at the cost. of reduced expressiveness and, particularly for nonlinear methods, diminished interpretability. However, since most baseline manifold learning methods we compare rely on nonlinear dimensionality reduction, we include clustering results based on the embedding space of the Graph Scattering Transform (Gama et al. (2019)) for completeness. The Graph Scattering Transform, which is primarily used in graph classification and semi-supervised learning tasks, can be adapted for un- supervised learning via simple vector- concatenation. We note that the Graph Scattering Transform offers less interpretability than MS-IMAP, which preserves a one-to-one correspondence between input fea- tures and embedding dimensions, Additionally, the increased dimensionality of the embeddings produced by the Graph Scattering Transform can lead to higher computational costs. As shown in Table 1, MS-IMAP consistently delivers stable and competitive performance, ranking best or second-best across all benchmarks. In contrast, state-of-tie ohett--rtrerttte methods often underperform or show instability at least one dataset.
Runtime and Computational Complexity The computational complexity of Multi-Scale IMAP is of Oog IN)) Experiments were performed on Virtual Server containers with 32 cores Intel Xeon 8259CL running at 2.50Ghz and 256GB of RAM. Our method is practical for handling datasets with millions of points and dozens to hundreds of features within a few hours, making the method well-suited for real-world applications. For example, on the Cancer QC data-set of 48,969 samples and 306 features took 27 mins. Additional details and rumning times on all datasets are provided in the Appendix. Note that the increase in running times compared to methods such as UMAP is primarily due to our extension of manifold learning to include interpretability.
Data
HAR
Census
Zilionis
AWA
Cors
Method
ARI
AMI
ARI
AMI
ARI
AMI
ARI
AMI
ARI
AMI
Features
0.45++0000
0.58+00.0
0.11-0.00
0.15++0.00
0.48-+0.04
0.63++0.01
0.07-00.00
0.16++0.00
0.01-0.02
0.07+00.02
ISOMAP
0.50++0.00
0.65+0000
0.18=0.00
0.09+00.00
0.44++0.03
0.56++0.00
0.52=00.0
0.60++0.00
0.12=0.02
0.17+0.02
DM
0.301.0.00
0.40
0.01
0.02.10.00
0.22.10.00
0.42:10.00
0.06
0.07:10.00
UMAP
0.61110000
0.
0.23=0.00
0.15::0.00
0.52:10.00
0.711:0.00
0.74 10.00
0.8210.00
0.33:00.1
0.36110.2
G-NNE
0.6210.00
0.70.10.00
0.20-0.00
0.24.10.0
0.39-10.01
0.60-10.00
0.71.10.00
0.80-10.00
0.20
0:2410.00
Scattering.T
0:15200000
0.58:0.00
0.12=0.00
0,16::000
0.46:20001
0.612:0002
0.67=0000
0.75:b0.00
0.0620.03
0.1420.08
PHATE
0.35:10.02
0.6620.02
0.20=0.00
0.091:001
0.65:10.01
0.72::0.00
0.68:2003
0.73:10.00
0.12=0.01
0.22:2001
TriMap
0.59+:0.2
0. 0.72=0.01
0.20=0.00
0.14+0000
0.62::0.08
0.77=0.01
0.74=0001
0.81#0001
0.15=0000
0.20=0.00
PaCMAP
0:62::001
0.7522.08
0.24=0000
0,15::0.00
0.35:20.03
0.762:0.00
0.75:00.1
0.81=0.01
0.21=0.01
0.2620.00
HeatGoo
0:60+-0.00
0.89 0.69+0.00
0.15=0.00
0. 0.10+0000
0.65=20.00
0.74+0000
0.17=0.00
0.20++0.00
MS-IMAP
0.67:10.0
0.80
0.23
0.15
0.70
0.76.10.01
0.74:10.00
0.81.10.00
0.33 0.00
0.36
Table I: Clustering results comparison using ARI and AMI on the Census, Zilionis, AWA, and Cora datasets.
Each entry report the mean clustering accuracy and standard deviation. The best performance is bolded.
6 Theoretical results: Sampling set for Smooth Manifolds with functions defined over Paley-Wiener Spaces
In this section, we characterize the theoretical properties of the representation power of the SGW operator by considering functions sampled from the Paley-Wienee spaces Pesenson (2008) combinatorial graphs. Pesenson introduced the Paley-Wiener spaces and demonstrated that Paley-Wiener functions of low type are uniquely determined by their values on certain subsets, known as uniqueness sets U. We show that the SGW operator can represent functions f within the Paley- Wiener space more efficiently than the graph Laplacian operator L. This efficiency is demonstrated by showing that the SGW operator is more effective in representing functions with larger bandwidth w in the Paley- Wiener spaces (i.e., with higher frequencies) using subsets of nodes from the uniqueness sets U. To characterize the representation properties of functions defined over PG)( we employ the notion of the A-set, introduced by Pesenson which is central to our investigation. Formally, the Paley-Wiener Wiener space of w-bandlimited signals is defined as PG.
{ {1l (())()) = A w
9

---

Published in Transactions on Machine Learning Research (04/2026)
The space Lz) is defined as the Hilbert space of all complex-valuee functions, and (y(S) is defined as the space of all functions from [G) with support in S: [2((S) = (o E [z2(G) (p(r) 0, E ((G) S}. The A-set is defined as follows: a set of vertices S C V is a -(S)-set if any 4 E (2( satisfies the inequality ||pl| Σ A||L for some constant A(S) > O. The infimumm of all A > 0 for which is a A-set is denoted by A. Our theoretical results analyze how efficiently the SGW operator can represent. band-limited functions on a graph in terms of how many sampled nodes are nceded for stable reconstruction. The Ag-set associated with the SGW operator &b quantifies this sampling requirement: a smaller Ag means that all band-limited graph signals can be recovered from fewer samples while preserving stability, and is therefore desirable. An alternative way to view the efficiency of the SGW representation is that because the Ag-set associated with the SGW operator is smaller than the corresponding A-set for the Laplacian, one can stably represent band- limited functions with a higher cutoff frequency w, since the admissible bandwidth is inversely proportional to the A-set value. The ability of the SGW operator &a to efficiently represent functions f E (G) can be summarized in the following theorem (the proof is provided the Appendix)., which highlights the role of the Ag-set with respect to the operator vb. We show that any signal I E (G), where X S ag for some Ωg AN, can be uniquely represented by its samples the uniqueness set U using the SGW operator. Under certain conditions related to the SGW operator, its associated Ag-set is smaller than the A-set associated with the Laplacian operator.
Theorem 1 Let G = (V, W) be connected graph with n vertices and f E (W,(G) for A1 S w < Amaz- The SGW operalor to can be construcled such that the set S is a Ag-sel and the set U = VVS is a uniqueness set for any space (G) with w < 1/A and Ng < A for any 4p E Tr2(S), where A is the A-set of the Laplacian operator.
7 Discussion
Identifying the key factors in high-dimensional datasets is essential for leveraging unlabeled data in many practical applications. In this work, we introduce a novel contrastive learning framework for interpretable manifold learning via graph embeddings, leveraging both low- and high-frequency information to enhance representation quality. Our method employs Spectral Graph Wavelets (SGW) to construct a multi-scale graph representation of the input feature space. This representation is optimized using a stochastic gradient descent (SGD)-based scheme, combined with novel 3D tensor encoding strategy that enriches the embedding process.
MS-IMAP differs fundamentally from traditional dimensionality reduction or visualization techniques. While many manifold-learning methods (e.g., Laplacian Eigenmaps, Diffusion Maps, UMAP) aim to reveal intrin- sic structure through low-dimensional embeddings, MS-IMAP instead constructs structure-preserving inter- pretable embeddings for reanstream tasks by enforcing a one-to-one correspondence bet ween input features and embedding dimensions. The advantage of MS-IMAP does not come from dimensionality alone: simply increasing the number of embedding coordinates often does not lead to improved performance, as shown both in our comparison with graph scattering transforms and in the additional high-dimensional baselines reported in Appendix A.2.
To enhance the embedding with structural information, we introduce a 3D tensor that captures SGW coeffi- cients across nodes, features, and frequency bands. This enables the model to encode both global and local graph structure. Importantly, although SGW captures information across the graph, each embedding coor- dinate remains feature-aligned (one per input feature) while its update incorporates joint structural coupling through the shared graph and SGW operator. This preserves the interpretability and separability of features throughout training. We demonstrate the utility of these interpretable embeddings for both feature selection and traditional clustering tasks. Specifically, we assess clustering performance based on features seleeted from both the input space and the embedding space, using two unsupervised methods: Laplacian Score and Mutual Information- selection. A key advantage of our framework is its inherent linkage between the original and embedded fcature spaces, a property rarely offered by existing nonlinear manifold lcarning techniques. Furthermore, we show that the resulting embeddings are highly competitive with state-of thett-r graph embedding methods.
10

---

Published in Transactions on Machine Learning Research (04/2026)
To further justify the use of SGW for representation, we analyze its theoretical properties by studying functions in Paley-Wiener spaces on combinatorial graphs. We show that the SGW operator enables more expressive representations, characterized through the concept of the A-set.
Limitations and future directions include several aspects. First, our approach assumes that the input features -ad the similarity measures used- adde--ddeuallly capture the underlying manifold structure, particularly geodesic distances. Second, extending the current SGD-based optimization framework to handle out-of-sample generalization remains a challenge. Finally, while the bebbeddiig-featurre correspondence facil- itates global feature importance analysis, further improvements may be achieved by incorporating localized importance estimation techniques.
Acknowledgments
We thank the anonymous reviewers and the Action Editor for their helpful comments and suggestions, which improved the quality of this paper.
References
Ehsan Amid and Manfred K. Warmuth. Trimap: Large-scale dimensionality reduction using triplets. arXiv preprint Xrx:::1910. 00200,, 2019.
Davide Anguita, Alessandro Ghio, Luca Oneto, Xavier Parra, and Jorge L. Kleyes-Ortiz A public domain dataset for human activity recognition using smartphones. 21st European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning, pp. 437 442, 2013.
Mikhail Belkin and Partha Niyogi. Laplacian eigennaps for dimensionality reduction and data representation. Newral Computation, 2003.
R. R. Coifman, Lafon, A. B. Lee, M. Maggioni, F. Warner, and S. Zucker. Gcometric diffusions as tool for harmonic analysis and structure definition of data: Diffusion maps. In Proceedings of the National Academy of Sciences, 2005.
Sebastian Damrich and Fred A. Hamprecht. On umap's true loss function. In Neural Information Processing Systems, 2021.
Shay Deutsch and Gerard Medioni. Unsupervised learning using the tensor voting graph. In Scale Space and Variational Methods in Computer Vision (SSVM), 2015a.
Shay Deutsch and Gerard Medioni. Intersecting manifolds: Detection, segmentation, and labeling. In Proceedings of the Twenty- Fourth International Joint Conference on Artificial Intelligence, 2015b.
Shay Deutsch and Gérard G. Medioni. Learning the geometric structure of manifolds with singularities using the tensor voting graph. Journal of Mathematical Imaging and Vision, 47(3):402-422 2017.
Shay Deutsch, Antonio Ortega, and Gérard G. Medioni. Robust denoising of piece- piece-wise smooth manifolds. In 2018 IEEE International Conference Acoustics, Speech and Signal Processing, 1CASSP 2018, Calgary, AB, Canada, April 15-20, 2018, pp. 2786-2790, 2018.
Dheeru Dua and Casey Graff. Uci machine learning repository. https://archive. ics uci edu, 2017.
Christopher Frye, Damien de Mijolla, Tom Begley, Laurence Cowton, Megan Stanley, and Tlya Feige Shapley explainability on the data manifold. In International Conference on Learning Representations, 2021.
Fernando Gama, Joan Bruna, and Alejandro Ribeiro. Stability of graph scattering transforms. Proceedings of the 33rd International Conference Newral Information Processing Systems, 2019.
Feng Gao, Guy Wolf, and Matthew J. Hirn. Geometric scattering for graph data analysis. In International Conference on Machine Learning, 2018.
David K. Hammond, Pierre Vandergheynst, and Rémi Gribonval. Wavelets graphs via spectral graph theory. Applied and Computational Harmonic Analysis, 2011.
Xiaofei He, Deng Cai, and Partha Niyogi. Taplacian score for feature selection. In NIPS, pp. 507 514, 2005. Matthias Hein and Markus Maier. Manifold denoising. In Advances in Newral Information Processing Systems, 2006.
Babak Hosseini and Barbara Hammer. Interpretable discriminative dimensionality reduction and feature selection on the manifold. Springer International Publishing, 2020.
11

---

Published in Transactions on Machine Learning Research (04/2026)
Guillaume Huguet, Alexander Tong, Edward De Brouwer, Yanlei Zhang, Guy Wolf, Ian Adelstein, and Smita Krishnaswany. A heat diffusion perspective geodesic preserving dimensionality reduction. Tn shirty-seventh Conference on Neural Information Processing Systems, 2023.
Thomas N Kipf and Max Welling. Variational graph auto-encoders. arXiv preprint arXit:6611.0330 2016. Sadamori Kojaku, Jisung Yoon, Isabel Constantino, and Yong-Yeo Alm. Residual2vec: Debiasing graph embedding with random graphs. In NewrIPS, 2021.
Andrew Kachites McCallum, Kamal Nigam, Jason Rennie, and Kristie Seymore. Automating the construc- tion of internet portals with machine learning. Information Retrieval, 3:127-163, 2000.
Leland Mclunes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction, 2018. URI. wourraalofopessuuucresttttttoo tware.
Kevin R Moon, David Van Dijk, Zheng Wang, Scott Gigante, Daniel B Burkhardt, William S Chen, Kristina Yim, Antonia van den Elzen, Matthew J Hirn, Ronald R. Coifman, et al. Visualizing structure and transitions in high-dimensional biological data. Nature biotechnology, 37(12):1482-14999 2019.
Isaac Z. Pesenson. Sampling in saley-wiener spaces on combinatorial graphs. Transactions of the American Mathematical Society, pp. 5603 5627, 2008.
Mauro Maggioni Ronald R. Coifman. Diffusion wavelets. Applied and Computational Harmonic Analysis, 2006.
S. Roweis and L. Saul. Nonlinear dimensionality reduction by locally linear embedding, SCIENCE, 290: 2323 2326, 2000.
Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Fliassi-Rad. Collec- tive classification in network data. AI magazine, 2008.
Jian Sun, Maks Ovsjanikov, and Leonidas Guibas. A concise and provably informative multi-scale signature based on heat diffusion. Computer Gruphics Forum, 2009.
Atsushi Suzuki, Atsushi Nitanda, jing wang, Linchuan Xu, Kenji Yamanishi, and Marc Cavazza. Generaliza- tion bounds for graph embedding using negative sampling: Linear vs hyperbolic. In Advances in Newral Information Processing Systems, volume 34, pp. 1243 1255, 2021.
J. Tenenbaum, V. de Silva, and J. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 2000.
Laurens van der Maaten and Geoffrey Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research, 2008
Yingfan Wang, Haiyang Huang, Cynthia Rudin, and Yaron Shaposhnik. Understanding how dimension reduction tools work: An empirical approach to deciphering t-sne, umap, trimap, and pacmap for data visualization. Journal of Machine Learning Research, 2021.
Yongqin Xian, Bernt Schiele, and Zeynep Akata. Zero-shot learning- - the good, the bad and the ugly. 2017 IEEE Conference Computer Vision and Pattern Recognition (CVPR), pp. 3077-3086, 2017. Yuchen Yan, Baoyu Jing, Lihui Liu, Ruijie Wang, Jinning Li, Tarek Abdelzaher, and Hanghang Tong. Reconciling competing sampling strategies of network embedding. Tn NewrIPS, volume 36, 2023.
Arman S. Zharmagambetor and Miguel A. Carreira-Perpinaa Learning interpretable, tree-based projection mappings for nonlinear embeddings. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, pp. 9550-9570, 2022.
Rapolas Zilionis, Camilla Fngblom, Christina Pfirschke, Virginia Savova, David Zemmour, Hatice Duygu Saatcioglu, Indira Krishnan, Giorgia Maroni, Claire V Meyerovitz, C. Michael Kerwin, Sun Choi, William G Richards, Assunta De Rienzo, Daniel Geoffrey Tenen, Raphael Bueno, Elena Levantini, Mikael J. Pittet, and Allon M. Klein. Single-cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species. Immunity, 2019.
12

---

Published in Transactions on Machine Learning Research (04/2026)
Table Summary of main notation.
Symbol
Description
N
Number of nodes (samples)
D
Number of features (imput dimensions)
K
Number of SGW scales
E (l..... ......,] N
Node (sample) indices
de (1..... D
Feature indices
& {L.....
K}
Scale indices Input data matrix (nodes x features)
x Nx
L
Graph Laplacian
L UAUT
Taplacian eigendecomposition
U [uo .... UN-1
Matrix of Laplacian eigenvectors
A - diag( Ao 1N-1)
Diagonal matrix of eigenvalues
(s n)
SGW coefficient at scale Sk. feature d, node n
(Sk :;: E DxN
SGW tensor slice at scale 8% (all features, all nodes) SGW coefficients of feature d across nodes at scale 8g SGW coefficients of all features node at scale 8 Optimized embodding at scale 8
xx (8k; d, R.
xx (Sk: ., E R
xx (sk: :,:)
1x D
Final MS-IMAP embedding (aggregated over scales)
Wrm
Edge weight between nodes and m
6olx (8k)
SGW-basod structural similarity between and at scale 8k
A
Appendix
A.1 Theoretical results: Sampling set for Smooth Manifolds with functions defined over Paley-Wiener Spaces
In this section, characterize the thcoretical properties of the representation power of the SGW operator by considering functions sampled from the Faley-Wience spaces (Pesenson (2008)) combinatorial graphs. The Paley-Wis Wiener spaces were introde combinatorial graphs in (Pesenson (2008)) and corresponding sampling theory was developed which resembles the classical one. Pesenson proved in (Pesenson (2008)) that Paley-Wiener functions of low type are uniquely determined by their values on certain subgraphs (which are composed from a set of nodes known as the uniqueness sets) and can be reconstructed from such sets in a stable way. We demonstrate that the SGW operator can represent functions f that reside the Palcy-Wience space on combinatorial graphs more efficiently than the graph Laplacian operator L. The effectiveness of the SGW operator representation in this case can be understood in several ways. In way, by the ability of the SGW operator to accurately represent functions with larger bandwidth, i.e., f E PG) where w < w' We first summarize the main notions and definitions. The space L2(G) is the Hilbert space of all complex-valued functions f V + C with the following inner product (S.g) = EveG (V)g(() and the norm
|I (
1/2
Σ
(8)
The Laplace (normalized) operator L defined by the formula (Pesenson (2008)):
f(v) Vd(u)
f(u) Vd(u)
Cf(e) / dd(v)
L2(G) L2(G)
E
(9)
ane
In order to prove these results, we will first need the following definitions:
Definition 1 The Paley- Wiener space of -bandlmmtted signals f E Iz(lG) is defined follows:
13

---

Published in Transactions on Machine Learning Research (04/2026)
- = (1((0)
(10)
We consider a simple, undirected, unweighted, and connected graph (V. W), where is its set of N vertices and W is its set of edges. The degree of is number of vertices adjacent a vertex is and is denoted by d(v). We assume that degrees of all vertices are bounded by the maximum degree denoted
d(G) = max.evd(r).
(11)
The following definition (Pesenson (2008)) explains the uniqueness set.
Definition 2 A set of vertices U C V is a uniqueness set for space (G) if for every two functions from (G) that coincide on U, then they coincide on V.
Definition 3 For subset S C V, denote L2(S) as the space of all functions from L2(G) with support in S:
L2(5)= = fop ((()) = 0,000 -
Definition 4 Pesenson (2008) We say that set of vertices S c V is a A -set if for any E L2(S) it admits a Poincore inequality with a constant A = A(S)= > 0
||ll| S Al|Coll, & E (S).
(12)
The infimuum of all A(S) > for which S is A -set will be called the Poincore constant of the set S and denoted by A.
The definition above provides tool to determine when bandlimited signals in Paley- Wileer spaces PW.(G) (G) can be uniquely represented from their samples given set. The role of A-sets the following Theorem by Pesenson (Pesenson (2008)), that shows that if S C V then any signal f PW) can be uniquely represented by its function values in the complement set U = V(G) \ S:
explained and proved in
Theorem 2 Pesenson (2008) If c V is a A- set, then the set U = V(G) \ is uniqueness set for any space (G) with w < x X-
Remark: Note that non-trivial uniqueness sets can not exist for functions from any Paley-Wiener Wiener subspace (G) with any A S w S AN, but they can exist for some range An S w < Q. as was shown in (Pesenson (2008)).
We state one of our main results, in which we employ the SGW operator v to characterize the uniqueness set using the set therefore extending the A-set concerning the graph Laplacian operator L.
Theorem 3 Let G = (V, W) be a connected graph with N vertices. Assume that there exist a set of vertices S V for which the conditions (1)-(2) in Lemma 1 below hold true. Let tp be the SGW operator using a polynomial p(C)) with the coefficients {ax K such that tbs = o f. Then, for any ip I2(S), we have that the following inequality holds:
||aا|| < A l1b2ll
(13)
and thus the set S is Ay-set for the operator with Ag
Σ
D A2k
We recall the following results from (Pesenson (2008)). Note that (Pesenson (2008)) established the con- struction of A-set by imposing specific assumptions the sets S and U. Our result Theorem holds similar assumptions.
Lemma 1 (Lemma 3.6 of Pesenson (2008)) Given a connected graph G (V, W), set of vertices S c V, is complement U V V.S, for which the following conditions hold true.
14

---

Published in Transactions on Machine Learning Research (04/2026)
1. For every S there crists U that is neighbor of i.e., w(u, 8) 0.
2. Por every E S there exists at least one neighbor node u E U whose adjacency set intersecls S only over 8.
Then there erist set of vertices S C V which is A-set, with A = d(G).
We require one more property the Laplacian operation, Lemma below. To this end, we use a lemma from (Pesenson (2008));
Lemma 2 (Lemma 3.9 of Pesenson (2008)) If S is A-sel, then for any E Lz(S) and all 2 0, k == 2', 1= 0,1,
A pll s | pll, op L2(S)
Now we are ready to prove modification of Lemma 3.9 from (Pesenson (2008)):
Lemma 3 Let S be A-sel, v L2(S), and let m be on inleger. Then
|| [["*"ll..
Am
Because every positive integer has binary representation, then we have for some js E 20..,2....3, 1,2, N, that
Proof:
24.
m
Using Lemma 2, let h: = 2 and t = 2
2
= k + t. Then we have
|| = +|
2 A ||
|| 4..44...00 2
2
bll
2
We can now let h2 2 and to 20'
+ 2 that,
|| 2
[
2
1
2
1
2 2 A
||
2
s0|l-
2-
Continuing in this way, eventually get,
[ pl 2
Am ||
completing the proof.
In the next theorem, we expand the characterization of the uniqueness set using the A-set concerning the graph Laplacian operator L to include cases where we employ the SGW operator z.. We thereby characterize using the Ag-set for the SGW operator. We now turn to prove Theorem 3, which was stated earlier. Proof of Theorem 3:
the uniqueness
15

---

Published in Transactions on Machine Learning Research (04/2026)
Let S C V be set satisfying the conditions in Lemma 1, S is A-set. Then by Lemma 3, for any integer k 2 !
A |l s p|l
The SGW operator ip has the form,
M/
0g
5 k=0
for some 2 1, ag I are coefficients, and where the summation over i and represents the summation the nodes of the graph. Now we calculate,
1/2
Σ i.j
Σ k=0
ax (C s)
=
1/2
Σ
[<* ?
1/2
Σ
(()s)))s)
0
K
1/2
Σ
<+1l=+>p1?
K
1/2
1
2 |14p11
0
1=
K
1/2
= ||
1:
2
k=0
Letting
Ap
Ek=
A
Then we have ||p||| s Aollvoll so S is also a setset
From the assumptions and proof of Theorem 3, we get the immediate result:
Corollary 1 Let c V be a A-set. Then is also a Ageset.
Remark 1: An important property which can be observed from Theorem 3 is the following: given the SGW y as a spectral representation operator applied on f. we can choose coefficients {au]k=o such that we obtain Ap set associated with the operator v., which is smaller than the A set associated with the Laplacian operator L. This implics that the operator & provides more flexibility and better control over smoothness properties in comparison to the Laplacian operator.
Remark 2: Since the Ag set can be chosen be smaller than A
(for proper choice of the coelficients { using the operator v/) then the SGW operator & provides a more efficient representation for f e PW (G) with w < w' using the same subsets of nodes from the uniqueness set U in comparison to the Laplacian operator L.
Remark 3: Note that the characterization of the the uniqueness set does not rely on reconstruction
method of the graph signal values of f(S) from their known values U.
16

---

Published in Transactions on Machine Learning Research (04/2026)
The next Theorem demonstrates the role of the Ag -set with respect to the operator e where show that any signal I E (G) can be uniquely represented by its samples on the uniqueness set U. This results resembles the role of A-sets with respect to the graph Laplacian operator L, yet with different bound then Lemma
Theorem 4 Let G (V, W) be connected graph with N verlices and E PW)(G) for A <w < Xmar- Given the SGW operator 1,, and a set S which is Ag set. Then the set U = V is uniqueness set for any space (G) with w < 1/Ap.
Given f.f E (G), then f - f (G). Assume that f * If f.f coincide U = VVS. then f f E (2(S) and therefore
Proof:
|lf - =
(14)
Since by E R., we have that by properties of a vector space in RN, using the Cauchy- Schwarz inequality and assuming layl < 1 Vk, we have:
, S wllf fll-
(15)
Combining the inequalities above and using the inequality Ayw < 1 we have that
||f - fll S A S Agwlls fll < IIf - fll
(16)
which is a contradiction to the assumption that f * Thus, the set U = V S is a uniqueness set for any space (G) with w < 1/Ag-
Remark 1: Note that Aw < 1 implies that Agu 1 given A < A, then
increase the size number of nodes in the uniqueness set U for (G) (for a proper choice of the coefficients {ax using the operator v). In other words, we may increase the size of (thus reducing the size of U) and still obtain a uniqueness set with a smaller size for the graph signals in PG)
Remark 2: We note that the results of Theorem 4 concerning the uniqueness are independent from the stability properties of the representation. In order to achieve stability which is important for reconstruction. it is required to construct a wavelet operator using multiple scales fj.j = 1,...,, as proposed in (Hammond et al. (2011)). More specifically, we can express the function &o using multiple scales lj (here we replace the previous notation of scale with L = 1,..J, not to confuse with nodes 8 E S). Then for a fixed scale ij 1/2
we have that the SGW is given by [1)b,(1))11 = (E Ek l (C*)p(())))
and
1/2
N K ΣΣΣ
1k
Σ
a (C (()))
(17)
In similar way to the arguments provided in Thcorem 3 choose coefficients aty, associated with the Laplacian polynomial such that the inequality ||s]|| S Appl is satisfied.
Remark 3: Optimal Choice of coefficients for Minimizing Ag for optimal choice of Ap in the multi-scale setting.
To find the coefficients a that minimize Ay where:
Ap
(18)
V V
I
A
i we optimize the coefficients at, which control how much of the signal at scale 1 and power k is captured by the SGW operator. We the Lagrangian F.
subject to the constraint: Σ/=1 Ek o
17

---

Published in Transactions on Machine Learning Research (04/2026)
incorporating the energy constraint into the optimization problem and Define:
J
K
+ P
Σ
oi kmma
-1
1
(19)
Ex k=
wherep is the Lagrange multiplier that ensures the total energy of the coefficients is equal to 1. Taking partial derivalives, solving for 0r and differentiating L with respect to each coefficient dty.
OF
at, Rt..k// 2k
+ 2pda, & = O.
(20)
Dat,k
o? f Ex
3/2
Eラ
k=
n
rearranging the equation above we obtain:
3/2
J K
at
Σ
2
A 2
A
2
j=l
R==
assuming 0t,,k # 0, we can divide both sides of the equation by aty. yiclding:
3/2
0 *
2p
ΣΣ
(21)
2
2
Z is proportional to the inverse of A implying that the optimal wavelet thus the sum Σ/=1 Ek A coefficients should decay proportionally to the inverse powers the Laplacian eigenvalues at,,k x Ne.
Next, applying the energy constraint: ΣE=1 Σk o and substituting 0t AEr obtain:
Σ
(22)
A
j=1 k=0
Thus, the normalized optimal coefficients at are:
0t
(23)
Σ Σ A2F
Note that the wavelet energy is distributed across all scales and powers of the Laplacian, with a decay that prioritizes low-frequency components (small k). Using the form of the coefficients ar above, the expression for A
A
(24)
V MK k= A
which is the minimized value of Ag, showing that the optimal choice of coefficients reduces the overall constant, improving the localization and efficiency of the SGW operator. Thus, the optimal choice wavelet coefficients a that minimizes Ay in the multi-scale SGW operator is proportional to 1 low-frequency components.
x prioritizing
18

---

Published in Transactions on Machine Learning Research (04/2026)
A.2 Additional analysis of embedding dimensionality
This section provides additional clarification and experiments regarding the role of embedding dlimensionality and the comparison between MS-IMAP and high-capacity baselines.
Dimensionality comparison. As discussed in Sec. 5.2, MS-IMAP is not a dimensionality reduction or visualization method. It is an interpretable graph embedding that preserves a one-to-one correspondence between input features and embedding coordinates by design, and therefore retains the original feature di- mensionality rather than compressing 2-3 dimensions. The comparison in Sec. 5.2 is intended to situate MS-IMAP with respect to standard manifold learning and visualization baselines and to show that, even under the additional constraint of preserving one-to-one feature- cabbedi-n correspondence, its perfor- mance remains comparable to these methods. The goal is not to altribule performance solely to higher dimensionality; in fact, increasing dimensionality alone often does not guarantee improved performance, as illustrated both in our comparison with graph scattering transforms and in the additional results reported here.
To more directly assess the effect of representational capacity, we include comparisons to higher-capacity un- supervised representation methods (e.g,, autoencoder self-supervised baselines) and to higher-dimensional variants of manifold learning methods using multiple choices of embedding dimensionality. This allows to test whether MS-IMAP's benefits persist when competing methods are not constrained to 2-3 dimensions and to make the role of dimensionality more transparent.
We specifically examine:
Unsupervised deep representations (GAE). We add a Graph Auto-Encoder (GAE) baseline following (Kipf & Welling (2016)), implemented via PyTorch Geometric's GAE class, with 2-layer GCN encoder (Kipf and Welling, 2017) and an inner-product decoder trained in self-supervised fashion to reconstruct the INN adjacency matrix built from the input features. We tune the em- bedding dimensionality over
d {2, 5, 10, 20,50, D),
where D is the input dimensionality, and report in Table the best clustering performance (ARI/AMI) achieved by GAE over these values of d for each dataset.
2. Higher-dimensional manifold learning basclines. We likewise re-run the manifold learning methods UMAP, HeatGeo, and PHATE with the same set of oulput dimensions,
d (2, 5, 20, 50, D),
and again report in Table the best performance over d for each method and dataset.
The experimental results show that, in most cases, increasing the embedding dimension beyond the standard 2-3 dimensions does not yield substantial performance gains for these methods. Across datasets, MS-IMAP remains competitive with state-of-the-att baselines, typically achieving the first- or second-best ARI/AMI For GAE in particular, even with careful tuning of d, its clustering performance on these tabular &NN-graph datasets is comparable to or below that of simpler manifold-learning baselines. This is consistent with prior observations that unsupervised deep graph encoders can struggle to produce effective embeddings for clustering when the underlying graphs are irregular or constructed from tabular features rather than arising from natural network structure.
Remark: Please note that the discrepancy between Table and Table 3 comes from differences in the evaluation protocol rather than from MS-IMAP itsclf. In Table 1, cach bascling (including UMAP and PHATE) was run with its own hyperparameter sweep tuned for the standard low-dimensional use case (e.g.. recommended configurations at 2-3 dimensions), and we reported the best result from that sweep. In Table 3, the goal specifically to study the effect of embedding dimensionality. To keep this experiment tractable, we fixed the remaining hyperparameters across the different values of (d 2,5, 10, 20, 50, D) (and, for some methods, averaged over multiple random seeds) rather than performing full hyperparameter search for
19

---

Published in Transactions on Machine Learning Research (04/2026)
Table 3: Clustering performance (ARI / AMI) for the dimensionality experiment. For each method and
dataset we report the best score over d {2, 5, 10, 20, 50, D}.
HAR
Census
Zillionis
AWA
Cora
Method
ARI
AMI
ART
AMI
ARI
AMT
ART
AMI
ART
AMT
UMAP
0.64
0.76
0.22
0.14
0.61
0.76
0.8
0.82
0.35
0.38
HeatGeo
0.63
0.71
N/A
N/A
N/A
N/A
0.77
0.82
0.18
0.22
PHATE
0.56
0.67
0.2
0.092
0.27
0.71
0.76
0.79
0.16
N/A
GAE
0.39
0.32
0.05
0.09
0.53
0.64
0.7
0.77
0.252
0.30
MS-IMAP
0.67
0.80
0.23
0.15
0.70
0.76
0.74
0.81
0.33
0.36
(3.a) HAR
(3.b) Zillion
(3.c) Cora
Figure 3: Silhouette score as function of the number of selected features k on HAR, Zillion, and Cora, comparing the label-free clustering quality obtained from the raw features (blue) and from the MS-IMAP embedding (orange), both ranked by Laplacian Score. For each dataset and cach k, we run k-mcans with several choices of the number of clusters C and report the best silhouette score over C.
every method, d and dataset combination. As consequence, the exact configuration that achieved the best scores in Table is not guaranteed to appear among the Table 3 runs, so the "best over (d)" values in Table can be slightly lower than the Table numbers for some methods and datasets. Thus Table 3 show that under a consistent and reasonably tuned setup across dimensions, increasing the embedding dimension beyond 3 does not systematically improve the performance of these baselines, and that MS-IMAP remains competitive this higher-capacity regime.
Finally, the choice to retain the input dimensionalitty in MS-IMAP is primarily motivated by interpretability: most nonlinear embeddings sacrifice the direct, feature-wise correspondence that is central to the feature- importance analyses in Secs. 4-5.1. The results in Table 3 indicate that this interpretability does not come at the cost of substantial performance loss.
A.3 Silhouette-basee analysis of feature selection
To further analyze the effect of feature selection on clustering structure, we report the internal clustering metric of the silhouette score on the selected feature subsets. This label-free mcasure does not use ground- truth labels and therefore does not replace our primary ART/AMT evaluation, but provides an additional unsupervised check that the selected subsets induce coherent clusters the embedding space. For cach dataset and each number of selected features k, we run k-means on the selected features for several values of the cluster number C' (e.g., C' [4.9) for HAR, C e [7, 13] for AWA, C' E [2,4] for Census) and record the best silhouette score over C.
Figures 3 4 show the resulting curves when Laplacian Score is applied either to the original features or to the MS-IMAP embedding. For HAR, Zillion, AWA, and Census, the MS-IMAP embedding consistently achieves higher and more stable silhouctte scores (typically around 0.6-0.7) than the original fcatures, whose scores degrade steadily as increases. For Cora, the original feature space exhibits a very high silhouette for the top one two features but then rapidly collapses toward values near zero
features added,
20

---

Published in Transactions on Machine Learning Research (04/2026)
Test Sinccette
Fest ditcueete us
Kii
Bost baelle wer CEO
Streeteeccm 60
Dcd altccdiicccc oud K Ensoocinng
(4.a) AWA
(4.b) Census
Figure 4: Silhouette score as function of the number of selected features k on AWA and Census, comparing the label-free clustering quality obtained from the raw features (bluc) and from the MS-IMAP embedding (orange), both ranked by Laplacian Score. As Fig. 3, for each dataset and each k we report the best silhouette score over several choices of the number of clusters C.
whereas the MS-IMAP embedding maintains moderate but stable silhouette (around 0.3) over wide range of k. Overall, MS-IMAP based feature selection trades narrow silhouette peak for small number of raw features for more robust and consistently higher cluster quality once k exceeds few features. A systematic study of model-selection strategies for choosing C' or the embedding dimension (e.g., eigengap silhouette-based selection) in combination with MS-IMAP is left for future work.
A.4 Alternative Encoding: Encoding By Simple concatenation
An alternative to the 3D tensor-basod representation approach presented in the main body of the paper, one can simple concatenation, where all features are transformed using SGW and each data point is represented by vector formed by concatenating all corresponding SGW coefficients. For completeness, we also provide simple encoding method that involves concatenating the multi-scale representation of all features and filters (associated with different scales) into a single vector representation for each point. This results in a matrix representation denoted as vix We detail the encoding and optimization employed for this method. Note that we designate the concatenation using ||. with
c(wy, v., (sk,:) = why (sk, (8k,,)))))))
denoting the concatenation of the vectors corresponding to the multi-scale representation v (sk,:) and 1 (s).:). For method all features and all scales are concatenated together, the resulting matrix thx can be represented as
thx (x
(sKs,))
where tix E (XKD))
1. Optimize Embedding Method 1:
Given the encoded multi-scale representation vix E (X(D))
Ux E RKD N perform optimization in the SGW domain, using the following fuzzy cross entropy loss function:
C W) Σ
Wnm
Wam V'e
+ (1 Wmm blgg
1 Wnnm
(25)
1xx
Uno
Unm
where Cpm
1+ 1+10xxn Vxm.
21

---

Published in Transactions on Machine Learning Research (04/2026)
Feature space
Initialization Manifold latent
Regularized Manifold space
(f1)
422
V2
4/ (nwi)
f
Figure 5: The encoding step of our framework illustrated as mapping between each coordinate of the input features and corresponding dimension in the embedding space. This mapping is facilitated through our proposed approach, which aligns the input dimensions of the features with the latent embedding space. By establishing this correspondence, we acquire the capability to interpret each dimension in the embedding space. Such interpretation allows for various analyses, including providing feature importance relative to the latent space.
AMI vs Number Features Kept
AMI vs Number of Features Kept
0.35
Adjusted Mutual Information (AMI)
0.7
0.30
Information
0.5
Original CLSI
0.5
Original ILSI
MS-4NAP Embeiding (1:)
MS-NMAP Embeiding 8.5)
Original OMD
Original (MI1
(MS-IMAP Embeoding (MI
M1--NAA Embeiding (o)
adiusted
0.3
0.13
0.05
0.2
0.00
Number Features Kept
Number of Features
(6.a) Cora dataset (AMT)
(6.b) Zillions dataset (AMT)
Figure 6: Clustering results Cora and Zillions datasets using subset of features based feature importance, comparing the MS-IMAP embedding space and the input features. Plots (a)-(b) show Adjusted Mutual Information (AMT) results.
Dropping terms that do not depends on t's,, the gradient of the loss is approximated by:
0c W)
Wrm v. nm
(ram)
Σ m
the Prm (vxx., Vixm )
(26)
|[am 12
m
where Fnm r/x. - xxm. and v'x. corresponds to the vector associated with node n multi-scale representation concatenated for all the graph signals and all spectral bands.
22

---

Published in Transactions on Machine Learning Research (04/2026)
Algorithm 2: MI-based Unsupervised Feature Selection
Input: MS-IMAP embeddings dxx, number of clusters kt
Output: MI scores for each feature
Step 1: Apply k-means clustering on the embedding space vhe to oblain cluster labels
clusters = k-means((yy) (1x, k)
Step 2: Use the resulting cluster labels as pseudo-labels
pseudo_labee labels c(txx)
Step 3: Initialize an empty list to store mutual information (MT) values for each feature
mi pcores=|]
Step 4: For each feature d in the original data or embedding space:
Step 4.1: Compute mutual information between feature (d) and the pseudo- labels
Step 4.2: MI score = (I(xxxd) o(t)x)
Step 4.3: Append MI score to mi scores
Step 5: Output the list of MI scores
Return mi scores
A.5
Additional details: Embedding Feature Importance using Mutual Information and Laplacian Score
A.5.1
Computing the Laplacian Score with respect to the embedding features
The Laplacian score is calculated with respect to embedding features dix using the Laplacian graph L and the degree matrix D. To compute the Laplacian Score, we first subtract the mean and then calculate it as follows: I, (vx)aa (w) L(WW) Smaller scores indicate that the embedding feature (Vx)) has a greater projection D
onlo the subspace of eigenvectors corresponding to the smallest eigenvalues, signifying higher importance concerning the global graph structure. Therefore, the feature importance of (vx)a can be directly interpreted as the importance of the corresponding original feature. Although the graph is typically constructed using local neighborhoods (e.g., via k-nearest neighbors), the Laplacian operator reflects global smoothness over the entire data manifold.
A.5.2 Embedding Feature Importance using Mutual Information
We summarize the MI-based unsupervised feature selection procedure in Algorithm 2., as described in Scction MI based unsupervised feature selection self-supervised (using the psoudo-labels), capturing which feat are useful for distinguishing between clusters. MI scores are also directly interpretable, where high MI score indicates that the feature contributes to separating groups.
In comparison, applying the Laplacian Score(He et al. (2005)) to the embedding may be less effective for feature importance, as the embedding process can mix original features in ways that obscure their individual contributions. While the Laplacian Score prioritizes local neighborhood preservation, this does not necessar- ily correspond to discriminative power for clustering. In contrast, the proposed MI-based approach directly identifies features that are informative for distinguishing clusters. However, relies on the quality of the pseudo-labels derived from clustering, which may not always reflect true cluster structure.
We also include additional experiments that demonstrate the effectiveness of MS-IMAP for feature impor- tance by evaluating clustering performance function of the number of selected features features-complenee complementing the results presented in Section 4. Specifically, we assess the most significant features identified using Laplacian Score and Mutual Information-bast selection, applied to both the MS-IMAP embedding and the original feature space.
Clustering results, evaluated using Adjusted Rand Index (ART) and Adjusted Mutual Information (AMT), are obtained by applying k-means different subsets of embedding features selected by each method, and compared to the corresponding subsets selected directly from the raw features.
23

---

Published in Transactions on Machine Learning Research (04/2026)
(7.a) HAR dataset (ARI)
(7.b) HAR dataset (AMI)
Figure 7: Clustering results (AMI) on the HAR datasets using a subset of features based on feature impor- tance, comparing MS-IMAP embedding space and the input features.
(8.a) AWA datasel (ARI)
(8.b) AWA dataset (AMI)
Figure 8: Clustering results ARI and AMI on the AWA dataset using a subset of features based on feature importance, comparing MS-IMAP embedding space and the input features.
On the HAR and AWA datasets, Figures 7, 8, show that, for both feature selection methods, clustering accuracy consistently higher when using the MS-IMAP embedding than when using the original fea- tures-acros range of feature subset sizes and for both evaluation motries. This improvement highlights the robustness of the MS-IMAP embedding, while maintaining interpretability through its one-to-one corre- spondence between embedding dimensions and original features.
A.6
Dataset Details
The two moons dataset depicts two interleaving half-circles. We sampled N 600 points and used Gaussian noise level having standard deviation 0.12. An example is show in Figure 9. To produce these points we use Sci-kit Learn's sklearn. sklearn.dataseen wwtasets........... make socee__mons.._o In Figure 9, we show example clustering for each of the methods mentioned in the Experimental Results section of the main paper.
Two Moons.
Census. From the UCI Machine Learning Repository Dua & Graff (2017), this dataset contains 14 features that are mix of categorical, numerical, and binary. Such features include ago, marital status, sex, etc. The goal is to predict whether sample makes less than equal to $50, 000, strictly more. We use 32, 561 samples in our dataset.
Lung Cancer. The Zilionis dataset is widely used, and consists of single-cell RNA sequencing data. It has 306 features, and 48, 969 samples. The data has 20 classes corresponding to cell type. More can be found at Zilionis et al. (2019) Animals with Attributes. The Animals with Attributes (AWA) dataset, contains 5,000 data points corresponding 10 unseen classes, where the testing image features are obtained from the pre-trained ResNet architecture, with D 2, 048 dimensions, and the semantic features are provided with D 85 dimensions. More information be found in Section 4.1 of Xian al. (2017).
24

---

Published in Transactions on Machine Learning Research (04/2026)
14-
10
147
105
20-
0.
-15
-05
32
-03
20
(a) Two moons
(b) UMAP
(c) t-SNE
(d) HeatGeo clustering
(e) Proposed Clustering
input
Clustering
Clustering
Figure 9: Clustering results using two moons
Cora. The Cora dataset has 2,708 scientific publications categorized into seven classes. The network has 5,429 links and each publication is described as a binary word vector indicating the presence or absence of a word. The dictionary consists of 1,433 words. In our experiments, we do not use the given graph, and only use the features. Sen et al. (2008); McCallum et al. (2000).
HAR. The Human Activity Recognition (HAR) dataset consists of sensor data collected from 30 subjects performing six different actívitios: walking, walking upstairs, walking downstairs, sitting, standing, and lying. The data was recorded using accelerometer and gyroscope sensors embedded in a smartphone worn on the waist. Each activity is represented as a feature vector of 561 attributes and 7352 samples, extracted from the raw sensor signals through signal processing techniques. In our experiments, we do not use the temporal dependencics and only utilize the extracted fcatures.
A.7 Experimental Results on two moons Synthetic Dataset
Here show further experiments more data, where assess the robustness of our method using the two moons dataset, which is 2D manifold depicting two interleaving half-circles. We sampled N 600 points, and set the Gaussian noise standard deviation to 0.12. While spectral based methods such as UMAP are effective under relatively "modest" noise levels, their performance deteriorates in the presence larger amounts of noise. As shown in Table I under the Moons column, our approach robust and correctly clusters most points despite the large noise level, and competitive with the competing methods.
Data
Two Moons dataset
Method/ //ccur/y
ARI
AMI
Features
0.24
0.18
UMAP
0.54
0.51
t-SNE
0.42
0.35
ISOMAP
0.36
0.3
Scattering Geometric.1 T
0.26
0.19
Diffusion Maps
0.25
0.19
PHATE
0.42
0.35
HcatGoo
0.54
0.52
MS-IMAP
0.89
0.87
Table 4: Comparison of clustering performance on the two moons dataset.
A.8 Hyperparameter Tuning Details
We note that for UMAP, t-SNE, PHATE, PaCMAP, TriMap and Heat we hyperparameter tune the performance on each dataset, rumning the method five times for each hyperparameter configuration and taking the average. Then we report the best average score. Additional details regarding hyperparameter
25

---

Published in Transactions on Machine Learning Research (04/2026)
are provided below.
We tune UMAP on each dataset in the Experimental Results section, by tuming over the parameter space in Table
Hyperparameter
Set of values
neighbors
2, 10, 15, 20, 30, 50, 100 0, 0.1, 0.5, 0.99
min dist
components
2, 3
Table 5: Space of parameters in which we tuned UMAP.
A.9 TSNE Hyperparameter Tuning Details
We tune TSNE on each dataset in the Experimental Rosults section, by tuning over the paramcter space in Table 6.
Hyperparameter
Set of values
perplexity
15, 30, 60, 200
exaggeration
12, 24
components
2, 3
Table 6: Space parameters in which we tuned TSNE.
A.10 PHATE Hyperparameter Tuning Details
We tune PHATE on each dataset the Experimental Results section, by tuning over the parameter space
in Table
Hyperparameter
Set of values
knn
3, 5, 10
decay (i.e. a)
6, 10, 14, 18, 40 auto, 50, 100, 150, 200
t
Table 7: Space of parameters which we tuned PHATE.
A.11
HeatGeo Hyperparameter Tuning Details
We tune HeatGeo on each dataset in the Experimental Results section, by tuning over the parameter space in Table
A.12 Ablation Studies of Hyperparameters for MS-IMAP
Here we study what effect the hyperparameters the number of nearest neighbors, and the number of filters have on the clustering performance of MS-IMAP. In table 9, we demonstrate the effect of varying the number of neighbors, between 10, 15, and 20 neighbors. We see that the results are mostly the same, thus showing MS- 1MMPP is robust on real datasets.
We also study the effect of varying the number filters. In table 10, we also see similar results when using different filters, showing the stability of MS-IMAP.
26

---

Published in Transactions on Machine Learning Research (04/2026)
Hyperparameter
Set of values
knn
5, 10, 15
lap type
normalized, combinatorial
harnack regul
0, 0.5, 1
Table 8: Space of parameters in which we tuned HeatGeo.
Dataset
Census
Zilionis
AwA
Number of Neighbors
Accuracy
ART
AMI
ARI
AMI
ARI
AMT
10
0.22
0.15
0.71
0.78
0.73
0.80
0.23
0.15
0.70
0.76
0.74
0.81
20
0.22
0.15
0.70
0.77
0.71
0.79
Table 9: Ablation study MS-IMAP with varying the number of neighbors. The number of filters is kept.
at 5.
A.13 Ablation Studies of Hyperparameters for t-SNE, Isomaps
We study the effect of varying the number of neighbors for the methods: t-SNE, Isomaps, and Diffusion Maps. For t-SNE this would be the perplexity hyperparameter, for Isomaps this would be the mumber of neighbors, and for diffusion maps, this would be the parameter that affects the width the Gaussian kernel, i.c. exp(-/a).
In Table 11, we choosing smaller perplexity of 15 does worse than the perplexity of 30, 60 and 200.
In Table 12, see similar pattern to t-SNE, where choosing too low causes reduction in performance. But performance stabilizes around choosing the number of neighbors as 5-10+.
Dataset
Census
Zilionis
AwA
Number of Filters/
Accuracy
AR1
AMI
ARI
AMI
ARI
AMI
5
0.23
0.15
0.70
0.76
0.74
0.81
6
0.22
0.15
0.70
0.76
0.73
0.81
7
0.22
0.15
0.72
0.78
0.74
0.81
Table 10: Ablation study MS-IMAP with encoding method 2, varying the number of filters. The number of neighbors is kept at 15.
27

---

Published in Transactions on Machine Learning Research (04/2026)
Dataset
Census
Zilionis
AwA
Perplexity
Accuracy
ART
AMI
ARI
AMI
ARI
AMI
15
0.03
0.04
0.37
0.67
0.69
0.76
30
0.15
0.15
0.38
0.68
0.73
0.80
60
0.17
0.17
0.38
0.69
0.73
0.80
200
0.2
0.24
0.39
0.69
0.71
0.80
Table 11: Ablation study on t-SNE, varying the perplexity.
Datasct
Census
Zilionis
AwA
Neighbors
Accuracy
ARI
AMI
ARI
AMI
ARI
AMI
2
-0.05
0.02
0.41
0.55
0.43
0.58
5
0.18
0.09
0.44
0.56
0.52
0.60
10
0.18
0.09
0.44
0.57
0.51
0.62
Table 12: Ablation study on Isomap, varying the number of neighbors.
A.14 Fast computation using Chebyshev polynomials
We provide additional details regarding the fast computation of SGW coelficients Hammond et al. (2011).
Directly computing the SGW coefficients above requires calculating the entire eigensystem of the Laplacian, which is computationally intensive- O(N3) for N points. Instead, Hammond ct al. Hammond et al. (2011) suggested computing the SGW using a fast algorithm based on approximating the scaled generating kernels through low-order polynomials. The wavelet coefficients at each scale are then computed as a polynomial of L applied to the input data, using approximating polynomials given by truncated Chebyshev polynomials.
The Chebyshev polynomials (y) are computed using the recursive relations: Te(y) 2yTTx Tx-2(y) for & 2 2, where To 1 and Ti y-
The SGW coefficients are then approximated using wavelet and scaling function coefficients as follows:
K (5((5,)) ~ (Sonss Σ EenTa(()))
(27)
where cj,bj > O are the Chebyshev coefficients and T3, are the shifted Chebyshev polynomials (x(x) = Tx (==a) for x E [0, Amax], where x = aly + 1), a Amax- The scaling function coefficients, which are corresponding to a low-pass filter operation, are approximated in a similar way using Chebyshev polynomials. Nole that the scaling kernel function a low pass filter h satisfying h(0) > and h(x)) - 0 when 00. the graph is sparse, we obtain a fast computation of the matrix-vector multiplication TyA (L),, where the computational complexity scales linearly with the number of points, resulting in complexity of O(N) for inpul signal f E N The SGWs elficiently map an inpul graph signal (a vector of dimension N) to NK K scaling and wavelet coefficients.
A.15 Runtime and Computational Complexity
We have performed experimental runtime studies on empirical datascts. The execution time of our method with Python code implementation with 32 cores Intel Xeon 8259CL rumning at 2.50Ghz and 256GB of RAM on the Cancer QC data-set of 48,969 samples and 306 features took 27 mins, on Cora 5 min, on the AWA dataset 1.35 mins and on Census 7 min. In ficlds such as finance and genomics/bioi sizes range from under 100,000 to several million data points, analysis is often conducted offinee once the graph embedding computed. Runtimes of up to few hours generally acceptable and do not present
where typical datasct
28

---

Published in Transactions on Machine Learning Research (04/2026)
significant challenge, making the method well-suited for real-world applications. The computational com- plexity of MS-IMAP is of O(ND log( N)) for construction of the multi-scale represental ions which includes the k nearest neighbor graph using k-d tree, the SGW transform which is O(N) for each dimension of the manifold for sparse graphs. The optimization stage has a complexity which scale with the number of edges in the graph, which has complexity of O(kDN). Note that our approach generally requires the embedding dimensionality to be least equal to the input dimensionality of the features, which adds computational complexity.
29


Reviews:
[{"review_id": "R9xiJYgG52", "reviewer": "Reviewer_rd3L", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces the method MS-IMAP, a new representation learning method. It uses spectral wavelets at multiple scales and learns an embedding via contrastive learning. Notably, the dimensionality of the learned embedding is the same as that of the input data. This one-to-one correspondence between embedding and original features allows the definitions of importance scores for the original features. They evaluate their method quantitatively by the clustering performance on various datasets and compare to dimensionality reduction / visualization methods and a graph scattering approach.\n\n**Strengths:**\n\nS1 The one-to-one correspondence between original and embedding features allows to compute meaningful feature importance scores.  \nS2 The method explicitly incorporates several frequency scales and can thus capture both broad and fine data structure.  \nS3: The authors optimize hyperparameters for each competing method.  \n\n**Weaknesses:**  \n*Major:*   \nW1 Comparison mostly to visualization methods  \n\nMS-IMAP is not a dimensionality reduction or visualization method. Indeed, it maintains the same number of dimensions as the input had. Therefore, it is not meaningful to almost exclusively compare to dimensionality reduction / visualization methods, which only have 2-3 dimensions, when MS-IMAP has hundreds to thousands. If anything, the often slim margins by which MS-IMAP outperforms methods with order of magnitude less capacity cast doubt on the effectiveness of MS-IMAP. As the authors mention in Sec 5.2. the low embedding dimension of visualization methods leads to reduced expressiveness. As a result, I do not find the results in section 5.2 meaningful even if they were to show a clear performance increase for MS-IMAP. However, in many columns of Table 1 MS-IMAP gets outperformed by or draws with a method of significantly less capacity. \n\nI think a much more meaningful comparison would be to compare to unsupervised deep representations, e.g., via auto-encoders or self-supervised graph neural networks (keeping one or two dimensionality reduction methods in the comparison is fine, I think). Several of the dimension reduction methods allow embedding to higher-dimensional spaces: E.g. Diffusion Maps approximate the diffusion distance which is exactly represented by $N$-dimensional embeddings. Similarly, PHATE defines a graph distance, the potential distance, which is realized as Euclidean distance between points in $N$ dimensions, and then runs MDS to reduce the dimensionality. One could simply omit the MDS step, or use a much higher embedding dimension. One could run UMAP (with random initialization to avoid long Laplacian Eigenmaps computation) to much higher dimensions. Possibly, the methods from Gao et al 2018  2019 and Sun et al 2009 would also be relevant competitors.\n\nW2 Isolation of different features  \nI appreciate that during the optimization, only values for the same feature get combined. However, statements such as \n\n>  each embedding coordinate is updated solely based on the SGW coefficients of its corresponding input feature\n\nin the Discussion are too strong because the structural term $v_{ij}^{\\psi_x(s_k, :, :)}$ depends on all features jointly. In a similar vein, the graph structure itself depends on all features jointly. \n\n*Minor:*  \nW3 Expand on feature importance by using MS-IMAP selected features in other methods  \nAs I understand figures 2 and 3, they show the performance on a subset of the original features. This subset is selected by importance scores either computed on the original or the MS-IMAP features. For the MS-IMAP results, the selected features are again fed into MS-IMAP. I think it would be good to highlight how useful these features are more generally. For instance, one could compute the clustering performance of the MS-IMAP features without additional application of MS-IMAP or as input to some of the other methods considered in the paper (graph scattering, Diffusion Maps, t-SNE, UMAP). In addition, it would be great if the authors could comment qualitatively on the selected features. Which features get selected? Is it plausible that these carry semantic relevant to the clustering task?\n\nW4 Section 6 hard to understand and disconnected  \nI am not an expect on spectral graph wavelets and did not check the correctness of the theoretical claim Theorem 3. However, I can say that a non-technical, high-level, natural language explanation of why Theorem 3 exhibits desirable properties of MS-IMAP is needed. What do $\\Lambda$-sets mean for representation properties of an operator. Is a small / low $\\Lambda$-set desirable? If so, why? Overall this section feels somewhat disconnected from the rest of the paper.\n\n\nW5 Missing competitors   \nThe authors describe Gao et al. 2018, Gama et al. (2019), Sun et al. (2009) as related methods, because they compute higher-dimensional graph features. However, they only compare against Gama et al.'s method. Why are the other two not relevant for the empirical comparison?\n\nW6 No code was submitted.\n\n**Typos:**  \n- Full stops should be after the equations, not on the next line on page 4.\n- UMAP's ARI score on the Cora dataset is as good as MS-IMAP's in the shown digits, but not bold in Table 1\n- TriMap's and PaCMAP's ARI scores on the AWA dataset are equal to or better than the bold ones in that column of Table 1.\n- Appendix A.1 ends in the middle of a sentence (that was not capitalized)\n- \"Multi-Scale UMAP\" --> \"Mulit-Scale IMAP\" on page 21.\n- UMAP uses an exponential kernel, not a Gaussian kernel (Sec 2.1)\n\n**Questions:**  \n\nQ1 What does the acronym MS-IMAP stand for?\n\nQ2: What doe the following notations mean?\n\n$\\hat{f}$ in Eq (1)\n\n$\\tilde{\\psi}_x$ in Eq (2)\n\n$\\psi_{x_i}$ after Eq (3) (is it the same as $\\psi_x(:, :, i)$)?\n\nQ3: What is meant by weighing each scale by its energy concentration? Where can I find this in the equations describing the model?\n\nQ4: Are there any hyperparameters of MS-IMAP that get optimized per dataset? Which number of neighbor was used for the graph construction? Which wavelet scales were used?", "claims_and_evidence_yn": "No", "claims_explanation": "I think that several claims in the paper have problematic supporting evidence. The biggest issue to me is that MS-IMAP, which does not reduce the dimensionality at all, is almost exclusively compared to methods that reduce the dimensionality by several orders of magnitude to 2 or 3 dimensions (W1, W5). This is an apples-and-oranges comparison. Second, I think that there is more interaction between the different features (through the graph and the structural term), than the authors claim (W2). The theoretical section could be explained more accessibly and connected more closely to the proposed method (W4).", "audience_yn": "Yes", "audience_explanation": "The feature selection approach seems promising and interpretability is a pressing issue. The spectral graph wavelet approach is conceptually interesting, but need to be validated empirically in a more convincing way.", "requested_changes": "C1: Include several higher capacity competitors (e.g. Gao et al., (graph) autoencoders, deep (graph) SSL methods, high-dimensional UMAP, DiffMaps, PHATE variants) to determine if MS-IMAP outperforms models that are not handicapped by extremely low-dimensional representations (W1, W5).\n\nC2: Adapt the description of feature isolation (W2).\n\nC3: Expand the feature selection section (W3).\n\nC4: Explain the relevance of the theoretical section more clearly (W4).\n\nC5: Address typos and questions above.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "0Eccv2vOtM", "reviewer": "Reviewer_1k5z", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "# Summary\n\n**MS-IMAP** applies Spectral Graph Wavelets with contrastive learning to produce interpretable embeddings that preserve a one-to-one mapping between input features and embedding dimensions. It encodes multi-scale graph structures with 3D tensors, optimizes using SGD and cross-entropy loss, and provides theoretical grounding through Paley–Wiener spaces.  \nExperiments on five datasets demonstrate strong clustering performance and enable feature importance analysis using Laplacian Score (LS) and Mutual Information (MI).\n\n---\n\n# Strengths\n\n### Interpretability Contribution\nMS-IMAP’s one-to-one feature mapping fills a real gap in manifold learning. Unlike UMAP, t-SNE, or Laplacian Eigenmaps, it allows direct tracing from embedding dimensions to input features, which is especially valuable in domains like finance and biology.\n\n### Comprehensive Evaluation\nThe study includes thorough experiments with well-tuned baselines, diverse datasets across multiple domains, detailed ablation studies, and competitive results that are best or second-best on all benchmarks.\n\n### Dual Feature Importance Methods\nUsing both LS (unsupervised smoothness) and MI (task-aware discriminative power) provides complementary insights into feature relevance, as clearly discussed in Section 4.3.\n\n---\n\n# Weaknesses\n\n### Unfair Dimensionality Comparison\nMS-IMAP retains the full input dimension \\(D\\) (from 85 to 2,048), while baselines are reduced to two or three dimensions. This gives an inherent advantage in clustering performance and contradicts the manifold learning assumption that data lies in lower-dimensional spaces.\n\n### Circular Feature Importance Validation\nPseudo-labels are created using k-means on MS-IMAP embeddings, and MI is then computed between these embeddings and their own derived labels. This forms a circular validation process, meaning improvements mainly reflect self-consistency rather than independent verification.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The theoretical claims are rigorously proven (Theorems 2-3 with full proofs), clustering performance is demonstrated across 5 diverse datasets with proper statistical reporting and extensive baseline comparisons (9 methods with hyperparameter tuning), and the feature importance advantage is consistently shown in Figures 2-3, 5-6.", "audience_yn": "Yes", "audience_explanation": "The paper addresses interpretable manifold learning, a timely problem given increasing demand for explainable ML. The novel one-to-one feature correspondence, rigorous theoretical foundation connecting harmonic analysis to graph learning, multi-scale representation approach, and results on practical datasets across multiple domains (biology, finance, vision) would appeal to interpretable ML researchers, graph learning community, and practitioners in high-stakes applications", "requested_changes": "Please refer to the weakness.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "0cP0ygF7If", "reviewer": "Reviewer_6YrC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "1. This paper introduces a framework for multi-scale graph representation using a 3D tensor and contrastive learning.\n2. This paper proposes a manifold-based interpretation method for feature importance.\n3. This paper uses spectral graph wavelets to capture both global and local graph structures.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The paper provides extensive experimental validation on multiple public datasets to support its claims. It compares the performance of MS-IMAP against several baselines. The results show that MS-IMAP consistently outperforms the baselines, with the most significant improvements noted in feature importance evaluation. It also includes a theoretical justification for its approach, particularly for the use of spectral graph wavelets over traditional Laplacian operators.", "audience_yn": "Yes", "audience_explanation": "This paper introduces MS-IMAP, a graph embedding framework that addresses the common interpretability problem in manifold learning by combining spectral graph wavelets and contrastive learning. The core contribution is its ability to establish a one-to-one correspondence between embedding dimensions and original features, allowing for direct assessment of feature importance. Additionally, MS-IMAP uses spectral graph wavelets to capture both global and local graph structures, generating more expressive embeddings, while an efficient 3D tensor optimization method enhances its computational efficiency.", "requested_changes": "The paper could be improved by revising its language, which at times reads like it was generated by a large language model and lacks depth. The authors should also pay attention to small formatting details, such as correctly punctuating inline formulas as part of a sentence.\n\nRegarding the common challenge of selecting the number of dimensions for spectral clustering without prior information, the authors could strengthen their work by discussing and applying well-established model selection techniques. Promising methods include the eigengap heuristic, which finds the largest gap between consecutive eigenvalues of the graph Laplacian, and the silhouette score or Calinski-Harabasz index, both of which evaluate cluster quality for a range of k values. Incorporating an experiment using these methods would demonstrate a more robust and practical approach to their framework.", "broader_impact_concerns": "I have no concern.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)