Distributional Reduction: Unifying Dimensionality Reduction and Clustering with Gromov-Wasserstein
Authors: Hugues Van Assel
Cédric Vincent-Cuaz
Nicolas Courty
Rémi Flamary
Pascal Frossard
Titouan Vayer
Venue: Accepted by TMLR ()

Abstract: Unsupervised learning aims to capture the underlying structure of potentially large and high-dimensional datasets. Traditionally, this involves using dimensionality reduction (DR) methods to project data onto lower-dimensional spaces or organizing points into meaningful clusters (clustering). In this work, we revisit these approaches under the lens of optimal transport and exhibit relationships with the Gromov-Wasserstein problem. This unveils a new general framework, called distributional reduction, that recovers DR and clustering as special cases and allows addressing them jointly within a single optimization problem. We empirically demonstrate its relevance to the identification of low-dimensional prototypes representing data at different scales, across multiple image and genomic datasets.

Published in Transactions on Machine Learning Research (06/2025)
Distributional Reduction: Unifying Dimensionality Reduc- tion and Clustering with Gromov-Wasserstein
Hugues Van Assel
hugues. van asselßens-lyon.fr
ENS de Lyon, CNRS, UMPA UMR 5669
Cédric Vincent-Cuaz
cedric. cadrrc.oincenttccuuzeee.......
EPFL, LTS4
Nicolas Courty
courty@hnio-uube.u
Université Bretagne Sud, TRISA UMR 6074
Rémi Flamary
École polytechnique, Institut Polgtechnique de Paris, CMAP UMR 7641
Pascal Frossard
frssalffoossard@eepff..t
EPFL, LTS4
Titouan Vayer
titouan vatiannvayer@inrii...
Inria, ENS de Lyon, CNRS, Université Claude Bernard Lyjon LIP UMR 5668
Reviewed on OpenReview: https: // opwpn/vveeeee...ee./ net/ formm cm683354
Abstract
Unsupervised learning aims to capture the underlying structure of potentially large and high-dimensional datasets. Traditionally, this involves using dimensionality reduction (DR) methods to project data onto lower-dimensional spaces or organizing points into meaningful clusters (clustering). In this work, revisit these approaches under the lens of optimal transport and exhibit relationships with the Cromov-Wasserstenn problem. This unveils new general framework, called distributional reduction, that recovers DR and clustering as special cases and allows addressing them jointly within a single optimization problem. We empirically demonstrate its relevance to the identification low-dimensional prototypes representing data different scales, across multiple image and genomic datasets.
1 Introduction
One major objective of unsupervised learning (Hastie et al. 2009 is to provide interpretable and meaningful approximate representations of the data that best preserve its structure í.e. the underlying geometric relationships between the data samples. Similar essence to Occam's principle frequently employed supervised learning, the preference for unsupervised data representation often aligns with the pursuit of simplicity, interpretability or visualizability the associated model. These aspocts are determinant in many real-world applications where the interaction with domain experts is paramount for interpreting the results and extracting meaningful insights from the model. For instance, the design of tissue atlases the inference of single-cell trajectories, which are essential for addressing various biological challenges (Rao et al. 2021 Saelens et al. 2019), depends on the analysis of metacells, i.e. granular and interpretable representations of cells that enable separating sampling effect from biological variance (Baran et al. 2019)
Dimensionality reduction and clustering. When faced with the question of extracting interpretable representations, from a dataset community has proposed variety of methods. Among them, dimensionality reduction (DR) algorithms
XN) Nxx of N samples in R", the machine learning

---

Published in Transactions on Machine Learning Research (06/2025)
Eaisediing
Iupur X
Emiishhing 2
Z (auzean
Caiz Lorentz)
Cx(3) (SEA)
Cx(X) (SEA)
Compling
Coupling
Figure 1: Illustration of our DistR method on two toy examples. For each example (left and right), the data pairs X and Z represent the original data and their embedding as found by our method DistR. These pairs are associated to similarity matrices Cx(X) E NxN and Cz(Z)= E guxx (n < N), which encode the pairwise similarities in their respective spaces. In the left example, we consider Gaussian similarity for Z and Symmetric Entropic Affinity (SEA) (Van Assel et al, 2023) for X. In the right example, we consider SEA for X and similarity in the Lorentz model for Z (Nickel and Kiela 2018). DistR also provides a coupling T between the points in X and Z, illustrated in purple, with its marginals.
have been widely used to summarize data in a low-dimensional space Z = (Z), Zn E Rx with d & p. allowing for visualization of every individual points for small enough d (Lee et al.) 2007: Van Der Maaten ct al. 2009). Another major approach is to cluster the data into n groups, with n typically much smaller than N, and to summarize these groups through their centroids (Bishop, 2006; Von Luxburg 2007.. Clustering is particularly interpretable since it provides a smaller number of representative points that can be more easily inspected. The cluster assigmments can also be analyzed. Both DR and clustering follow a similar philosophy of summarization and roduction of the dataset using a smaller size representation.
Two sides of the same coin. As matter of fact, methods from both families share many similitudes, including the construction of similarity graph between input samples. In clustering, many popular approaches design a reduced or coarsened version of the initial similarity graph while preserving some of its spectral properties (Von Luxburg, 2007 Schaeffer 2007). In DR, the goal is to solve the inverse problem of finding low-dimensional embeddings that generate a similarity graph close to the one computed from input data points (Ham et al. 2004; Hinton and Roweis 2002). Our work builds on these converging viewpoints and addresses the following question: can DR and clustering be expressed in common and unified framework ?
A distributional perspective. To answer this question, we propose to look at both problems from distributional point of view, treating the data an empirical probability distribution μ NV Σ, 6x;- This enables to consider statistical measures of similarity such as Optimal Transport (OT), which is at the core of our work. On the one hand, OT and clustering are strongly related. The celebrated k-means algorithm can be soen as a particular case of minimal Wasserstein estimator where a distribution of n Diracs is optimized w.r. their weights and positions (Canas and Rosasco 2012). Other connections between spectral clustering and the OT--baedd Gromov- Wassersten (GW) distance have been recently developed in Chowdhury and Needham (2021): Chen et al. (2023) : VVincent-Cuaz et al. (2022a). On the other hand, the link between DR and OT is less explored. DR. methods, when modeling data as distributions, usually focus joint distribution between samples within each space separately, see e.g. [an Assel et al. (2023) or Tu et al. (2019). Consequently, they do not consider couplings to transport samples across spaces of varying dimensions. At the time of this paper's submission, we note that other authors have independently developed a similar line of work, as detailed in Clark et al. (2024); Murray and Pickarski (2024). Their studies focus on "muliidimensional scaling" problems and establish connections to the semi-rclaxed GW framework but do not explore alternative DR approaches or clustering aspects.
Contributions. In this paper, we develop general distributional framework that encompasses both DR and clustering as special cases. We notably cast those problems as finding a reduced distribution that minimizes the GW divergence from the original empirical data distribution. Our method proceeds by first constructing an input similarity matrix Cx(X) that is matched with the embedding similarity Cz(Z) through an OT coupling matrix T. The latter establishes correspondences between input and embedding samples. We illustrate this principle in Figure where one can notice that (Z) preserves the topology of Cx CX) with a reduced number of nodes. The adaptivity of our model that can select an effective number of cluster n,
2

---

Published in Transactions on Machine Learning Research (06/2025)
is visible in the bottom plot, where only the exact number of clusters in the original data (9 out of the 12 initially proposed) is automatically recovered. Our method can operate in any embedding space, which is illustrated by projecting in either 2D Euclidean plane or a Poincaré ball.
Outline. We show that this framework is versatile and allows recovering as special cases many popular DR methods such as the kernel PCA and neighbor embedding algorithms, but also clustering algorithms such as weighted k-means and its kernel counterpart including spectral clustering (Chan et al. 2004 Dhillon et al. 2007). We first prove in Section 3 that DR can be formulated as GW projection problem under some conditions the loss and similarity functions. We then propose in Section 4 novel formulation of data summarization as a minimal GW estimator that allows selecting both the dimensionality of the embedding d (DR) and the cardinality of the support n (Clustering). Finally, we show in section the practical interest of our approach on images and genomics datasets.
Notations. The 8t entry of a vector V denoted as either v, or [v], Similarly, for a matrix M, My and My both denote its entry (i.j). SN is the set of permutations of [N] Pv((R) refers to the set of discrete probability measures composed of N points of Rd. Σn stands for the probability simplex of size N that is IN := {h E R s.t. E, h4 1}. IMg(M) and exp(M) are to be understood element-wise. For x E IR diag(x) denotes the diagonal matrix whose elements are the 24.
2 Background on Dimensionality Reduction and Optimal Transport
We first review the most popular DR approaches and introduce the Gromov- Wasserstein problem.
2.1 Unified View Dimensionality Reduction
Let X (X1,..,,xy) X E Nxp be an input dataset. Dimensionality reduction focuses on constructing a low- dimensional representation or embedding = (Z,, =N)
zN) T RNxxd, where < p. The latter should preserve prescribed geometry for the dataset encoded via symmetric pairwise similarity matrix Cx R XN obtained from X. To this end, most popular DR methods optimize Z such that a certain pairwise similarity matrix in the output space matches Cx according to some criteria. We subsequently introduce the functions
C + RNNN C2 Nxd RNNN
(1)
which define pairwise similarity matrices in the input and output space, from the dataset X and the embedding Z respectively. The DR problem can be formulated quite generally the optimization problem
min xd
Σ ([Cx(X))s,[C2( (K))E[NT
(DR)
ZERvxd
where L: X IR + R is a loss that quantifics the discrepancy between similaritics between points in the input space IRP and in the output space Rd. Various losses are used, such as the quadratic loss L2(2, y) = (x - y or the Kullback-Leibler divergence T (I. y) = log(æ/u) - x + y. Below, we recall several popular methods that can be placed within this framework.
recovers spectral methods by choosing the quadratic loss I, = T2 and Cz(Z) ((z;, 21)) (E the matrix of inner products in the embedding space. Indeed, in this case, the objective value of eq. (DR) reduces to
Spectral methods. When Cx(X) is a positive semi-definite matrix, eq. DR
Σ L2(1Cx(X)ss (z, = [X) - ZZ T 17: (NJ)
where ||=]|= is the Frobenius norm. This problem is commonly known as kernel Principal Component Analysis (PCA) (Schölkopf et al. 1997) and an optimal solution given by Z = =(ViVI,,... V ( where X, is the i-th largest eigenvalue of (X) with corresponding cigenvector v, (Eckart and Young 1936). As shown by Ham et al. (2004); Ghojogh et al. (2021), numerous dimension reduction methods can be categorized in this manner. This includes PCA when Cx (X) =XXT is the matrix of inner products in the input.
3

---

Published in Transactions on Machine Learning Research (06/2025)
space; (classical) multidimensiooall scaling (Borg and Groenen 2005 when Cx Cx(X)= HD with Dx the matrix of squared euclidean distance between the points in I&? and H = IN N I Allly is the centering x H with D the geodesic distance matrix;
matrix; Isomap Tenenbaum et al. 2000 with Cx(X) = - HD (1) Laplacian Eigenmap (Belkin and Niyogi 2003), with Cx(X) L& the pseudo-inverse of the Laplacian associated to some adjacency matrix Wx but also Locally Linear Embedding (Roweis and Saul 2000), and Diffusion Map (Coifman and Lafon 2006) (for all of these examples we refer to Ghojogh et al. 2021 Table 1).
-
Neighbor embedding methods. An alternative group of methods relies on neighbor embedding techniques which consists in minimizing in Z the quantity
Σ (is))[N]
Lxl((Cx(X)); (C2()))))))))()))))))
(NE)
Within framework, this corresponds to eq. DR with L LKL. The objective function of popular methods such as stochastic neighbor embedding (SNE) (Hinton and Roweis 2002 t-SNE (Van der Maaten
and Hinton 2008) can be derived from eq-
NE with particular choice of Cx. Cz. For instance SNE and t-SNE both consider in the input space a symmetrized version of the entropic affinity (Vladymyrov and
Carrcira- Perpinan 2013 Van Assel ct al. 2023). In the embodding space, (Z) is usually constructed from
8. "kernel" matrix Kz which undergoes a scalar
Van der Maaten and Hinton 2008), row-stochastic (Hinton
and Roweis 2002 or doubly stochastic (Lu et al.
2019 Van Assel et al. 2023) normalization. Gaussian kernel
Kz v exp( z. 2, 11,,, or heavy-tailed Student-t kernel [Kz]aj = (1+ Iz - z =,115) (113)-1, are typical choices
Van der Maaten and Hinton 2008). We also emphasize that one can retrieve the UMAP objective (Melnnes
et al. 2018) from eq. (DR) using the binary cross-entropy loss T.CCE(2, v)
log = + (1 x.) log I- A
R
comprehensive overview and probabilistic analysis of these methods can be found in Van Assel et al. (2022)
Remark 2.1. The usual formulations of neighbor embedding methods rely on the loss L(x, y) = x log(z/
instead of LKL. However, due to the normalization, the total mass Ex,[C2 (Z)]i constant (often equal to 1) in all of the cases mentioned above. Thus the minimization in Z with the LKt. formulation is equivalent to the usual formulations of neighbor embedding objectives.
Non-Euclidean geometries. Most DR methods can also be extended to incorporate -on-Euclidean geometries. For instance, Hyperbolic spaces (Chami et al. 2021 Fan et al. 2022 Guo et al., 2022 Lin et al., 2023) are of particular interest as they can capture hierarchical structures more effectively than Euclidean spaces and further mitigate the curse of dimensionality. These methods introduce additional hyperparamelers and optimization difficulties that can limit their applicability. For this reason, we report some numerical experiments with such kernels in Appendix F.10 as proofs of concopt supporting the versatility of our method.
2.2 Optimal Transport Across Spaces
Optimal Transport (OT) (Villani et al., 2009 Peyré et al., 2019 is a popular framework for comparing probability distributions and is at the core of our contributions. We review in this section the Gromov-
Wasserstein formulation of OT aiming at comparing distributions across spaces.
Wassert-nss (GW). The GW framework (Mémoli 2011 Sturm 2012) comprises collection of OT methods designed to compare distributions by examining the pairwise relations within each space. For two matrices C E RNNN C E Rax*, and weights h E ΣN, h E Σ,, the GW discrepancy with inner loss L. (Peyré et al. 2016 is defined as
CW,(C,C,h,b)
min C,(,, Tel((hh)
(GW)
with (T(C,,,T)
Σ LCC,,Co)TTa ajkl
and U(, h) {T E RN T1, h, h} is the set of couplings between h and h. In this formula- tion, both pairs hC and hC,h) can be interpreted as graphs with corresponding connectivity matrices C,C, C, and where nodes are weighted by histograms h, (with implicit supports). Equation (GW) is thus a quadratic problem (in T) which consists finding a soft-assignment matrix T that aligns the nodes of the two graphs in a way that preserves their pairwise connectivities.
4

---

Published in Transactions on Machine Learning Research (06/2025)
From distributional perspective, GW defines distance between distributions that do not belong to the same metric space. For two discrete probability distributions μxx = -i== ( E PPy(R") Mz = ( E (IR)) and pairwise similarity matrices Cx(X) and Cz associated with the explicit supports X == (x1.." xx)) and Z = (Z1.-... . the quantity GC (X), Cz( hx , hz) is a measure of dissimilarity or discrepancy between μx, μzz- Specifically, when == L2g, and Cx (X), Cz(Z) are pairwise distance matrices, GW defines proper distance between jux and μzz with respect to measure preserving isometries (with weaker assumptions on Cx, Cz, GW defines a pseudo-metric w.r.t. a different notion of isomorphism (Chowdhury and Mémoli 2019), detailed in Appendix C
Due to its versatile propertics, notably in comparing distributions over different domains, the GW problem has found many applications in machine learning, e.g., for 3D meshes aligmment (Solomon et al. 2016 Ezuz
et al. 2017). NLP [Alvarez-Melis and Jaakkola 2018). (co-)clustering (Peyré et al. 2016 Redko et al. 2020).
single-cell analysis (Demetei et al. 2020), neuroimaging (Thual et al. 2022), graph representation learning
Xu 2020: Vincent- Cuaz al. [021 Liu et al., 2022b Vincent- Vincent-Cuaz et al.) 2022b Zeng et al.) 2023 and
parlitioning (Xu et. al. 2019; Chowdhury and Needham 2021]. In this work, we leverage the GW discrepancy to extend classical DR approaches, framing them as the projection of distribution onto space of lower dimensionality.
3 OT Formulation of DR
In this section, we outline the strong connections between the classical eq.
(DR) and the GW problem.
Gromov-Monge interpretation of DR. As suggested by eq. (DR dimension reduction seeks to find embeddings Z so that the similarity between the (i.j) samples of the input data as close as possible to the similarity between the (i,5) samples of the embeddings. Under equivariant assumptions on Cz, this also amounts identifying the embedding and the best permutation that realigns the two similarity matrices Cx(X) and CZ) Recall that the function Cz is equivariant by permutation, if, for any N x N permutation matrix P and any Z, Cz(PZ) PCz(()) (Bronstein et al. 2021). This type of assumption is natural for Cz: if we rearrange the order of samples (i.e. the rows of Z). we expect the similarity matrix between the samples to undergo the same rearrangement.
Lemma 3.1. Let Cz be a permutation equivariant function and L any loss. The minimam of eq. equal to
DR is
min min ZER GSN
(2)
Also, any sol. Z of eq. (R) is such that (Z. id) is solution of eq. And conversely any (Z.o) sol. of eq. is such that Z is a solution of eq. (R) up to a.
Lemmn proven in Appendix B.1 establishes an equivalence between eq. (R) and the optimization the embedding Z w..t a Gromov-Monge discrepancy (Mémoli and Needham 2018 given in eq. (2 which seeks to identify the permutation a that best aligns two similarity matrices, by solving a combinatorial quadratic assignment problem (Cela 2013). We can delve deeper into these comparisons and demonstrate that the general formulation of dimension reduction is also equivalent to minimizing the Gromov- Wasserstein objective, which serves as a relaxation of the Gromov-Monge problem Mémmoli and Needham 2022).
DR as GW Minimization. We suppose that the distributions have the same number of points (N = n) and uniform weights (hz = hx 1 Ix). We recall that a matrix C E NNNN is conditionally positive definite (CPD), resp. conditionally negative definite (CND), if it is symmetric and Vx R. s.t. x =11N have Cx 0, resp. < 0.
We thus extend Lemma 3.1 to the GW problem in the following (proof in Appendix B.2
Theorem 3.2. The minimum eq. DR is equal to minz GW Cz(Z), IN, Ix) in the following cases:
(1) (spectral methods) When Cx (X) is any matria, L = La and Cz -ZZT
5

---

Published in Transactions on Machine Learning Research (06/2025)
(ti) (neighbor embedding methods) When Im(Cx) C R L LKL, the matrir Cx(X) is CPD and, for any Z,
Cz(Z) = diag( Kz diag(3z))
for some ozz, Bz E R and K2 E R >0
Nx. such that log K2) is CPD.
Remarkably, the first item of Theorem 3.2 shows that all spectral DR methods can be seen as OT problems in disguise, as they all equivalently minimize a GW problem. The second item also provides some insights into this equivalence in the case of neighbor embodding methods that require more assumptions. For instance, the Gaussian kernel Zz, used extensively in DR (Section 2.1 , satisfies the hypot hesis log( (gg( --11:-2,,11ā); is CPD (see e.g- Maron and Lipman 2018)- The terms oz 8z also allow for considering all the usual normalizations of Kz: by scalar so as to have [ [Z((Z) = but also any row doubly stochastic normalizations (with the Sinkhorn-Knopp algorithm Sinkhorn and Knopp 1967).
Matrices satisfying log (z) being CPD are well-studied in the literature and are known as infinitely divisible matrices [Bhatia 2006). It is noteworthy that the t-Student kernel does not fall into this category. Moreover, in the aforementioned neighbor embedding methods, the matrix Cx(X) is generally not CPD. The intriguing question of generalizing this result with weaker assumptions on Cz and Cx remains open for future research. Interestingly, have observed in the numerical experiments performed in Section 5 that the symmetric eutropic alfinity of Van Assel et al. (2023) was systemalically CPD.
Remark 3.3. In Corollary B.3 of the appendix we also provide other sufficient conditions for neighbor embedding methods with the cross-entropy loss L(x, y) = x log(æ/y). They rely on specific structures for Cz but do not impose any assumptions on Cx. Additionally, in Appendix B.3] we provide necessory condition based a bilincar relaxation of the GW problem. Although its applicability is limited due to challenges in proving it in full generality, it requires minimal assumptions Cx, Cz and L.
In essence, both Lemma 3.1 and Theorem 3.2 indicate that dimensionality reduction can be reframed from a distributional perspective, with the search for an empirical distribution that aligns with the data distribution in the sense of optimal transport, through the GW framework. In other words, DR is informally solving min GW(, n D 8x.9 n E''
2i=16,)) 62,).
4 Distributional Reduction
The above perspective DR allows for the two following generalizations. Firstly, beyond solely determining the positions z, of Diracs (as in classical DR) we can now optimize the mass of the distribution μzz. This is interpreted as finding the relative importance of each point in the embedding Z. More importantly, due to the flexibility of GW, we can also seek a distribution in the embedding with a smaller number of points n < N. This will result in both reducing the dimension and clustering the points in the embedding space through the optimal coupling. Informally, we thus propose a new Distributional Reduction (DistR) framework that aims at solving min (2P. GW() Σ 8x,,112)
4.1 Distributional Reduction Problem
Precisely, the novel optimization problem that we lackle this paper can be formulated as follows
min GWz(Cx(X), Cz(), bx, hz) z hz
(DistR)
This problem comes down to learning the graph (Z), hz) parametrized by Z that is the closest from (Cx (C), in the GW sense. When n < N, the embeddings then act as low-dimensional prololypical exomples of input samples, whose learned relative importance hy accommodates clusters of varying proportions in the input data X (soc Section I We refer them as prototypes. The weight vector hx typically assumed be uniform, that is hx \/2 1N, in the absence of prior knowledge. As discussed in Section traditional DR amounts to setting n N, hz IN

---

Published in Transactions on Machine Learning Research (06/2025)
PCA
srGW on grid 16*16
srGW on grid 32*32
srGW on grid 64"64
srGW on grid 96*96
y-0
Y-1
Figure 2: GW projections of a genomics dataset Chen et al. 2019 on regular grids with increasing resolutions, respectively encoded as Cx(X) XX" and (z ZZ Pixels on cropped grids are colored by interpolating samples' colors according to the transport plan and their intensity is proportional to their mass.
One notable aspect of our model is its capability to simultaneously perform DR and clustering. Indeed, the optimal coupling T E [0, N of problem eq. (DistR is, by construction, a soft-assignment matrix from the input data to the embeddings. In Section 4.2 we explore conditions where these soft assignments transform into hard ones, such that each point is therefore linked to unique prototype/cluse /cluster.
allows each point x; to be linked to one or more prototypes žj (clusters).
A semi-relaxed objective. For given embedding Z and L= L2, it is known that minimizing the DistR objective w.r.t by is equivalent to problem that is computationally simpler than the usual GW one, namely the semi-relaxed GW divergence srGW (Vincent-Cuaz et al. 2022a given by
min Cz((x(X),C2(2))17 (hx)
(srGW)
where U₁ (hx) = (T E RN*x" T1, = hx}. To efficiently address eq. DistR we first observe that this equivalence holds for any inner divergence L. Additionally, we prove that srGW,, remains divergence soon as L is itself a divergence. Consequently, srGW, vanishes iff both compared measures are isomorphic in a weak sense (Chowdhury and Mémoli 2019). We emphasize that taking a proper divergence L is important (and basic assumptions on X), as it avoids some trivial solutions. These results are detailed in Appendix C
Interestingly, srGW projections, i.e. optimizing only the weights h2 over simple fixed supports Z, have already remarkable representational capability. We illustrate this in Figure 2 by considering projections of a real-world dataset over 2D grids of increasing resolutions. Sctting Cx (X) XX and C2(X) ZZ we can see that those projections recover faithful coarsened representations of the embeddinggs learned using PCA. DistR aims to exploit the full potential of this divergence by learning a few optimal prototypes that best represent the datasct.
Computation. R is a non-convex problem that we propose to tackle using a Block Coordinate Descent. algorithm (BCD, Tseng 2001) guaranteed to converge local optimum (Grippo and Sciandrone 2000; Lyu and Li 2023). The BCD alternates between the two following steps. First, optimize in Z for fixed transport plan using gradient descent with adaptive learning rates (Kingma and Ba 2014). Then solve for srGW problem given Z. To this end, we extended both solvers proposed in Vincent-Chure Cuaz et al. (2022a) to support losses L2, LKL and LBCE (see Appendix E.I. Namely, the Mirror Descent promoting smoothness along iterations depending on a hyparameter E > and the Conditional Gradient seen as the edge-case E = 0. (Peyré et al.) 2016), vanilla implementation leads to M(nN++n2N operations to compute the loss or its gradient. However in many DR methods, Cx(X) or Cz(Z), or their transformations within the loss L, admit explicit. low-rank factorivations. Including e.g. matrices involved in spectral methods and other similarity matrices derived from squared Euclidean distance matrices (Scetbon et al. 2022). In these settings, we exploit these factorizations to reduce the computational complexity of our solvers down to O(Nn(p + d) + (N + n)pd + min(n", (ndd) when L = L2, and O(Nnd + + when L = LKLL We refer the reader interested in these algorithmic details to Appendix E As shown in Appendix F.9 DistR already runs
Following Proposition
in a few seconds on datasets with several thousands of samples without leveraging low-rank properties.
Related work. The CO-Optimal Transport (COOT) clustering approach proposed in Redko al. 2000 is the closest to our work. It simultaneously estimates sample and feature clustering and thus be
7

---

Published in Transactions on Machine Learning Research (06/2025)
directly applied to perform joint clustering and DR. Despite SOTA performances on co-clustering tasks, COOT-clustering corresponds to linear DR method which is limiting for joint clustering and DR tasks. In Appendix F.2 we show that COOT leads to a low average homogeneity of the prototypes {zs he indicating that it often assigns points with different labels to the same prototype. In contrast, DistR exploits the more expressive non-linear similarity functions offered by existing DR methods and leads to significantly improved results. Other joint DR-clustering approaches, such as Liu et al. (2022a), involve modeling latent variables by a mixture of distributions. In comparison, DistR more versatile, as it easily adapts any (L, Cx. Cz) of existing DR methods.
4.2 Clustering Properties
We elaborate now on the links between DistR and clustering methods. In what follows, we call coupling T E [0, 1 with a single non-null element per membership matriz. When the coupling is membership matrix each data point is associated with a single prototype thus achieving hard clustering of the input samples. We will see that a link can be drawn with kernel k-means using the analogy of GW barycenters. More precisely the srGW barycenter Vincent-Chun et al. 2022a) seeks for a closest target graph (C, h) from (Cx, hx) by solving
min min CEREXX (hx)
(X),
(srGWB)
We stress that the only (important) difference between eq. srr SrGB and eq- (Dist is that there is no constraint imposed on C in srGWB. In contrast, eq. (DistR) )
looks for minimizing over C E (Zz(Z): Z E
RNxd). For instance, choosing Cz(Z)- = ZZ T in eq. Dist is equivalent to enforcing rank CC) I d in
oq. (rr
We establish below that srGWB is of particular interest for clustering. The motivation for this arises from the findings of (Chen et al. 2023), which demonstrate that, when Cx(X) is positive semi-definite and T is constrained to belong to the set of membership matrices (as opposed to couplings in U2 (h)), eq. 4srGWB) is equivalent to kernel k-means whose samples are weighted by hx (Dhillon et al. 2004 2007). Interestingly, these additional constraints are unnecessary. Indeed, we show below that the original srGWB problem admits membership matrices as the optimal coupling for a broader class of Cx input matrices for L = L2 (see proof in Appendix D).
Theorem 4.1. Let hx E EN and L = L2. Suppose that for any X E NxP the matriz Cx is CPD or CND. Then the problem eq. (rGW admits a membership matria as optimal coupling, i.e. there is a Ey (Xx(),,, C) with only one non-zero value per row.
minimizer of T E (hxxh)) -
The implications of this theorem are as follows. First, as shown Appendix D given an optimal plan T, the coefficient (i,5) of the barycenter C can be written as n.n; [p [ex[[C(xx) (X) pg pn TpTTa where n4 = Ep Tpi- Consequently, when T is membership matrix, C represent graph weights, where each element (i.j) is the weighted sum of the original graph weights Cx(X) for nodes in the clusters of i and j. following standard edge contraction methods (Loukas, 2019). Second, Theorem 4.1 and relations proven in Chen et al. (2023) state that eq. (srGWB) is equivalent to the aforementionel kernel k-means when Cx is positive semi-definite This equivalence also emphasizes spectrum spectrum--peeservatin properties of the srGWB problem, whose optimal value satisfies FI, (C,(CX(X), (X),
EN Xp - Σリ v2, where { and {uz}ī= are the eigenvalues and C respectively, and satisfy the Pointcaré
sorted in descending order of diag(hx) -1/2 (ddxgXxigg(bxx"""""" Separation Theorem. Finally, as the (hard) clustering property holds for more generic types of matrices, namely CPD and CND, srGWB stands out as a fully-fledged clustering method. Although these results do not apply directly to DistR., except if e.g. the dimension is set to the unknown rank of the srGW barycenter, we argue that they further legitimize the use of GW projections for clustering. Interestingly, we also observe in practice that the couplings obtained by DistR are always membership matrices, regardless of Cz.
5 Numerical Experiments
In this section, we demonstrate the relevance of our approach for joint clustering and DR, over 8 labeled datasets detailed in Appendix F including: 3 image datasets (COTL-20 Nene et al. 1996. MNIST & fashion-
8

---

Published in Transactions on Machine Learning Research (06/2025)
DistR (ours)
DR C
C- DR
COIL MNIST
Fashion-MNIST
SNA1 SNA2
ZEISEL ZEISEL 2
PBMC
C, ber 1
SEA/ Student
1a ba
SEA Student
104
G 9
x* 4T
norm. homogeneity
norm. Romogenaa
norm, silhouette
norm. silhouette
norm. k-means NMI
norm. k-means NMI
Figure 3: Best trade-off between homogeneity vs silhouette (2 first plots), and homogeneity vs NMI (2 last plots). Scores are normalized in [0, 1] via min-max scaling over a dataset. Small markers are scores per dataset for 5 runs, while big ones are their mean. We illustrate the 20-80% percentiles of scores per method as a colored surface.
MNIST Xiao et al. 2017 and 5 genomic ones (PBMC Wolf et al. 2018 SNA & 2 Chen et al. 2019 and ZEISEL Zeisel et al. 2015). To this end, we examine how effectively DistR prototypes offer discriminative representations of input samples various gramularities, using evaluation system described below. Code is provided at https:///gtthul cow cwm
Benchmarked methods. Let us recall that any DR method presented in Section 2.1 is fully characterized by a triplet (L, Cz) of loss and pairwise similarity functions. Given such triplet, we compare our DistR model against sequential approaches of DR and clustering, namely DR then clustering (DR-C) and Clustering then DR (C-DR). These 2-step methods are representative of current practice in the field of joint clustering and DR (Baran et al. 2019). In both cases, the clustering step is performed with spectral clustering (a benchmark with Kmeans clustering is provided in Appendix F.7). The three methods produce sample-to-prototype association matrix (formally described in Appendix F.1 further used to evaluate performance. For the choices of (L, Cx, Cz), we experiment with both spectral and neighbor embedding methods (NE). For the former, we consider the usual PCA setting with d = 10. Regarding NE, we rely on the Symmetric Entropic Affinity (SEA) (T] from Van Assel et al. (2023) for Cx and the maraa-normalized student similarity for C2 (Van der Maaten and Hinton] 2008), setting the dimension d = as coincides with real-life visualization purposes while outperforming PCA-based approaches. We report. validated hyperparemeters (e.g perplezity, number of prototypes n) and implementation details in Appendix F.] and provide our code with the submission. Morcover, we report additional experiments using the vanilla tSNE and UMAP kernels in Appendix leading to similar behaviors than kernels reported in the main paper.
Evaluation metrics. We consider several evaluation metrics. Firstly, we wonder whether the prototypes are individually discriminating input samples {x:) w.r.t their class label. We consider the homogeneity score [Roeebberg and Hirschberg] (007) that quantifies to which extent each prototype z is associated with points in the same class. Secondly, we aim to evaluate whether the embedding space is also globally discriminant, i.e. that prototypes organize themselves to form clusters that match the classes, measured by two metries. i) We associate with each protolype a label using majority voting from the associated samples. Then we make use of the popular silhouette score (Rousseeuw) 1987) to evaluate whether prototypes' positions are correctly alignod with these labels. ti) As commonly done to assess DR representations (Huang et al., 2022). we also consider the downstream lask of clustering the prototypes using the k-means algorithm. We then assign to each (x,} the cluster of its associated prototype, acting as a predicted label, and compute the Normalized Mutual Information (NMI) [Kvålseth| 2017] between these labels and {v:}. More details about evaluation metrics can be found in Appendix F
Results. First, we illustrate in Figure 3 the best trade-off between the aforementioned metrics achieved by all methods. For each method and dataset, we considered the model maximizing the sum of the two
1 An enhanced version of the original ISN affinity that preserves entropy normalization during symmetrization.
9

---

150 Published in Transactions on Machine Learning Research (06/2025)
MNIST
Fasition-MNN -INIST
COIL
PBMC
2
3
2
2
2
8
Figure 4: Examples of embeddings produced by DistR, using the SEA similarity for Cx and the Student's kernel for Cz, respectively in I²² for the first three datasets and the Poincaré ball for the last one. Displayed images are medoids for each cluster i.e. arg max,[Cx (Cx(X)xxx)) for cluster k. The area of image is proportional to [hz]a-
COIL
MNIST
SNAREseq
30
Silhouet
20-
20
DistR (ours)
DR-cC
10
20
10
C->DR
0-
50
100
150
200
50
100
150
200
0
50
100
150
200
k-meaus NMI
85
60
55-
40
50
100
150
200
50
100
150
200
0
50
100
150
200
100
100
100
Homogenetty
95-
95
90
85
30
100
150
200
50
100
150
200
0
50
100
200
Number of prototypes
Number of prototypes
Number of prototypes
Figure 5: Scores (x100) w.r. with respect to the number of prototypes (in RR' using the t-SNE model: SEA similarity for Cx (Vn Assel al. 2023), Student's kernel for Cz and loss LKL.
normalized metrics to account for their different ranges. DistR, being present the top-right of all plots, provides on average the most discriminant low-dimensional representations endowed with simple geometry, seconded by C-DR. The significant variance of these dynamics for all methods emphasize the difficulty of performing jointly clustering and DR. Interestingly, we show in Appendices |.5 that DistR is the best suited to describing most datasets at various granularities low-dimension. For a given configuration, DistR leads the most consistent performance over n on average across datascts. We also check that this consistency holds for different dimensions d in Appendix F.S We then display in Figure the scores obtained with the neighbor embedding model for several datasets and
10

---

Published in Transactions on Machine Learning Research (06/2025)
across all considered number of prototypes Other kernels and datasets are illustrated in Appendix F.4 Interestingly, looking at the homogeneity score (third row), one can notice that DistR is at least as good as C-DR for grouping points of the same label, even though C-+RR performs clustering in the high dimensional input space. DistR is also generally better than sequential approaches at representing a meaningful structure of the protolypes in low-dimension, measured by the silhouette and k-means NMI scores. Therefore, our approach DistR scems to effectively achieve the best equilibrium between homogeneity and preservation of the structure of the prototypes. For qualitative illustration, also plot some embeddings produced by our method in Figure
6 Conclusion
By making connection between the GW problem and popular clustering and DR algorithms, we proposed a unifying framework denoted as DistR. DistR enables transforming any DR algorithm into a joint DR-clustering method that produces embeddings (or prototypes) with chosen gramularity. We believe that the versatility of the GW framework will enable new extensions in unsupervised learning. For instance, the formalism associated with (semi-relaxed) GW barycenters naturally enables addressing multi-view scttings with multiple unaligned inputs of different sizes. Other promising directions involve, better capluring multiple dependency scales in the input data by hierarchically adapting the resolution of the embedding similarity graph, enabling batch optimization of embeddings to operate over much larger datasets.
Acknowledgments
The authors are grateful to Aurélien Garivier for insightful discussions. This project was supported in part by the ANR project OTTOPIA ANR-20-CHIAA-000 and through the MATTER project A-R-23- 20---23---00000000 We gratefully acknowledge the support of the Centre Blaise Pascal's IT test platform at ENS de Lyon (Lyon, France) for Machine Learning facilities. The platform operates the SIDUS solution (Quemener and Corvellec 2013). This work benefited from state aid managed by the Agence Nationale de la Recherche under the France 2030 programme, reference ANR-23-TACT-00055 Finally, it received funding from the Fondation de l'École polytochnique and the Personalized Health and Related Technologies PHRT project 2022-644.
References
D. Mellis and T. S. Jaakkola. wassssewwasersteii alignment of word embedding spaces. arXiv preprint Xrx: 00013, 2018.
Y. Baran, A. Bercovich, A. Sebe-Pedros, Y. Lubling, A. Giladi, E. Chomsky, Z. Meir, M. Hoichman, A. Lifshitz, and A. Tanay. Metacell: analysis of single-cell rna-seq data using k-nn graph partitions. Genome biology, 2011):1-19, 2019.
G. Bécigneul and O.-F. Ganea. Riemannian adaptive optimization methods. arXio preprint urXiv:110.00066. 2018.
M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data represcntation. Neural computation, 15(6):1373 13991196,,, 2003.
R. Bhatia. Infinitely divisible matrices. The American Mathemalicul Monthly, 113(3):221 -23,2,,,,,, 2006.
G. Birkhoff. Tres observaciones sobre el algebra lineal. Univ. Nac. Tucuman, Ser. A, 5:147-14 154, 1946.
C. M. Bishop. Pattern recognition and machine learning. Springer google schola, 2:5-43, 2006.
I. Borg and P. J. Groenen. Modern multiddimensional scaling: Theory and applications. Springer Science & Business Media, 2005.
M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. Xixi preprint SSt:2104.13477 2021.
11

---

Published in Transactions on Machine Learning Research (06/2025)
G. Canas and L. Rosasco. Learning probability measures with respect to optimal transport metrics. In Aduances in Newral Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
L. Cao, D. McLaren, and S. Plosker. Centrosymmetric stochastic matrices. Linear and Multilinear Algebra, 70(3):449 1644 2022.
E. Cela. The quadratic assignment problem: theory and algorithms, volume Springer Science & Business Media, 2013.
I. Chami, A. Gu, D. Nguyen, and C. Ré. Horopca: Hyperbolic dimensionality reduction via horospherical projections. In International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 1419 1429. PMLR, 2021.
E. Y. Chan, W. K. Ching, M. K. Ng, and J. Z. Huang. An optimization algorithm for clustering using weighted dissimilarity mcasures. Pattern recognition, 37(5):943-995 952, 2004.
S. Chen, B. B. Lake, and K. Zhang. High-throughput sequencing of the transcriptome and chromatin accessibility in the same cell. Nature biotechnology, 37(11):1452-1453 2019.
Y. Chen, R. Yao, Y. Yang, and Chen. A gromov-wasserstee geometric view of spectrum-preeervvigg graph coarsening, orXiv preprint ar Xi0:22306008844 2023.
S. Chowdhury and F. Mémoli. The gromov wasserstein distance between networks and stable network invariants. Information and Inference: A Journal of the IMA, 8(4):757 787, 2019.
S. Chowdhury and T. Neodham. Generalized spectral clustering via gromov-wasserstein learning. International Conference on Artificial Intelligence and Statistics, pages 712-720. PMLR, 2021.
R. A. Clark, T. Needham, and T. Weighill. Generalized dimension reduction using semi-relaxed gromov- wasserstein distance. arXiv preprint arXiv:2405.1500 2024.
R. R. Coifman and S. Lafon. Diffusion maps. Applied and computational harmonic analysis, 21(1):5-30, 2006.
P. Demetci, R. Santorella, B. Sandstede, W. S. Noble, and Singh. Gromov-wasserstein optimal transport align single-cell multi-omics data. BioRaie, pages 2020 04, 2020.
N. S. Detlefsen, J. Borovec, J. Schock, A. H. Jha, T. Koker, L. Di Liello, D. Stancl, C. Quan, M. Grechkin, and W. Falcon. Torchmetrics measrriin reproducibility in pytorch. Journal of Open Source Software, 7 (70):4101, 2022.
S. Dhillon, Y. Guan, and B. Kulis. Kornel k-mcans: spectral clustcring and normalized cuts. In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 551-556, 2004.
I. S. Dhillon, Y. Guan, and B. Kulis. Weighted graph cuts without eigenvectors multilevel approach. IEEE transactions on pottern analysis and machine intelligence, 29(11):1944-1957, 2007.
C. Eckart and G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1()): 211 218, 1936.
D. Ezuz, J. Solomon, V. G. Kim, and M. Ben-Chen. Gwenn: A metric alignment layer for deep shape analysis. Tn Computer Graphics Forum, volume 36, pages 49-57. Wiley Online Library, 2017.
X. Fan, C.-H. Yang, and B. C. Vemuri. Nested hyperbolic spaces for dimensionality reduction and hyperbolic nn design. In Conference on Computer Vision and Pattern Recogmition, pages 356 365, June 2022.
Flamary, N. Courty, A. Gramfort, M. Alaya, A. Boisbunon, Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, al. Pot: Python optimal transport. The Journal of Machine Learning Research, (111))333133355 3578, 2021.
12

---

Published in Transactions on Machine Learning Research (06/2025)
B. Ghojogh, A. Ghodsi, F. Karray, and M. Crowley. Unified framework for spectral dimensionality reduction, maximum variance unfolding, and kernel learning by semidefinite programming: Tutorial and survey, arviv preprint Xri::22106. 15379, 2021.
L. Grippo and M. Sciandrone. On the convergence of the block nonlinear gauss seidel method under convex constraints. Operations rescorch letters, 26(3):127-136, 2000.
Y. Guo, H. Guo, and
Yu. CO-SNE: Dimensionality reduction and visualization for hyperbolic data. In Computer Vision and Pattern Recognition (CVPR). pages 11-20, Los Alamitos, CA, USA, jun 2022.
J. Ham, D. D. Lee, S. Mika, and B. Schölkopf. A kernel view of the dimensionality reduction of manifolds. In Proceedings of the twenly-first international conference on Machine learning page 47, 2004.
T. Hastic, R. Tibshirani, and J. Friedman. Unsupervised Learning. Springer New York, New York, NY, 2009.
G. E. Hinton and S. Rowcis. Stochastic neighbor embedding. Advances in neural information processing systems, 15, 2002.
H. Huang, Y. Wang, C. Rudin, and E. P. Browne. Towards a comprehensive evaluation of dimension reduction methods for transcriptomic data visualization. Communications biology, 5(1):719, 2022.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. ar preprinl X12::1412.66950, 2014.
M. Kochurov, R. Karimov, and S. Kozlukov. Geoopt: Riemannian optimization pytorch, 2020.
H. Konno. A cutting plane algorithm for solving bilinear programs. Mathematicaal Programming, 11(1):14-27 1976.
T. O. Kvålseth. On normalized mutual information: measure derivations and properties. Entropy, 19(1,),331, 2017.
S. Lacoste-Julien. Convergence rate of frank-wolfe for non-convex objectives, ar preprint arXiv: wrXiv:1607.003// 2016.
J. A. Lee, M. Verleysen, et al. Nonlinear dimensionality reduction, volume 1. Springer, 2007.
Y.-. E. Lin, R. R. Coifman, G. Mishne, and R. Talmon. Hyperbolic diffusion embedding and distance for hicrarchical reprosentation learning. In International Conference Machine Learning. 2023.
W. Liu, X. Liao, Y. Yang, H. Lin, J. Yeong, X. Zhou, X. Shi, and J. Liu. Joint dimension reduction and clustering analysis single-cell rna-seq and spatial transcriptomics data. Nucleic acids research, 50(12): e72 e72, 2022a.
W. Liu, J. Xie, C. Zhang, M. Yamada, N. Zheng, and H. Qian. Robust graph dictionary learning. In The Eleventh International Conference on Learning Representations, 2022b.
A. Loukas. Graph reduction with spectral and cut guarantees. Journal of Machine Learning Research, 20 (116):1-42, 2019.
Y. Lu, J. Corander, and Z. Yang. Doubly stochastic neighbor embedding on spheres. Pattern Recognition Letters, 128:100 106, 2019.
H. Lyu and Y. Li. Block majorization-miiiiizztiii with diminishing radius for constrained nonconvex optimization. 08 2023.
S. Marcel and Y. Rodriguez. Torchvision the machine-vision package of torch. In Proceedinggs of the 18th ACM international conference on Multimedia, pages 1485-1488, 2010.
H. Maron and Y. Lipman. (probably) concave graph matching. Advances in Neural Information Processing Systems, 31, 2018.
13

---

Published in Transactions on Machine Learning Research (06/2025)
L. McInnes, J. Healy, and J. Melville. Umap: Uniform manifold approximation and projection for dimension reduction. or Xio preprint orXiu: 10n.1802.00 03420, 2018.
F. Mémoli. Gromov- wasserstein distances and the metric approach to object matching. Foundations of computational mathematics, 11:417-487, 2011.
F. Mémoli and T. Needham. Gromov-monge quasi-metrics and distance distributions. arXi 2018, 2018.
F. Mémoli and T. Needham, Comparison results for gromov-wassersein and gromov-monge distances. arXiv preprint 2rx::2212. 14128, 2022.
R. Murray and A. Pickarski. On probabilistic embeddings in optimal dimension reduction. arXiv preprint arXiv:2488.02483, 2024.
Y. Nagano, Yamaguchi, Y. Fujita, and M. Koyama. wrapped normal distribution on hyperbolic space for gradient-based learning. In International Conference on Machine Learning, volume 97, pages 4693 -4702. PMLR, 09-15 Jun 2019.
S. A. Nene, S. K. Nayar, H. Murase, et al. Columbia object image library (coil-20). 1996.
M. Nickel and D. Kiela. Learning continuous hierarchies in the Lorentz model of hyperbolic geometry, In International Conference on Machine Learning, volume 80, pages 3779-3788, 10-15 Jul 2018.
A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and
A. Lerer. Automatic differentiation in pytorch. 2017.
F. Podregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, al. Scikit-learn: Machine learning in python. the Journal of machine Learning research, 12:2825-288 2830, 2011.
G. Peyré, M. Cuturi, and J. Solomon. Gromov-wasserstein averaging of kernel and distance matrices. In International conference on machine learning, pages 2664-2672. PMLR, 2016.
G. Peyré, M. Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trendso in Machine Learning, 11(5-6):355 607, 2019.
E. Quemener and M. Corvellec. Sidus-the solution for extreme doduplication of an operating system. Linar Jowrnal, 2013(335):3, 2013.
A. Rao, D. Barkley, G. S. França, and I. Yanai. Exploring tissue architecture using spatial transcriptomics. Nature, 596(7871):211-295)) 2021.
1. Redko, T. Vayer, R. Flamary, and N. Courty. Co-optimal transport. Advances in Neural Information Processing Systems, 177779),1 2020.
A. Rosenborg and Hirschberg. V-mcasure: A conditional entropy-based external cluster ovaluation measure. Proceedings of the 2007 joint conference on empirical methods in natural language processing and compulational natural language learning (FMNL.P-CNT.I), pages 410-420, 2007.
P. J. Rousseeuw. Silhouettes: graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics, 20:53-65, 1987.
S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. science, 290 (5500):2323 2326, 2000.
W. Saclens, R. Cannoodt, H. Todorov, and Y. Sacys. A comparison of single-cell trajectory inference methods Nature biolechnology, 37(5):5477555 554, 2019.
M. Scetbon, G. Peyré, and M. Cuturi. Linear-time gromov wasserstein distances using low rank couplings and costs. In International Conference Machine Learning, pages 19347-19382 19365. PMI.R, 2022.
14

---

Published in Transactions on Machine Learning Research (06/2025)
S. E. Schaeffer. Graph clustering. Computer science reciew, 1(1):27-64, 2007.
Schölkopf, A. Smola, and K.-.. Müller. Kernel principal component analysis. In International conference on artificial neural networks, pages 388-588. Springer, 1997.
R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, (1 -38, 1967.
J. Solomon, G. Peyré, V. G. Kim, and S. Sra. Entropic metric alignment for correspondence problems. ACM Transactions Graphics (ToG), 35(4):1-13, 2016.
K.-T. Sturm. The space of spaces: curvature bounds and gradient flows on the space of metric measure spaces. Xiee preprint or0i::1208.0 0434, 2012.
J. B. Tenenbaum, V. d. Silva, and J. C. Langford. A global geometric framework for nonlinear dimensionality reduction. science, (9055)))) (315000):311-2222 2323, 2000.
A. Thual, Q. H. TRAN, T. Zemskova, N. Courty, R. Flamary, S. Dehaene, and B. Thirion. Aligning individual brains with fused unbalanced gromov wasserstein. Advances in Newral Information Processing Systems, 35: 21792 21804, 2022.
P. Tseng. Convergence of block coordinate descent method for nondifferentiable minimization. Journal of optimization theory and applications, 109:475-44 494, 2001.
Van Assel, T. Espinasse, J. Chiquet, and F. Picard. probabilistic graph coupling view of dimension reduction. Advances in Neural Information Processing Systems, 35:10696-10708, 2022.
H. Van Assel, T. Vayer, R. Flamary, and N. Courty. SNEkhorn: Dimension reduction with symmetric entropic affinitics. In Thirth-seventh Conference on Neural Information Processing Systems, 2023.
L. Van der Maaten and G. Hinton. Visualizing data using t-snc. Journal of machine learning research, 9(11), 2008.
L. Van Der Maaten, E. Postma, J. Van den Herik, et al. Dimensionality reduction: a comparative. Mach Learn Res, 10(66-71), 2009.
T. Vayer, L. Chapel, R. Flamary, R. Tavenard, and N. Courty. Optimal transport for structured data with application on graphs. arXiv preprint arXiv:1805.091 091144 2018
C. Villani et al. Optimal transport: old and new, volume 338. Springer, 2009.
C. Vincent-Cuaz. Optimal transport for graph representation leurning. PhD thesis, Université Côte d' A'Auur, 2023.
C. Vincent-Cuaz, T. Vayer, R. Flamary, M. Corneli, and N. Courty. Online graph dictionary learning, In International conference on machine learning, pages 10564-10574. PMLR, 2021.
C. Vincent- Vincent-Cuaz, R. Flamary, M. Corneli, T. Vayer, and N. Courty, Semi-relaxed gromov wassersten divergence and applications graphs. In International Conference on Learning Representations, 2022a.
C. Vincent-Cuaz, R. Flamary, M. Corneli, T. Vayer, and N. Courty. Template basod graph neural network with optimal transport distances. Advonces in Neural Information Processing Systems, 35:11800 11814, 2022b.
M. Vladymyrov and M. Carreira- Perpinaa.. Entropic affinities: Properties and efficient numerical computation. In International conference on machine learning, pages 477-485. PMTRR 2013.
U. Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17:395-414 2007.
F. A. Wolf, Angerer, and F. J. Theis. Scanpy: large-scale single-cell gene expression data analysis. Genome biology, 19:1 5, 2018.
15

---

Published in Transactions on Machine Learning Research (06/2025)
H. Xiao, K. Rasul, and R. Vollgraf. Fashion-mnist: novel image dataset for benchmarking machine learning algorithms. Xiv preprint arXio:1008.077 07747, 2017.
Xu. Gromov Gromov-wasserster wasserstein factorization models for graph chustering. In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 6485, 2020.
H. Xu, D. Luo, and L. Carin. Scalable soanoovwwssersteed learning for graph partitioning and matching. Advances in neural information processing systems, 32, 2019.
A. Zeisel, A. B. Munoz-Manchadoo S. Codeluppi, Lönnerberg, G. La Manno, A. Juróus, Marques, H. Munguba, L. He, C. Betsholtz, et al. Cell types the mouse cortex and hippocampus revealed by single-cell rna-seq. Science, 347((222))11135 1142, 2015.
Z. Zeng, R. Zhu, Y. Xia, H. Zeng, and H. Tong. Generative graph dictionary learning. In A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volune 202 of Proceediiggs of Machine Learning Research, pages 40749 40769. PMLR, 23-29 Jul 2023.
A Appendix
Outline of the appendix:
Appendix provide the proofs for the results stated Scction 3, namely the one for Lemma 3.1
Appendix B. for Theorem 3.2 in Appendix [B.2 and additional necessary and sufficient conditions
under different assumptions as discussed in Remark 3.3 developed in Appendix B.3
Appendix C provide the definition of weak isomorphism in the GW framework, proofs regarding the characterization of the generalized srGW discrepancy as a divergence, mentioned in Section of the main paper.
Appendix proof of Theorem the clustcring propertics of srGW barycenters.
Appendix Algorithmic details for DistR under generic and low-rank settings.
Appendix Several additional information and results regarding the experiments detailed Section of the main paper.
provides details on methods' implementation, validation of hyperparameters, datasets and metrics.
Appendix
Appendix
compares DistR with COOT clustering.
Appendix 3
Best trade-off between metrics using t-SNE with la symmetrization and UMAP.
Appendix
4 reports complete scores on all datasets.
Appendix
and study homogeneity silhouette and homogencity kmcans NMI scores for various numbers of protolypes.
Appendix study the sensitivity of benchmarked methods to the embedding dimension d using ral DR methods.
Appendix compares computation time across methods.
Appendix
110
detail our proofs of concepts with hyperbolic DR kernels.
B Proof of results in Section 3
B.1 Proof of lemma 3.1
We recall the result.
16

---

Published in Transactions on Machine Learning Research (06/2025)
Lemma 3.1. Let Cz be permatation equivariant function and L any loss. The minimam of cq. equal to
DR is
min min ZER TENN
(2)
And conversely any (Z, sol. of Also, any sol. Z of eq- (DR) is such that (Z, id) is solution of eq- 2
eq. is such that Z is a solution of eq. (R)) up to
Proof. By suboptimality of a = id we clearly have
min min ZERNIO JESN
(C2 [C([)(()),))122)))))))))
S
min ZERNxd
(4)
For the other direction, take an optimal solution (Z, of eq. (Cz(Z))o())))))) [ZCx(Z)PT). (Cz(PZ)) for some permutation matrix P. But Pz is admissible for problem eq.
2 Using the permutation equivariance of Cz,
DR Hence
Σ (C22)))))))))))))) 2 ZERaxe dESx
min min
min
Σ [x(( (C2())
(5)
ZERN
B.2 Proof of theorem 3.2
In the following DS is the space of N N doubly stochastic matrices. We begin by proving the first point of We will rely on the simple, but useful, result below.
theorem 3.2
Proposition B.1. Let Q C R and Im(Cx) S NNN,, Suppose that L(a.-) is convex for any a E Q and
min ZERNXA
min ZERNECTEDS
Σ [(( (T22(
(6)
Then the minimam
DR is equal to minz (Cx)C Cz 1 IN, N 1x).
Proof. Consider any doubly stochastic matrix T and note that (TCz(Z)TT) lij = [x2(((())))TTT))) Using the convexity of L(a,.) for any a E Z and Jensen's 's inequality we have
Σ (((xx(X))s (T2 1s)) y
_ ((CT(T))TTTT)))))
RN
(7)
Σ ijki
L(Cx(())s, (Z)
In particular
min min Z TEDS
([(C 14) I min
min
Σ
(8)
2
TEDS ijkl
Hence, using eq. G
min z
Σ [(((x(X)) (C2(Z)))) min min TEDS ijkt
(9)
But the converse inequality is also true by sub-optimality of T In for the problem minz minTeDs Eijut L([Cx(X));; (Z) Overall
Now we conclude by using that the RHS of this equation is equivalent the minimization in of (W)(Cx(X), CZ))), I N IN) (both problems only differ from constant scaling factor N2). L (10)
min min min Σ [ (Z)(Z)lg) Z
(X(Cx(()) CC (222))
Z
TEDS 8ykl
2
17

---

Published in Transactions on Machine Learning Research (06/2025)
As consequence we have the following result. Proposition B.2. Let Ω R and Im(Cx) DN The minimum eq. DR is
equal to
minz GW (CW)(C(X), N 1N, N 1N) when:
(1) L(a,)) conver for any a E Ω and the image of Cz stoble by doubly-stochaste stochastic matrices, i.e.
ZZ R Nxd VT DS, R Nxd TC>
(11)
(ti) () conver and non non-decreasing for any a E Qt and
Z R Nxd
VT DS R Nxd
Cz(Y) TC2(Z))T,
(12)
where is understood element-wise i.e. A < B
v(i.), Ag Bag.
Proof. For the first point it suffices to that the condition eq. RNNd T E DS} S (Z) Z E Nxd (in fact we have equality by choosing T = IN) and thus eq. 6 holds and we apply proposition B.1
11 implies that (TC2(Z)TT Z
For the second point we will also show that eq. 6 holds. Consider Z T minimizers
of minZ,TeDS Eg [(((x(X)))s (TC2 Is). By hypothesis there exists Y E Rxx such that v(i.)) INT (Z ly 2 [C2(Y))s Since ((((x(X)),)))) is non-decreasing
for any (1,5) then Eij ((Cx(X) (Z) Ju] 2 Eij M[xx(((())) (Cz(()) and thus
Eg I([Cx(X)]a (Z'Z(Z "ly) 2 minz Ei [[Cx(X)] (Cz(Zilu) which gives the condition eq. 6
and we have the conclusion by proposition B.1
We recall that function R Rxd - R is called permutation invariant if R(PZ) R(Z) for any Z E Nxd and N x N permutation matrix P. From the previous results we have the following corollary, which proves, in particular, the first point of theorem 3.2
Corollary B.3. We have the following equivalences:
(i) Consider the spectral methods which correspond to Cz(Z)= ZZT and L T.z. Then for any Cx
Σ L2((Cx(X)]s (2;, ZERNX4 xd
min
(13)
and
min (Cx((C ZZ ZERaxd
N 1N. 1N)
(14)
are equal.
(i) Consider the cross-entropy loss L y) x log(x/v) and Cx such that Im(Cx) C RNN Suppose that the simvllaril in the oulput space can be wrillen as
v(i.)) E [NT (C f(2; - 22))R(()).
(15)
for some. logarithmicclly concave function f Rd - R2 and normalizing factor R Rxd + which is both concer and permutation invoriant. Then,
R
min Σ (x), 84
(16)
ZER
and
(Cx) 1 xxd
1 IN, IN)
min
(17)
equal.
18

---

Published in Transactions on Machine Learning Research (06/2025)
Proof. For the first point show that the condition eq. 11 of proposition B.2 is satisfied. Indeed take any Z, T then TC2(2)T' TZZ T = (TZ) = Cz( TZZ).
For the second point if we consider L(a, b) = a x b then we use that the neighbor embedding problem eq. is equivalent to minZERNNa - [X)] (oC (Z)ly) minzeRRxa Eig I([Cx(X)]s; (Cz (C2 (Z)lg) where (z(Z))l g(z, -zj) + log(R(Z)) with g log of. Since f is logarithmically concave 9 is convex. Moreover we have that L(a,.) is convex (it is linear) and non-decreasing since R4 in this case (Cx(X) is non- negative). Also for any Z E RNxd and T E DS we have, using Jensen's 's inequality since T is doubly-stochastic,
16
TCz(Z) 10 = Σ g(zz 21 TT + Log
2 g
(Zk: Tx)T log(R(Z))
(18)
Zi TiT Ex(Tg) log log(((()))
Now we will prove that. log(R(Z)) 2 (TT(R. Using Birkhofl's theorem [Birkhol] 1946) the matrix T can be decomposed a convex combination of permutation matrices, í.e. T Σx P where (Pε) are permutation matrices and A E R4 with Σx A = 1. Hence by convexity and Jensen's incquality R(TZ) ROEx Zz < Ex (). Now using that R is permutation invariant we get. R(P;Z) R(Z) and thus R(TZ)= & E RZ = R(Z). Since the logarithm is non-decreasing we have log( ([()) 2 log (T())
and, overall,
(C22 by 0 Σ 2x
Σ (TTT)) + log(R(TZ)) (C2(TZ))
(19)
Thus if we introduce Y. TZ we have (TC li fied. Thus we can apply proposition B.2 and state that minzER Enj [((Cx(X) (Z)())) and minzeRNXa (Cx Cz(Z), N + are equivalent which concludes as minzeRNxd (Cx)Cx(X) (Z). 1N, 1x) minzERNxe GW (Wx(Cx(X), C2(Z), 11N, 11N).
CZ(()) and eq. 12 satis-
It remains to prove the second point of theorem
3.2 as stated below.
Proposition B.4. Consider Im(Cx) S R CPD and for any Z
NxN
L LKL- Suppose that for any X the matrix Cx (X) is
Cz(Z) diag diag(Bz).
(20)
where az ,2 R. and Kz E R is such that log(Kz) is CPD. Then the minimum eq. DR is equal to minz (Cx Cz(Z), N IN, N IN).
Proof. To prove this result will show that, for any Z, the function
1 Te u( IN, N IN) (xx Cz(Z),T),
(21)
is actually concave. Indeed, in this case there exists minimizer which is an extremal point of u( IN, N 1x). By Birkhoff's theorem BBrkkooff 1946) these extreme points are the matrices J where P is N x N permutation matrix. Consequently, when the function eq. (21) is concave minimizing G (Z)), N INN, 1N) in Z is equivalent to minimizing in Z
min PERNXN permutation
Σ (C
min Σ [[(xx(X)) (Cz(Z) 107))))))))1111))))))))))))) (22) oESN
kyke
which is exactly the Gromov-Monge problem described in lemma 3.1 and thus the problem equivalent eq. DDD by lemma 3.1
19

---

Published in Transactions on Machine Learning Research (06/2025)
First note that L(x, y) zlog(z) xxlog()) y for any T E 1 1ln the loss
Er. xx,(Cx(X), Cz(Z), T) is equal to
Σ L[(xx(X))xa (C2 = ax + bz tjk!
Σicx(()) (X) ik log (z) T
kikl
.x + lz
Σ
C oz)( 8z W W5
kjkl
ax + bz
N
Σ
(23)
kjk
Cx (X) ik log ()g)(3)
N
Cx (X)] log( (inTa
sjkl
where ox. bz are terms that do not depend on T. Since the problem is quadratic the concavity only depends Tr(Cx(X)T)) log(Kz)T). From (Maron and Lipman
the term - Eijaal [Cx (X)] log( [ (a) 2018) we know that the function + - Tr (XyT" log(Kz)T) is concave on U( 1 CPD and log( Wz) is CPD. This concludes the proof.
N 1N, N xlN) when Cx (X
B.3 Necessary and sufficient condition
We give a necessary and sufficient condition under which the DR problem is equivalent to a GW problem Proposition B.5. Let C₁ E NNN, L , x R + R and C C R NNN a subspace of x N matrices. We suppose that is sluble by permulations i.e. C E C implies that PCP E for any N N permulation motriz P. Then
Σ min CLEC (1.))E[N
L([Cllg,
min C.ce
min TER Ti T 1w=1y
E((Cllu (C2la) TaTTa
(24)
kjkt
if and only if the optimal assignment problem
min 01.0
[an,
min C.ec
L(CC))s, (C2lo1(0o11))
(25)
admits on optimal soluiion (ai, a5) with of 02.
Proof. We note that the LHS of eq. (24 is always smaller than the RHS since T IN is admissible for the RHS problem. So both problems equal if and only if
min C.EC (LN))NNT
Σ
L(CCly ([C)l))(CClo)) S min C.EC
min
Σ A([Caly,(C2lkl TaT
(26)
TER
4ykt
Tls=1x
1x=1x
Now consider any C2 fixed and observe that
min TEDS
[T((Ch)s, TueTi 2 S/Rl
min T
[(CC))y, T
(27)
kykt
The latter problem is a co-optimal transport problem Rodko et al.) 2020 and, since it a bilincar problem, there is an optimal solution (T(), T(2) such that both T and T are an extremal point of DS which is the space of N N permutation matrices by Birkhoff's theorem (Birkhoff] 1946). This point was already
20

---

Published in Transactions on Machine Learning Research (06/2025)
noted by Konnoo (1976) but we recall the proof for completeness. We note Light L((Cilaj. (2) and consider an optimal solution (T!") T!?) of minrow ( ED ((T()) (TT)) = Eight Ts T Consider the linear problem minTeDs (TT, T Since it is a linear over the space of doubly stochastic matrices it admits a permutation matrix p()) as optimal solution. Also ((P(1) T S <(T.") T!?) by optimality. Now consider the lincar problem minreDs @(P(1), T), in the same it admits permutation matrix P(2) as optimal which implies that solution and by optimality (((()), ((2)) I o(P(") T!?) thus ((P(), z oTT!", T!?)
(P(1) ((2) an optimal solution. Combining with eq. (27 we get
min min C2EC
EL((C (())))
(C2)a) T 2 C2eC 0102ESN
min min
Σ [(Cilss (Cloo))))))))
TEDS ijhl
(28)
min o
0(01
Now suppose that the optimal assignment problem mine, f(01, 02) admits optimal solution (ai o with of oº. Then with eq. (31)
min min
MICilo (C2)a) TuIa 2 min Σ
L(Cilo (C2lo100r())
C TEDS
C,EC
ijki
(29)
min min
Σ
oESN CLee
if
Now since stable by permutation then ((()o(())n)) (eseetT= C C, Sx} = C and consequently mingess minc,ec Eij L((Cily) (C2)a())()) minc,ec Enj L(CClly) [C2]y). Consequently, using eq. 22
min CEec TEDS
min
Σ [1(Cc,s, Tx 2 min Σ I(Ci),, (Caly),
(30)
C2EC
15h1
and thus both are equal.
Conversely suppose that eq. 24 holds. Then, from eq. 31 have
min e1N
0()
min mi CgEc
min
I min Ciec
min TEDS
E (C2la) ([2]l) TaTT
(31)
ijke
= min
Σ
I(Ci), (Cali) = s(id, id),
C2EC
(i.))|N
which concludes the proof.
Remark B.6. The condition on the set of similarity matrices C is quite reasonable: it indicates that if C is an admissible similarity matrix, then permuting the rows and columns of C results in another admissible similarity matrix. For DR, the corresponding C is C = (Cz (Z) Rxd). In this case, C2 (Z) is of the form
(Cz(Z))s h(f(z, 2z), (Z))
(32)
where f 12d ddd R, R E and g R R which is permutation invariant (Bronstein ct al. 2021). then C is stable under permutation. Indeed, permuting the rows and columns of C2 (Z) by a is equivalent lo considering the similarity Cz(Y), where Y == -==0(1).'* " Moreover, most similarities in the target space considered in DR take the form eq. (32 (#(2,). ((2,)) )) (kernels such as in spectral methods with t id), f(z,, 2///E //m f(z,, 2m) (normalized similaritios such as in SNE and t-SNE). Also note that the condition on (Cz(Z) Z R Nxd of proposition B.5 is met soon as Cz R Nxd permutation equivariant.
R N is
21

---

Published in Transactions on Machine Learning Research (06/2025)
C Generalized Semi-relaxed Gromov-Wassersten is a divergence
Remark C.1 (Weak isomorphism). According the notion of weak isomorphism in Chowdhury and Mémoli (2019), for a graph (C, with corresponding discrete measure μ ". h2ō two nodes æ4 and æj are the "same" if they have the same internal perception i. Cs Cj Ci Cj and external perception wk (i,j), Ca Cj., Cla Ckj- So two graphs (C1, h1) and (C2, h2) are said to be weakly isomorphic, if there exist canonical representation (C,, he) such that card(supp(h.)) - S n, m and Mi E {0, 1) (resp. M2 {0, x)) such that for le E {1,2]
C. -MM.M, C.M. and he M. h4
(33)
We first emphasize a simple result extending a proof in Vincent- Cuan et al. (002ahh.
Proposition C.2. Let any divergence. I: f x Ω + Rt for Ω S R. then for any (C, h) and (C, (C,h), we have GWr(C, h) = 0 if and only if GWt,, (C, h) 0.
Proof. If GW, h, 0, then there exists T E (e((,, h) such that
Ex(
Σ ((Ci,, ijkl
(34)
so whenever T3 + O, we must have [[Cij, Ca) = 0 i.e C.j = Cut L is a divergence. Which implics that (C, C, for any other divergence L' well defined on any domain Ω Q necessarily including L2. D
Lemma C.3. Let any divergence 12 Ω + R+ for S R. Let (C,h) E S2ax x En and (C, h E mm x S-m- Then srGW, CC, = if and only if there erists h E Em such that (C, h) and (C, h are weakly isomorphic.
Proof. (=) As stG = 0 there exists E U(h h) such that (C, C.T. = hence GW (C, C, h. h) = 0. Using proposition C.2 GWL, L2 = 0 hence using Theorem 18 in Chowdhury and Mémoli (2019), implies that (C, and (C, are weakly isomorphic.
(+) As mentioned, (C,h) and (C,h) h) being weakly isomorphic implies that GW f, 0. So there exists T e U(h h), such that Ey (C, C, T) 0. Moreover T is admissible for the srGW problem as T E U(h h) c U, (h), thus (CGW h, h) = 0.
C.1 About trivial solutions of semi-relaxed GW when is not a proper divergence
We briefly describe here some trivial solutions of srGW,, when not a proper divergence. We recall that
C,,
min TER TI,
(,,, Σ LC, UGC.CT)TTTT
(35)
sjkl
Suppose that h E Mi and that C has a minimum value on its diagonal i.e. min(s.)) Cy = ming Tn. Suppose also that Va, ((a,)) is both convex and non-decreasing. First have E,
Ty 1 for any i [N]. Hence
hr
22

---

Published in Transactions on Machine Learning Research (06/2025)
using the convexity of L, Jensen inequality and the fact that L(a,)) is non-decreasing for any
Tik Tj ku h4 h
Σ
LCC
Gj
ijkt
Tik Tst hhahi
I(Cij,
Cki
hi h
Tak T3t hhg hs hj
Σ I(Cijo Σ Ra
(min C Crm) mmm
(36)
1/
Σ L(C (min C )mm))),hh
Σ L(Cij, (min Cnn))
Now suppose without loss of generality that ming Cs C11 then this gives
min
Σ
Cm),,,C)TTTTT 2
(37)
TER
TI, =h
Kjkl
h1
h2
0
Now consider the coupling T
It is admissible and satisfies
hN
0
C, Σ L(Ci),C1))le S min TER Th.
C,(C,,,T).
(38)
Consequently the coupling T* is optimal. However the solution given by this coupling is trivial: it consists in sending all the mass to one unique point. In another words, all the nodes in the input graph are sent to a unique node in the target graph. Note that this phenomena is impossible for standard GW because of the coupling constraints.
We emphasize that this hypothesis on cannot be satisfied as soon as L is a proper divergence. Indeed when L is a divergence the constraint " -L, is non-decreasing for any a" is not possible as would break the b (at some point L. must be decreasing).
divergence constraints Va, b I(a,b) b) and I(a,
Clustering properties: Proof of theorem 4.1
We recall that a matrix C E NNN is conditionally positive definite (CPD), resp. negative definite (CND), if Vx E RN, x In = s.t. x Tx 2 0, resp. I 0. We also consider the Hadamard product matrices as A AOB= (Ar x Byj)y. The i-th column a matrix T is the vector denoted by T. For a vector x E R we denote by diag(x) the diagonal
matrix whose clements are the I;.
We state below the theorem that we prove in this section.
Theorem 4.1. Let hx E En and L Lz- Suppose that for any X Nxp the matric Cx(X) CPD CND. Then the problem eq. Grr admits a membership matris as optimal coupling, i.e. there is minimizer of T Uuxxx) min (xx C, T) with only one non-zero value per row.
In order to prove this result we introduce the space of semi-relaxed couplings whose columns are not zero
u (bx) = (T E (bx): : Vi E [n], T. + 0},
(39)
and we will use the following lemma.
23

---

Published in Transactions on Machine Learning Research (06/2025)
Lemma D.1. Lethy EN, L L2 and Cx symmetric. For any T E (xx),, the matrir C(T) R"*" defoned by T Cx 1w) for (i,j) such that T.s and T. + olherwise
C(T)
(40)
0
is a minimizer of G C Rex Ex(Cx(XX, For T E Un the expression of the minimum
= cle Ix IN) (C(T) o C(T)))
(N)
(41)
which defines a continuous function on (hx).. If T E U, (hx) it becomes
G(T) ctc- Σ Cx 5)2
(T) IN) ctc cte-l diag ( T 1N) Cx 1 1
(42)
Proof. First see that C(T) is well defined since T&f #0 = T + because T is non-negative.
Consider, for T (hx), the function
F(T,C) CTT
= Σ [(Cx((X)]ix -x C.s))
(43)
sjkl
A development yields (using T T'NN- = hxx
F(T, = (x + CCED Tel) - 2 ECs(Cx(X)))TTTTT
(44)
ak
31
We can rewrite Eijre C (X) = T(TTCx(X)TT) (TT 1NT(TTNT)) 1n "ha Thus
Also we have (b) TEa)
Σ R
Σ
Tb)( Σ
Thd TT(TT [C
(45)
Overall
F(T, = cte (r Lw) (CC)) "
(46)
Now taking the derivative with respect to C, the first order conditions
0₂ F2 C) 2 (T (N)(TT IN) TTTCx(()T)) 0.
(47)
For (i.g) such that T. and c. + 0 we have [(T ((T)), 0. For (i.j) such that T. or Tuy
have [ @ (T (T IN) 1g matrix C(T) satisfies the first order conditions. When L = Le the problem min cin (x), Eigte ([Cx (X)] 2 Tha convex in C. The first order conditions sufficient hence C(T) is minimizer.
and also [T' Cx(X)T), (X)T 0. In particular the
Also CxC((X]==CH) C(T) O (T IN)(T by definition of C(T) thus
T Tr([C(T): (T 1T) 1x) 1w) (CTT) (C(T) @ C(T))) = =Tr((T Ix) IN) (C(T)
(48)
Hence
F(T, C(T)) cte- - Tr ( (T T o C(T)))
(49)
Consequently for T E Un (bx) such that Vi E [n]] Tui # 0 we have
T Cx(TJ
F(T, C(T)) = cte
TT (T. IN)
(50)
= cte -I|
Cx IN) 1 T (X))T 11) -
24

---

Published in Transactions on Machine Learning Research (06/2025)
It just remains to demonstrate the contimuity of G. We consider for (x, E RR x R the function
[x (xCxx(Xyx
when x #0 and y + 0
g(x,y) = xxlillyl
(51)
otherwise
and we show that is continuous, For (x,y) + (0.0) this is clear. Now using that
0 y(x,)) Eig (x Cx (x) (, wz)² = || (X)]l ||
(52)
[E:;)(,; U1)
[E U1)
this shows limx-o y) = 0 = g(0, y) and limy-s y) = 0 = =g(x,0) Now for T E (hx) we have G(T) = ctc- Eig (T T.5) which defines a continuous function.
To prove the theorem we will first prove that the function G T CER (xx) C, T) is concave U, (hx) and by contimuity argument it will be concave Un (hx). The concavity will allow us to prove that the minimum of G is achieved in extreme point of Un (hx) which is membership matrix.
Proposition D.2. Let hx E Ey, L L2 and suppose that Cx(X) is CPD or CND. Then the function G:T millCeR.... (x(Cx(X), T) is concave on Un (hx). Consequently theorem 4.1 holds.
Proof. We recall that F(T, C) Ex(Cx(XX, T) and = F(T, C(T)). From lemma D.1 we know that C(T) is minimizer of C F(T, C) hence it satisfies the first order conditions 0T C(T)) 0. Every quantity is differentiable on U, (hx). Hence, taking the derivative of G and using the first order conditions
VG(T) FT, C(T)) PT, C(T)) FT, C(T)).
(53)
We will found the expression of this gradient. In the proof of lemma D.1 we have seen that
F(T,) = cte Tr((T 1x) (x) - 2Tr(T Cx(X)TC)
(54)
ctc Tr 1 2TT'' CX(XTTC). (XTT))
Using that the derivative of T + Tr(T ATB) is A TB + ATB and that (X(X) is symmetric we get
0 I Ix1 C) - (X)(T 2Cx(X) Finally, applying to the symmetric matrix C=C(T)
(55)
VG(T) = F -() 2 o O(T) (CX(X))
(56)
In what follows we define
D(T) (T 1N) R
(57)
when applicable. Using the expression of the gradient we will show that G is concave on Ut (hx) and will conclude by continuity argument Un (hx). Take (P, Q) Ut (hx) U, (hx) we will prove
G(P) - -(Q) (VG((Q) Q) 00.
(58)
From lemma D.1
we have the expression (since P E U (hx))
G(P) = cte - AP (x.
(59)
(same for G(Q)) and
((Q)) (
(60)
25

---

Published in Transactions on Machine Learning Research (06/2025)
We will now calculate (VG(Q), which involves Q(C() @ C(Q), P) and (Cx(()QC(Q)) P). For the first term we have
Q(C(Q) @ P)
Tr(P =T((P' - 62 p
= -(""L))(((()C()))"""
=Tr "11000/0550000000000000000000100
(61)
= XT0000/0000000000000000000000000
For the second term
(Cx(X)QC(Q))))
(62)
This gives
VV((Q), 2P) 210000110100000000X0000001111110
-ID(())*P -
(63)
22I =0000-=--00000000000000eeee-raa0
- DP)) 4 (1
From this equation we get directly that
(VG(Q), Q) -
(64)
and (VG(Q), P) -2
Hence
G(P) - C
(VG(Q), + (VG(Q))
64
-ID(P)) 1PP.
(65)
eg-
-ID(P) PTx -ttp:////D..----P..-...CCC.CCCC.
(D(Q)
Cx C2 (1
+ +2ID
We note = PDDPPP E N = QD(Q)QT E NNN, the previous calculus shows that
( (VG))))))(V P - (X=II(UX((X)UXxxX)) - (IVxCx(X)VXxxX))
(66)
26

---

Published in Transactions on Machine Learning Research (06/2025)
Now note that
1'X==P(P)P'IN = = diag((
IN) P
Plx hx
(67)
Since U is symmetric we also have UIN hx and similarly we have the same result for V. Overall
V IN = 1NiN and V1y = UIN. Since Cx(X) is CPD or CND we can apply lemma D.3 below
which proves that XXXUUXXXXUCXXX (X)) - -T(VCx(X)VxxXX
+ XTX(VCX(XUCCxXX)
= and
consequently that is concave on U (hx)- We now use the continuity G to prove that concave on U, (hx).
Take P E U (hx) and Q Un (hx) Ut (hx) iee there exists & [n] such that Q. = 0. Without loss
of generality suppose & = 1. Consider for E N' the matrix Q(m) (1 -. QQ.)). Then
((m) Q as m + +00. Also since Q (m) E U (hx) we have by concavity of
C((1-A)P A)P (Q(m)) 2 (1 - AJG(P) (G(((()).
(68)
for any X e [0, 1]. Taking the limit as m + co gives, by contimuity of G,
G((1 - 2(P1 + XQ) 2 (1
+
(69)
and hence G is concave on Un (hx). This proves theorem Indeed the minimization of T E Un (hx)
+
ming Ey (Cx (X), T) is minimization of a concave function over & polytope, hence admits an extremity
of Un(hx) minimizer. But these extremities are membership matrices they be described {diag(hx)) )P: PE E {0, 1) P 1n 1N3 (Cao et al. 2022).
Lemma D.3. Lel C E NRNN be CPD or CND matrix. Then for any (P, (P,Q)RRAA x RNNN such that 1_1 = Q11
1y and PIN = Q1N we have
Tr(P (PTCQC) S CTCP ( + TTQ CQC)).
(70)
Proof. First, since C is symmetric,
PTCQC Q Tr (P - Q) T(P C(P-Q)C) -P
C"CC+
(71)
- TTrPP TrQQ - CQCC)
We note U =P - Q. Since P Q11
IN we have U = 0. In the same way U1 = 0. We introduce H IN $ 1N Iv1 the centering matrix. Note that
HUH (U - UH U 1 (UIN)1] N
(72)
N
Also C is CPD if and only if HCH is positive semi-definite (PSD). Indeed if HCH is PSD then take such
that x 1x 0. We then have Hx x and thus x Cx = x (HCH)x > 0. On the other hand when C
is CPD then take any x and see that x HCHx = (Hx) T C(Hx). But (Hx) 1x = x (HTIN) = 0. So
(Hx) C(Hx) 2 0.
By hypothesis C is CPD so HCH is PSD and symmetric, so has a square root. But using eq.
72 we get
Tr ((P Q) C(P - Q)C) I((UTUU)) = TH(HU T(HHT"HCHUHC)) HCHUHC)
= Tr(U (HCH)U(HCH) = [[HCH) U(HCH)+ 20
(73)
For the CND case is suffices to use that C is CND if and only if -C is CPD and that
Tr (PP - C(P- Tr (IP (-C) Q)(-)) which concludes the proof.
E Algorithmic details
We detail in the following the algorithms mentioned Section | to address the semi-relaxed GW divergence computation in our Block Coordinate Descent algorithm for the DistR problem. We begin with details on the computation of an equivalent objective function and its gradient, potentially under low-rank assumptions structures C and C.
27

---

Published in Transactions on Machine Learning Research (06/2025)
E.1 Objective function and gradient computation.
Problem statement. Let consider any matrices C' R C Rmx= and a probability vector h Σn. In all our use cases, we considered inner losses L. R R - R+ which can be decomposed following Proposition in Peyré et al. (2016). Namely we assumee the existence of functions fI, f2, hi, h2 such that
va Ω2
L(a.b) fi(a) + f2(b) (1(2)22(b)
(74)
More specifically we considered
I2(a, b) (a b)² = fi a², f2(0) b² h1(a) = a, h₂(b) 2b,
[xL((,,)) b) = a log a +
= SI (a) aloge a a, f2(b) b, h1(a) a, l22(b) log b
- a
LBCE(a, b) = a log (1 a) log
(a) a a + (1 - a) log(1 1 a), f2(b) = loog(1 b).
b
h1(a) a, h2(b) log 1
(L2)
In this setting, proposed to solve for the equivalent problem to srGW,,
min F(T) (EUL.(h)
(srGW-2)
where the objective function reads
F(T) =
(75)
= (LLIITTAAA))TT--TCTT (C)
Problem SrGW-2 usually a non-convex QP with Hessian # (2C) & 11 h2(CC) Øx hC In all cases this equivalent form is interesting as avoids computing the constant term (Ci(C), hh T) that requires O(N2) operations in all cases.
The gradient of F w.r.t T then reads
=) (RC - C
(76)
When Cx and Cz(Z) are symmetric, which the case all our experiments, this gradient reduces VTF 2 F2).
Low-rank factorization. Inspired from the work of Scetbon et al. 2022 propose implementations of srGW that leverage the low-rank nature of (X) and C2 (Z). Let assume that both structures can be exactly decomposed as follows, Cx(X) = A where A1, E N and Cz = B, with B1, Bi E g=xx, such that r <<< N and 8 << can differ respectively from respective dimensions p and d (e.g. for used squared Euclidean distance matrices r = p + 2 and s = d + 2 For both inner losses L we make use of the following factorization:
L- Computing the first term Fi coming for the optimized second marginal can benefit from being factored if &f << n. Indeed, as fz(Cz(Z) Cy 2 (B, O (B, one use the flattened out product operator described in Scetbon et al. 2022 Section 5), to compute Cz(Z) 1N in O(min(n² ns?)). This way Fi (T) results from stacking N times x in O(1) operations for total number of computions of N + O(min(n² "(min("",n=))) (nss)).)nddd its scalar product with T to compute the loss comes down to O(N)) additional operations.
Then computing (T) and its scalar product with T can be done following the development of Scetbon et al. (2022. Section 3) for the corresponding GW problem, in O(Nn(r + s) + rs(N + n)) operations. So the overall complexity at is O(Nn(r + 8) rs(N + n) + min(n?, ns²)).
28

---

Published in Transactions on Machine Learning Research (06/2025)
Lkz: In this setting (z(Cz(Z)) = Cz(ZZ) and (1(Cx( Cx(X) naturally preserves the low-rank nature of input matrices, but (Z)) = log Cz (Z)) does not. So computing the first term Fi, can be performed following this pesanthesis order 1w (T]]))) torm F2 should be computed following this order A1 ( T) their respective scalar product can be computed in O(Nn). So the overall complexity O(Nur + rPr). Similar considerations can be applied to LBCE-
in (n s)) operations. While the second
log (2 (Z))) in O(Nnr + rn²) operations. While
Notice that in the gaussian kernel case for neighbor embedding methods, where [Cz(Z)]ij exp(-1Iz: 2y up to some normalization. We have [ -I|z, z,122 which admits low-rank factorization such that we can recover the complexity illustrated above for I= -.2.
E.2 Solvers.
We develop next our extension of both the Mirror Descent and Conditional Gradient solvers first introduced in Vincent-Cuaz et al. 2022a), for any inner loss L that dccomposes as in equation 74
Mirror Descent algorithm. This solver comes down to solve for the exact srGW problem using mirror-
descent scheme w.r.t the KL geometry. At cach iteration (t), the solver comes down to, first computing
the gradient VT)(T()) given in equation 76 evaluated in T(), then updating the transport plan using the
following closed-fore solution to KL projection:
((++1) & diag K0
K
(77)
where K(0 exp F
(T)) and & > 0 is an hyperparameter tune. Proposition and Lemma in Vincentt Cuaz (2023 Chapter 6) provides that the Mirror-Descent algorithm converges to a stationary point non-asy when L= L2. A quick inspection of the proof suffices to see that this convergence holds for any losses L satisfying equation 74 up to adaptation of constants involved in the
Lemma.
Conditional Gradient algorithm. This algorithm, known to converge to local optimum Lacoste-Julien
2016 iterates over the 3 steps summarized in Algorithm I
Algorithm CG solver for srGW,
repeat
F() Compute gradient w.r.t T of equation 76
X
X20 4: ((+1) +- r +7*X(0 with 7"' E [0, 1] from exact-line search.
until convergence.
This algorithm consists solving at each iteration (i) a lincarization (X, of the problem equation Sr where ((T()) the gradient of the objective in equation 76 The solution of the linearized problem provides a descent direction ( T(0, and line-scarch whose optimal step can be found in closed form to update
the current solution TÚ0, We detail in the following this line-search step for any loss that can be decomposed
as in equation 74 It comes down for any T U4 (h), to solve the following problem:
arg min g(7) F(T (X T))
(78)
>E(0.1]
Observe that this objective function be developed as second order polynom g(7) on² + by c. To find
optimal 7 it suffices to express coefficients and b to conclude using Algorithm 2 in Vayer et al. (0018)
Denoting x XT,, 1, qx and T 1, = qr and following equation 75 we have
a = (1n(qx- - - - (C) (C) -
(79)
- = -F((() Fx()) + (T),
29

---

Published in Transactions on Machine Learning Research (06/2025)
Finally the coefficient b of the linear term is
= (1
- -),T)
(80)
F Appendix of experimental section
We report in the following subsections of this section:
implementation details, validation of hyperparameters, datasets and metrics.
2
comparison with COOT chustering.
3
Best. trade-off between metrics using t-SNE with l2 symmelrization and UMAP.
.
complete scores on all datasets for all kernels relating to PCA, t-SNE and UMAP.
study homogeneity vs silhouette score for various numbers of prototypes.
F.6
study homogeneity k-means score for various numbers of protolypes.
F.7 7
Benchmark between Spectral and Kmeans clustering
F88
study sensitivity w.r.t the embedding dimension d on spectral methods.
F.
computation time study.
F.10
Proofs of concepts with hyperbolic DR kernels.
F.1 Experimental setting
Sequential methods. We detail in the following the sequential methods DR->C and C-DR considered in our benchmark. DR-+C representations are constructed by first runming the DR method (Section 2.1 associated with Cx, Cz) thus obtaining intermediate representation Ž N Then, spectral clustering (Von Luxburg 2007) on the similarity matrix Cz (Z) performed to compute a cluster assignment matrix T E RVxx The final reduced representation in gnxd is the average of each point per cluster, t.e. the collection of the centroids, which is formally diag(h) - T Z [gad where h = T 1N. For C-DR, a cluster assignment. matrix T E NNxx is first computed using spectral clustering on (X)... Then, the cluster centroid diag(h)-17) X, where h T 1N, 1NN, is passed input to the DR method associated with
(L, Cx. Cz).
Implementation. Throughout, the spectral clustering implementation of scikit-learn (Pedregosa et al., 2011 used to perform either the clustering steps or the initialization of transport plans. For all methods, Z initialized from í.idd sampling of the standard Gaussian distribution N 1) and further optimized using PyTorch's automatic differentiation (Paszke al., 2017) with Adam optimizer (Kingma and Ba 2014). OT-based solvers are built upon the POT (Flamary ct al. 2021 library. k-means is performed using the
scikit-learn (Pedregosa et al. 2011 implementation.
Validated hyperparameters. For the SEA and UMAP based similarities, we validated perplexity across the set {20, 50, 100, 150, 200, 250). For UMAP, we further optimized before learning prototypes, the cocfficients a and involved in the parameterized Student kernel for Cz using the algorithm provided by Mclnnes el al. (2018). For all kernels, the number of output samples n spans set of 10 values, starting at the number of classes in the data and incrementing in steps of 20. For the computation of T in DistR (see Section 4.1] we benchmark our Conditional Gradient solver, and the Mirror Descent algorithm whose hyperparameter is validated in the two first values within the set {10')__-3 leading to stable optimization.
Datasets. We provide details about the datasets used in our study. For image datasets, we use COIL-20 Nene et al. 1996), MNIST and fashion-MNIST Xiao et al. 2017 Regarding single-cell genomics datasets,
hattps://wwe www.w...wiaa...ooooooo..o......a htt.a///aaeaaeetttaaaatattttt//a
Laken from Torchvision (Marcel Rodriguez 2010)
30

---

Published in Transactions on Machine Learning Research (06/2025)
rely PBMC 3 Wolf et al. 2018), SNAREsec chromatin and gene expression (Chen al. 2019) and the acRNA-seg datase from Zeisel et al. 2015 with two hierarchical levels of label. Dimensions are provided in Table ] When the dimensionality of dataset exceeds 50, we pre-process it by applying a PCA in dimension 50, as done in practice
(Van der Maaten and Hinton 2008
Table 1: Dataset Sizes.
Number of samples
Dimensionality
Number of classes
MNIST
10000
784
10
F- F-MNIST
10000
784
10
COIL
1440
16384
20
SNAREscq (chromatin)
1047
19
SNAREseq (gene expressioon)
10
Zeisel
3005
5000
(8, 49)
PBMC
2638
1838
Scores. The homogeneity and NMI scores are taken from Torchmetrics scores used in the experiments are computed follows.
Detlefsen al. 2022) The other
i) Silhouete:: The first step is to do weighted majority vote and associate a label Uh to each prototype Mi.- This label is defined by the y maximizing ΣiErn] y= where yi is the true class label of the input data point x,. Then, to incorporate the relative importance wy [hzl of each prototype zj, we define the weighted mean intra-cluster and nearest-clusted distances that read
wk e [n]] az Y, w) Ejelo in wyd(zj, 2x)
Σje[m) 1 - uxdddd( 2x)
and by (Z, Y. w) min Σje('] 1 -ún Wj
k'fk
Ejen =
(81)
such that the silhouctte coefficient 8 (Z, Y, w)
b and the final silhouctte score reads as
max(br,ax)
S Y, w) = Σ kE[m]
WkSk (Z, w)
(82)
d) k-means: We first run the k-means algorithm the prototypes {zx [e(r) giving the predicted label Te for each h E [n]. The predicted label for each input data point E [N] is then given by its prototype's predicted label i.e. D: Yarg max, Tuk Then we compute the NMI score (Kválseth 2017) between (D:)ie[N and the ground-truth class labels (μi)ae[N)
F.2 Comparison with COOT clustering
The CoOptimll-Tranppoooo (COOT) problem, proposed in Redko et al. 2020), is as follows,
min T. R.) Th) iyki
min
E(xi= Zz1)² T 1y T Ja
(COOT)
where h, E ΣN, h, E En, he E Ep and he E Eg. One can seek to optimize the above objective with respect to Z to obtain a competitor method to DistR. This problem is called COOT clustering in Redko el al. (2020). the latter, T, then plays the role soft clustering matrix of its columns. The above is thus a lincar DR model.
soft clustering matrix the rows of X while T, can be seen as a
In Table 2 we display the homogeneity values obtained with COOT along with the methods described Section Precisely, it measures to what extent the clustering given by T, groups points with the same ground truth label. One can notice that COOT falls short compared to its competitors that leverage affinity matrices as in state-of-the-art (non-linear) DR methods.
F.3 Best trade-off between metrics using t-SNE and UMAP
*downloaded from Scanpy WWlf al. 2018).
http:////stuuuu.ccom//eenaaaaaas pytthoo
31

---

Published in Transactions on Machine Learning Research (06/2025)
methhods
Cx / Cz
COIL
MINST
FMNIST
PBMC
SNA1
SNA2
DistR (ours)
SEA St.
99.50 (0.60)
97.80 (0.00)
93.40 (0.00)
97.10 (0.00)
100.00 (0.00)
100.00 (0.00)
DR -C
98.10 (0.50)
93.30 (3.40)
94.30 (2.80)
96.30 (0.90)
100.00 (0.00)
100.00 (0.00)
C-+DR.
100.00 (0.00)
97.80 (0.00)
93.40 (0.00)
97.10 (0.00)
100.00 (0.00)
100.00 (0.00)
COOT
NA
43.90 (0.70)
(3.00)
8.50 (0.90)
15.90 (1.90)
44.00 (4.60)
49.70 (8.60)
DistR (ours)
6 hhr (1)=
96.80 (0.70)
97.00 (0.90)
93.40 (0.00)
97.10 (0.40)
80.50 (0.00)
100.00 (0.00)
DR-+C
73.30 (1.80)
98.30 (3.40)
93.20 (1.90)
90.20 (1.40)
72.70 (6.00)
(20.00)
C-+DR,
83.70 (0.00)
100.00 (0.00)
93.40 (0.00)
93.30 (0.00)
80.50 (0.00)
100.00 (0.00)
COOT
NA
45.50 (1.60)
(2.10)
(2.80)
(2.40)
(5.90)
(16.60)
Table 2: Best homogeneity scores for validated in a span up to 200 with increments of 20.
/
Student
UMAP
EA Student
UMAP
norm. homageneity
norm. homogeneity
norm homogeneity
norm, homogeneity
02-
1o
24
25
norm. silhouette
norm. silhouette
norm. k-means NMI
norm. k-means NMI
Figure 6: Complements to Figure 3 of the main paper for t-SNE with 12 symmetrization and UMAP Best trade-off between homogeneity vs silhouette (2 first plots), and homogeneity vs NMI (2 last plots). Scores are normalized in [0. 1] via min-max scaling over a dataset. Small markers are scores per dataset for runs, while big ones are their mean. We illustrate the 20-80% percentiles of scores per method as a colored surface.
F.4 Complete scores
We complete the results shown Scction by providing the scores obtainod on all datascts and models. Scores are plotted in Figure 7 for the PCA model, in Figure g for the t-SNE model with entropic symmetrization, then in Figure g with la symmetrization, and Figure 10 for the UMAP model.
32

---

Published in Transactions on Machine Learning Research (06/2025)
Fishion-MM MNIST
SNAREag!
17.
10-
17.
Silbometer
10-
330
-17
DR-+0
DotR (oua))
SNARFane2
PBMC
Trioel
Zrieol 2
15-
-5-
Silboaete
-15
39
-54
COLL
MNIST
MaISI
SNAREseq1
34-
42
10-20001 NAI1
56-
30+
84
10
05
10
150
24.
36
SNARE
PBMC
Zeioel
Zeieol 2
60-
mocoaaa NMI
GL
18
19
590
170 MNIST
280
190
190
590
28
28
COLL
Faeliio--MNSTT
SNAREsa11
75-
300
201
159
200
380
150
3.
190
150
2
Zeieel
Honogenetty
86-
Numbor of prototypes
Number prototypes
prototypos
Number prosorypes
Figure Scores (x100) with respect to the number of prototypes (in I0) produced by DistR using the PCA model: (. ) similarity for Cx (Van Assel al. 2023), ()) for Cz and loss Lz.
33

---

Published in Transactions on Machine Learning Research (06/2025)
COIT.
Fushion-MNIST MNIST
SNAREay!
15-
Silomette
C-Di DEAR (ous))
18-
SN
PBMC
Zeisel
Zejeel
Siloocte 3
-12
COLL
MNIST
MaNSI
SNAREss]1
65.0
02.5
lecanome NA11
375
530
32.
200
50
100
300
150
24
150
30
PIMC
Zeioel
Zeisel 2
totouue NMI
810
72.5
77.
712
67.5
10
190
290
180
120
3/6
190
290
28
COLL
MNIST
Faeliio--MNIST
SNAREsey1
97.5
55-
[lmnogreett
9255
12
34.9
55-
42-
550
240
230
03
159
200
300
150
180
150
2N
Zeisel
Zejeol 2
Honegeetty
120
Number of prototypes
Numaber prototypes
Number of prototypes
Number of procolypes
Figure 8: Scores (x100) with respect to the number prototypes (in IR?) produced by DistR using the t-SNE model: SEA similarity for Cx (Van Assel et al. 2023), Student's kernel for Cz and loss T.L.-
34

---

Published in Transactions on Machine Learning Research (06/2025)
Trahio-MNNIST
SNAREa0g!
Salbowett
13-
Disk lour))
5-
DR--CC
C 1DR
0-
-30
SNAREseq2
PBMC
Zejeel
Zeleel
37-
53x
15-
23-
COLL
MNIST
Paslio-MNNIST
SNAREsag1
65.0
625-
Secneenn NM
34-
23.5
32-
942
SNAREsq22
PBMC
Zrieel
Zeteol 2
61-
k-mnemn NMI
12
725
32-
335+
32-
22
675
32-
41-
169
COIL
MNIST
FeshionnMM MNIST
SNAREEsq1
101-
54-
Horaggrtnity
66-
K7.5
K53-
85
525
02
SNAREseg2
PBMC
Zeieel
Zeieel 2
63-
Hocrogrncity
112-
26
100
Number prolotypes
Number of peotolypes
Numbet of protodypes
Number of prototypes
Figure 9: Scores (x100) with respect to the number of prototypes (in R2) produced by DistR using the t-SNE model with le symmetrization: EA similarity for Cx (Van der Maaten and Hinton 2008), Student's kernel for Cz and loss LKL-
35

---

Published in Transactions on Machine Learning Research (06/2025)
COIT.
SNAREag!
20-
Siloometer
-00-
-30
-0
C-Di
20-
-12
-52
SNARFan:2
PBMC
Trioel
Zrieol 2
Silboaete
10-
-20
COLL
MNIST
MNIST
SNAREss]
Se-annon NAI1
2
43
30-
20-
10-
00
34
100
201
350
36
SNAREac;2
PBMC
Zeioe
Zeieol 2
00-
I: 11-ooaan NMI
40x
15
290
100
14
36
06
190
150
280
28
28
COLL
MNIST
Faelicn--MNIST
SNAREsey1
22
76-
330
201
30
3/
100
150
201
180
150
2NAREa<2
PBMC
Zeisel
Zejeol 2
180
Honogeneit
06-
Numbor of prototypes
Number prototyppe
Numbor prototypos
Number prosorypes
Figure 10: Scores (x100) with respect to the number of prototypes (in R2) produced by DistR using the UMAP model: t-conorm of the smoothed nearest neighbors Unj = - 1xxx,xx;, xx,12 1allz-puui pejos for Cx (McInnes al. 2018 E6.1 parameterized Student's kernel for Cz (MeInnes et al. 2018 Eq.17) and loss LBCE.
36

---

Published in Transactions on Machine Learning Research (06/2025)
F.5 Dynamics between homogeneity and silhouette scores across numbers of prototypes
For each method depending on some hyperparameters, models performing the best on average across all numbers of protolypes are illustrated in Figure for spectral methods and in Figure 12 neighbor embedding ones.
DistR (ours)
DR C
DR
COIL MNIST
Fashlon-MNIST PBMC
SNA1 SNA2
ZEISEL ZEISEL
)pe/ In-)
nelasses
Melasses
20
Netosses
40
norm, homageneity
Aicoosses 60
nelosses 80
Aelasses 100
60
norm, homageneity
4.4
Classes 120
Relasses +14
Cllasses 160
nom, homageneity
I
silhouette
silhouette
silhouette
Figure 11: Trade-off between homogeneity vs silhouctte scores using PCA model across various numbers of prototypes n. The illust follows the same principal than fig. 3
37

---

Published in Transactions on Machine Learning Research (06/2025)
DistR tours)
DR
DR
Fashion-MNIST PBMC
MNIST
ZEISEL
(SEA/ Student) Aclasses
nelesses + 20
Nelesses + 40
14.
norm, homogeneity
18-
26
42
044
Aelasses 60
Aleesses + 80
Aclssses
100
18-
+00
*
norm. homogeneity
0.3
12
0.4
82
6.0
0.1
Relasses
120
Reissses
140
neiasses
160
18.
norm, homogeneity
05
52
02
norm. silhouette
norm. silhouette
norm. silhouette
Figure 12: Trade-off between homogeneity vs silhouette scores using (-SNE model across various numbers of prototypes n. The illustration follows the same principal than fig.
F.6 Dynamics between homogeneity and NMI scores across numbers of prototypes
For each method depending on some hyperparameters, models performing the best on average across all numbers of prototypes are illustrated in Figure for spectral methods and in Figure for neighbor embedding ones.
38

---

Published in Transactions on Machine Learning Research (06/2025)
DistR tours)
DR
DR
Fashlon-MNIST PBMC
MNIST
ZEISEL
)xe Iee)
Aeisssss
nelosses + 20
Neiscces 40
+
norm, homageneity
0.
02
00
02
9.4
9.
60
0.2
0.
Aelosses
.00
Nelasses
80
Aeiasses
100
10
norm. homogeneity
52
0.
0/
00
06
0.0
0.1
nelasses
120
Reissses
140
Relasses
160
norm, homageneity
ne
15
16
0.2
en
64
k-means NMI norm.
norm. k-means NMI
norm. k-means NMI
Figure 13: Trade-off between homogeneity vs silhouette scores using PCA model across various numbers of prototypes n. The illustration follows the same principal than fig.
39

---

Published in Transactions on Machine Learning Research (06/2025)
DistR (ours)
DR
DR
Fashlon-MNIST PBMC
MNIST
SNA2
ZEISEL
(SEA/ Student)
Aclasses
olessses
+ 20
Nelasses + 40
1t.
**
norm, homogeneity
IE+
6.4-
8/3
0.
6.
Aclasses 60
clessses + 80
Ressssss 100
norm. homogeneity
30-0
1/5
15
03
544
3.0
0.4
0.5
Relasses 120
Reiasses
140
Delasses
160
16.
wxth
norm, homogeneity
02 en
62
36
norm. k-means NMI
norm. k-means NMI
nom. k-means NMI
Figure 14: Trade-off between homogeneity vs silhouette scores using (-SNE model across various numbers of prototypes n. The illustration follows the same principal than fig.
F.7 Benchmark between Spectral and Kmeans clustering
We compare in the following the performances of the joint DR/clustering algorithu across different clustering strategies. Namely for C + DR and DR - C, we compare our default choice of the Spectral Clustering (SC) algorithm with a Kmeans algorithm. Then for DistR, we compare the use of both algorithms to initialize the transport plans. Performances for both PCA and t-SNE kernels with KL projection are reported in Figure 15 On average across datasets, we observe that DistR (SC), the default for the method, achives better trade-off w.r.t the silhouette and Kmeans N.M scores than DistR (Kmeans), while maintaining similar homogencity scores. A similar behavior is observed for C-DR using PCA like kernels, whereas overall higher performances are observed for C-DD (SC) than C-DR (Kmeans). Therefore SC appears as a better clustering strategies to maximize performances for both DistR and C-DR, where DistR (SC) remains the most competitive methods. Finally, DR- (Kmeans) consistently leads higher homogeneity scores with comparable silhouette scores than DR-+C (SC), but their respective trade-off between homogeneity and Kmeans NMT are more erratic.
40

---

Published in Transactions on Machine Learning Research (06/2025)
DistR (SC) DR = (SC) DR (SC)
COIL MNIST
Fashion-MNIST
SNAT * ZEISEL
mean
DistR (Kmeans)
(Kmeans)
(Kmeans)
PBMC
ZEISEL
( 1a= ( )n=
( )R2 / ( lna
14.
norm. homogeneity
norm. homogeneity
6.5-
6.-
6.2-
60
ck
0/3
3/
2/
1.4
norm. silhouette
norm. k-means
SEA/ / Student
SEA / Student
144
10-
norm. homogeneity
homogeneity
X
X
X
60- 6
2.0
norm. silhouette
norm. k-means NMI
Figure 15: Best trade-off between homogeneity vs silhouette (left), and homogeneity vs NMI (right). Scores are normalized [0, 1] via min-max scaling over dataset. Small markers represent scores averaged over runs for given dataset, while big ones their mean over datasets. The illustration follows the same principal than fig. g
41

---

Published in Transactions on Machine Learning Research (06/2025)
F.8
Sensitivity analysis w.r.t the embedding dimension
ju:20,
CHARissel
FRME (ru8,
10
COLL 1w200. di 10
StAntscot lind, 2
HBME ludd 2w
con mnwe0, fw21
SwaRtess the 25
P096: (1ed8 P-
V
197-
182
150
COLL 20-20, 5)
SVARtsegl in:4, 1
TBMC h- 2=
SHARtae2 1=2)
PBNC (1443
120 $w5)
5
7BHC 21
100
COBLES 100 p
SMiRErng 13:166 20.
M
pw
MA4Bo pw
(re1s8
briocding 9 diuusionn
Enbolsing Soo
bebiiiing dincsson
Entedingg denemion
Emboliing dkmateen
Brbesding demenaon
Figure 16: Sensitivity analysis w.r.t the embedding dimension d for the spectral method analog to PCA with (x(X) = XX and Cz (Z) ZZ T. Over three datasets, we make vary the embedding dimension from d {12, 14, ...,, 20).
F.9 Computation time comparison
We compare in the following the computation time for all methods benchmarked in Table when using a spectral method and a SNE method. For fair comparison across methods, we do not apply here the low-rank factorizations used for experiments of the main paper. All experiments were done on server using GPU (Tesla VI00-SXMM2-2GIB) and composed of 18 cores Intel(R) Xeon(R) Gold 6240 CPU @ 2.60GHz. All DR steps benefit from GPU compatible implementation, while spectral clustering (SC) was performed CPU using scikit-learn implementation running CPUs.
Notice that to run our experiments we precomputed and saved SC steps for the maximum number of prototypes ran on the input structure Cx (X), for all benchmarked methods e.g CDR used for clustering or DistR and COOT used for initialization the transport plans. As DRC performs SC over learned embeddings using (Z) cannot be precomputed, so we include the time of performing SC for all n for all methods achieve fair comparison. Results for three datasets of increasing sample sizes SNA1, ZEISEL and MNIST (See Appendix F. reported in Figure 18 with pre-computation and in Figure 17 without.
42

---

Published in Transactions on Machine Learning Research (06/2025)
G In / 6 hes runtimes with precomputing
SNA1
ZEISEL
MNIST
10a
10"
10a
DIstR,
runtimes
10A
CDR
DRC
10
10'
101
10a
10P
44
84
124
164
204
48
80
128
168
208
10
50
90
130
170
210
number of protolypes
number of prototypes
number prototypes
SEA / Student runtimes with precomputing
SNAI
ZEISEL
MNIST
DistR
107
runtimes (5)
DistR,
10
DRC
CDR
COOT
100
10"
101
10-3
164
48
50
210
number of prototypes
number of prototypes
number of prototypes
Figure 17: Computation time comparison depending on the number of prototypes n for all methods over datasets, for 5 different initializations, while procomputing spectral clustering on the input structure Cx (X).
Jas/ runtimes without precomputing
10a
DissR
Disin
102
2.02
CDR
DOOT
10a
100
10'
84
124
164
204
48
88
128
168
208
10
50
90
130
170
210
number prototypes
number of prototypes
number prototypes
SEA/ Student: runtimes without precomputing
DistR
DRC
CDR
10a
COOT
1.0
101
124
208
10
130
170
number prototypes
number of prototypes
number prototypes
Figure 18: Computation time comparison depending the number of prototypes for all methods over datasets, for 5 different initializations, without precomputing spectral clustering the input structure Cx(X).
Overall our DistR models are competitive with all benchmarked methods terms of computation time too and can be run (and further validated) few seconds using a GPU after precomputing the spectral embeddings for medium size dataset (N=10000 for MNIST), by efficienlly leveraging low-rank kernels often used in the DR literature.
43

---

Published in Transactions on Machine Learning Research (06/2025)
F.10 Proofs of concepts with hyperbolic kernels
Hyperbolic spaces (Chami el al. 2021 Fan el al. 2022: Guo el al. 2022 Lin et al. 2023) are of particular interest as they can capture hierarchical structures more effectively than Euclidean spaces and mitigate the curse of dimensionality by producing representations with lower distortion rates. For instance, Guo et al. (2022) adapted t-SNE by using the Poincaré distance and by changing the Student's t-distribution with a more general hyperbolic Cauchy distribution. Notions of projection subspaces can also be adapted, e.g. Chami ee al. (2021) use horospheres as one-dimensional subspaces. To match our experiments with the neighbor embeddings in Euclidean settings, adapt the Symmetric Entropic Affinity (SEA) from Van Assel et al. (2023) for Cx and the scalar-normalized student similarity for Cz (Van der Maaten and Hinton 2008 by simply changing the Euclidean distance by an hyperbolic distance.
Implementation details. Computations in Hyperbolic spaces are done with Geoopt (Kochurov et al. 2020 and the RAAam optimizer (Bécigneul and Ganca 2018) replaces Adam. A Wrappod Normal distribution in Hyperbolic spaces (Nagano et al. 2019 is used to initialize Z in the hyperbolic setting. All the computations were conducted in the Lorentz model (Nickel and Kiela 2018), which is less prone to numerical errors. We used the distance function from Nickel and Kiela 2018) to form Cz- Specifically, we employed scalar-normalized Gaussian kernel, where the distances between points are computed in Nickel and Kiela (2018 Equation 1). After optimization, results are projected back to the Poincaré ball for visualization purposes. In this hyperbolic context, we adopted the formulation of Guo et al. (2022) which generalizes Student's 's t-distribution by Hyperbolic Cauchy distributions (denoted as H-Student in the results). Notice that Guo et al. (2022) considered weighted sums of DR objective depending respectively on L2 and LKL, including 2 additional hyperparamaters, plus various validated curvature levels for the inner hyperbolic distances. In the following experiments, we only kept LKL for comparison with the Euclidean SNE-based methods illustrated in Section 5 and previous Sections of while validating the same hyperparameters and setting the space curvature to 1. The silhouette score was adapted to this kernel considering the Hyperbolic distance instead of the Faclidean one, and we implemented a Hyperbolic Kmeans whose barycenters are estimated using the RAdam optimizer to compute the NMI scores.
We first report in Figure 19 a relative comparison of the best trade-off between local and global metrics achieved by all methods. Similarly to visualizations of the main paper, we considered for each method and dataset, the model maximizing the sum of the two normalized metrics to account for their different ranges. DistR, being once again present on the top-right of all plots, provides average the most discriminant low-dimensional representations endowed with a simple geometry, seconded by C-
Results.
DGtN tours)
COL
SNA 2310L 1
MRIST
ISMC
5MA2
231eL. 2
SEA/ H-Suudent
SEA/
e-Sddeent
norm, hamogeneity
norm. homogeneityy
02-
02
norm. silhouette
norm. k-means NMI
Figure 19: Best trade-off between homogeneity vs silhonette (2 first plots), and homogeneity vs NMI (2 last plots). Scores are normalized in [o. 1] via min-max scaling over a dataset. Small markers represent scores for 5 runs for a given dataset, while big ones are their mean. For each method we illustrate the 20-80% percentiles of normalized scores as a colored surface.
Then we report in Figure20jaa abbsllue performances for all methods and all datasets across various n, as done Appendix F.4 We can observe that both DistR. and C-DR achieve fairly high NMI and homogeneity scores across all settings, while DistR performs significantly better on average across the tested number of prototypes
44

---

Published in Transactions on Machine Learning Research (06/2025)
However, DR->C struggles significantly to learn both globally and individually discriminant prototypes. Notice that DR-C's homogeneity scores are significantly lower average than both benchmarked Fuclidean kernels. This mitigates drastically the significance of the silhouette scores computed for this method, letting essentially DistR and C-DR to compare. Even though DistR outperforms consistently C-DR w.r.t silhouette scores, these scores remain significantly lower than with the other t-SNE kernels. As the latter is equivalent. a null curvature, this indicates that further fine-tuning of the curvature within these hyperbolic kernels could be beneficial. Nevertheless, these results confirm the versatility of our DistR approach, capable of operating non -uclldenn geometries.
COIL
MINST
Tieton-MNIST
SNARBagg!
Sihooett
C-Bi
SNARR
Sllorote
100
COLL
MNNIST
Fadlin-MNIST
SKARBeeg!
Emoas N011
Lascus N011
100
COIL
MINIST
Fellin-MNSST
SNARB5el
Hounggastt
10
Nutder pootolpps
Nindies
Nunles
Figure 20: Scores (x100) with respect to the number of prototypes (in R00) produced by DistR using the Hyperbolic Student model.
45


Reviews:
[{"review_id": "Fu7COZFHHp", "reviewer": "Reviewer_QBZR", "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": "Overall, the paper is well-written, with novel ideas, and what appears to be sound statements. \n\n**comments (mixed strengths/weaknesses):**\n- I like section 2.1. One thing I am missing is that in PCA, kernel PCA, etc., it is possible to reuse the mapping from high dimension to low dimension on new test points. However, in tSNE this is not possible AFAIK. What is the reason for this, and is this reason visible from the perspective of equations (DR) and (NE)?\n- A nice direction for future work is to cast the widely popular t-SNE in terms of GW, using a result similar to Theorem 3.2 (ii). The authors note that t-SNE doesn't currently fit, because the pairwise similarity matrix in original dataspace is not CPD. The authors did a good job in pointing out the current limitation of their work.\n- Dimension reduction (which optimises equally weighted dirac masses) can be generalised to distributional reduction, which additionally optimises the weights. \n- DistR is a nonconvex problem, and the authors attempt to find local minima using a block coordinate descent, where one of the coordinate block steps is simply a semi relaxed DistR, which is \"computationally simpler\" (by the way, could you please indicate more precisely how it is computationally simpler, and what can be achieved?). \n- The text before Theorem 4.1 seems to speak in very general terms, however the theorem requires $L=L_2$ (and more mildly, CPD/CND similarity in input space). Shouldn't you add some sentences to $L=L_2$ in the preceeding text?\n- There is no direct link between DistR and and clustering, the authors do a good job in providing a connection between srGWB and hard clustering, and therefore provide an intuitive motivation for using GWB for clustering. The limitation but applicability of this result is appropriately discussed.\n- Figure 3 takes a little while to understand, but is overall very clearly presented, with a good caption.\n- I find it curious that the bottom-right subfigure in figure 4 has a constant homogeneity of 100 using distR. Could you explain this? \n\n\n\nMinor:\n- Figure 1 caption \"two toys examples\" should be \"two toy examples\" \n- \"the findings of (Chen et al., 2023), which\" should be \"the findings of Chen et al. (2023), which \"", "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": "Gromov-Wasserstein is able to define distance or compare data distributions over different (dimensional) domains, making it a good candidate for a fundamental technique in dimension reduction and clustering. The authors use GW to bridge many examples of clustering and dimension reduction, allowing for a joint dimension reduction and clustering technique. The approach (called distributional reduction) is broadly applicable and supported by theoretical results which recover special cases and other connections to existing methods. Demonstrations are provided on multiple popular datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address my comments above. Namely,\n- Discuss the ability of various methods to perform embeddings on test data\n- Please clarify computational aspects of your approach. In what sense is srDistR simpler?\n- Clarify text surrounding Theorem 4.1\n- Could you comment on the bottom-right part of Figure 4?\n\nPlease address the minor points too (no need to respond here).", "broader_impact_concerns": "This work presents general methodology and unifying theory, without a specific application in mind. The usual issues surrounding use of technology apply. As such, a broader impact statement is not necessary.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ms09p7yLDx", "reviewer": "Reviewer_3D1w", "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": "Strengths:\n- For the most part the paper is well written and relatively easy to read.\n- Although the connections between optimal transport and the problems of clustering and dimension reduction are not new, the paper does a good job of concisely characterising the relationships.\n- Theorem 3.2 (although strangely worded; see below) has some nice generality, and if novel is a nice contribution in itself\n- Theorem 4.1 (although needing some more discussion; see below) is also a nice contribution.\n\nWeaknesses:\n- There are some grammatical errors/quirks, some of which I indicate below in the following section.\n- Unless I am missing something Lemma 3.1 is completely trivial.\n- The experiments are not well designed, since \"overclustering\" using spectral clustering is nonsensical. In particular, the eigenvectors of the Laplacian(s) which are not associated with actual clusters will very easily lead to points arbitrarily far from one another being grouped together.", "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 describes a unified framework for \"data summarisation\", within which both the problems of clustering and dimension reduction lie,\nin terms of optimal transport and the Goromov-Wasserstein distance. \nSome theoretical results are given which characterise existing models within this framework. In addition the authors show that a relaxation\nof the problem, which does not tie the graph of the reduction to an explicit form, has a solution which induces a hard clustering.\nAn alternating algorithm for jointly reducing the dimension and granularity (the latter potentially for clustering) of the data is described.\n\nResults from experiments comparing their approach with the alternatives of performing the dimension reduction and clustering tasks separately, in both orders (dimension reduction then clustering, and clustering then dimension reduction) show some promise of the proposed approach.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "-- within Theorem 3.2 the way $\\alpha_Z, \\beta_Z$ are introduced is strange. Perhaps replace the word \"where\" with \"for some\"\n-- within the paragraph \"DR as GW Minimization.\" the final sentence is not making a statement, suggest to reword as \"... $\\forall x \\in \\R^N$ s.t. $x^\\top 1 = 0$ we have $x^\\top C x \\geq 0$...\"\n-- Having both notations C_X and C_{\\bf X} is not advisable, they are too similar.\n-- Can you elaborate on what is meant by \"reasonable assumptions\"?\n\n-- Figure 1 is not that easy to understand without a bit of guidance from the authors in terms of what is being shown. Even when it is discussed in the text, it was not completely clear to me. Can you give more detailed description?\n\n-- Some more discussion around Theorem 4.1 would be nice, especially whether the graph $\\bar C$ is meaningful in its own right when $T$ induces a hard clustering.\n\n-- I recognise the fact that using spectral clustering is the most \"obvious\" approach to use in DR-C and C-DR, however as mentioned above this will typically lead to meaningless results when producing a large number of \"prototypes\". Although arguably less natural, I strongly suggest using a more commonly used \"prototyping\" method. If the graph can be interpreted as/converted into a distance matrix then Partitioning Around Medoids (PAM) is the most obvious, otherwise, despite the fact that it doesn't directly connect with C_Z(Z) in DR-C, I would suggest kmeans.", "broader_impact_concerns": "None are apparent", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "GzoXiMzkF7", "reviewer": "Reviewer_YTGR", "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": "## Strength\n\n- Novel Theoretical Framework\n\nThe paper presents a new perspective that integrates dimensionality reduction (DR) and clustering into a single optimization problem using Gromov-Wasserstein (GW) optimal transport. The proposed framework is theoretically well-founded, with formal theorems and proofs that establish its validity. The framework can recover traditional DR methods (e.g., PCA, UMAP, t-SNE) and clustering methods (e.g., k-means, spectral clustering), making it highly flexible.\n\nDifferent from the sequential approaches, DistR learns both simultaneously, so that everything will be optimized at the same time. The kind of idea is better in identifying representative prototypes in a way that reduces information loss, potentially leading to more meaningful cluster assignments.\n\n- Efficient Optimization Approach\n\nThe paper uses an efficient computational strategy: Block Coordinates Descent (BCD) Optimization that allows for stable training. It uses semi-relaxed GW Divergence (srGW) that improves computational feasibility, making the problem more tractable.\n\n- Interesting perspective for working on multi-view learning\n\nIn the Conclusion section, the authors explicitly mention that their approach naturally enables addressing multi-view settings with multiple unaligned inputs of different sizes​. The framework’s versatility allows for addressing multi-view settings, accommodating multiple unaligned inputs of different sizes. This seems to be a new direction for unsupervised learning on multi-views.\n\n## Weakness\n\n- Lack of scalability evidence for large datasets. \n\nWhile there are couple of experiments for the scalability problem for different number of prototypes, there lacks sufficient evidences for DistR to applied in the high number of samples. Based on the Table 1 from Appendix F.1, the maximum number of samples used for validation is 10,000 samples. This is away smaller than most of the dimension reduction methods. For instance, the original MNIST Datasets contains 60,000 training points, and the 10,000 testing points, so I am wondering how the computational time would increase as the training dataset grows up (same for F-MNIST). I am interested to see how the relationship between the number of samples and the runtime would be like. Even though the authors have discussed about the low-rank factorization in their main paper, there is no empirical experiments to show how efficient it would improve.\n\n- The interpretability and explainability of the prototypes\n\nIn the paper, the authors mentioned that the prototypes that DistR have offered is a discriminative representation of the input samples at various granularities, and from the numerical experiments perspective, they are trying to use the homogeneity as a kind of explanations. However, the metrics seems to be still unclear to show whether the prototypes are meaningful or just artifacts of optimization. A higher number of fixed prototypes generally increases the likelihood of achieving higher homogeneity, as each prototype can better capture distinct class structures, reducing class mixing within clusters. When you increase the number of prototypes, you allow for finer granularity in clustering, meaning each class can be better separated. If the number of prototypes significantly exceeds the number of true classes, homogeneity might not continue to improve, as some prototypes could start capturing sub-clusters within a single class, potentially reducing purity. Since the subclusters numbers are hard to evaluate at this point, I am a bit concerned about if those metrics are proving the prototypes are meaningful. Moreover, Figure 4 is showing the images that are medoids for each cluster, which is the argmax point for each cluster instead of the standard deviation or similarity within the cluster, so I am doubt that if the diversity within each group shrink (for example of the standard deviation within one subcluster, and the standard deviation between two cluster I assume should have statistically difference between each other?) This problem might lead to a doubt in showing whether the prototypes are just artifacts of the optimization or not?\n\n- The stopping criteria/hyperparameter selection for the DistR\n\nBased on the results shown within the Figure 5 and also for the Appendix F.4, the results for the DistR varies a lot when the number of prototypes changed, I am wondering when we are selecting the proper hyperparameters for the DistR, how we are able to select the hyperparameter for them. It would be more likely to be applicable if those information can be added within the paper.", "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": "This paper introduces a new framework, Distributional Reduction (DistR), which is able to integrate dimension reduction and clustering together based on the Gromov-Wasserstein (GW) optimal transport problem. This proves the researchers a new perspective, trying to jointly learn the low-dimensional embeddings and representative prototypes. The paper provides formal proofs establishing connections between GW transport, spectral embeddings, and prototype-based clustering. It used the block coordinate descents algorithm for the semi-relaxed GW divergence computation to guarantee convergence to the local optimum. From the evaluation perspective, the authors have compared eight different datasets to demonstrate its effectiveness compared to the DR->C and C->DR methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "To strengthen the paper and make it more compelling for both theoretical and practical applications, here are a few changes that I would like to recommend the authors to do.\n\n- Improve Scalability Validation with Large Datasets\n\nWhile the paper discusses low-rank factorization as a way to improve efficiency, there is no empirical validation demonstrating that DistR can scale to large datasets. The current experiments use a maximum of 10,000 samples, which is significantly smaller than datasets typically used in dimensionality reduction and clustering benchmarks. It would be better to conduct scalability experiments on larger datasets (e.g., 100,000+ samples) to demonstrate that DistR remains computationally feasible, and provide runtime vs. dataset size comparisons, similar to what is done in existing DR methods papers. Empirically validate the efficiency of low-rank factorization, showing how much computational and memory savings it provides in practice.\n\n- Improving the convincing of the learned prototypes\n\nThe paper claims that DistR’s prototypes offer meaningful and structured representations, but no clear evidence is provided to confirm their interpretability. The homogeneity metric is used to evaluate the prototypes, but it does not prove that the prototypes represent distinct, meaningful classes rather than being artifacts of optimization. It would be better to conduct a stability analysis and an analysis on intra-clusters. For the stability analysis, I am wondering if the author can provide an evaluation on whether learned prototypes remain consistent across multiple training runs. For the intra-cluster analysis, aside from Figure 4, I am wondering if the authors can provide an analysis on the distribution of the intra-cluster variance to see if the clusters are truly distinct. This can also helps for explaining why the DistR would go over the true class labels in some sense.\n\n- Provide Clearer Guidelines for Hyperparameter Selection\n\nI would recommend the authors to discuss more on the hyperparameter sensitivity analysis explanations and the corresponding guidelines for how to select a proper hyperparameter with the metrics curve given. In the real world application, the instability of the trend in some of the datasets might lead to the confusion of the users. It would be better if the authors can have some discussion on this hyperparameter selection topic, so that the users are more capable of accepting this idea for combining the dimension reduction and clustering in the same place.", "broader_impact_concerns": "I think the ethical claims that the authors might need to made for the broader impact statement would be mostly from the detection of the typical prototypes which is also the reason why I requested for a few changes for the large datasets.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)