SEMI-RELAXED GROMOV-WASSERSTEIN DIVERGENCE WITH APPLICATIONS ON GRAPHS
Year: 2022

Abstract:
Comparing structured objects such as graphs is a fundamental operation involved in many learning tasks. To this end, the Gromov-Wasserstein (GW) distance, based on Optimal Transport (OT), has proven to be successful in handling the specific nature of the associated objects. More specifically, through the nodes connectivity relations, GW operates on graphs, seen as probability measures over specific spaces. At the core of OT is the idea of conservation of mass, which imposes a coupling between all the nodes from the two considered graphs. We argue in this paper that this property can be detrimental for tasks such as graph dictionary or partition learning, and we relax it by proposing a new semi-relaxed Gromov-Wasserstein divergence. Aside from immediate computational benefits, we discuss its properties, and show that it can lead to an efficient graph dictionary learning algorithm. We empirically demonstrate its relevance for complex tasks on graphs such as partitioning, clustering and completion.

heading: 1 INTRODUCTION
text: One of the main challenges in machine learning (ML) is to design efficient algorithms that are able to learn from structured data (Battaglia et al., 2018). Learning from datasets containing such nonvectorial objects is a difficult task that involves many areas of data analysis such as signal processing (Shuman et al., 2013), Bayesian and kernel methods on graphs (Ng et al., 2018; Kriege et al., 2020) or more recently geometric deep learning (Bronstein et al., 2017; 2021) and graph neural networks (Wu et al., 2020). In terms of applications, building algorithms that go beyond Euclidean data has led to many progresses, e.g. in image analysis (Harchaoui & Bach, 2007), brain connectivity (Ktena et al., 2017), social networks analysis (Yanardag & Vishwanathan, 2015) or protein structure prediction (Jumper et al., 2021).
Learning from graph data is ubiquitous in a number of ML tasks. A first one is to learn graph representations that can encode the graph structure (a.k.a. graph representation learning). In this domain, advances on graph neural networks led to state-of-the-art end-to-end embeddings, although requiring a sufficiently large amount of labeled data (Ying et al., 2018; Morris et al., 2019; Gao & Ji, 2019; Wu et al., 2020). Another task is to find a meaningful notion of similarity/distance between graphs. A way to address this problem is to leverage geometric or signal properties through the use of graph kernels (Kriege et al., 2020) or other embeddings accounting for graph isomorphisms (Zambon et al., 2020). Finally, it is often of interest either to establish meaningful structural correspondences between the nodes of different graphs, also known as graph matching (Zhou & De la Torre, 2012; Maron & Lipman, 2018; Bernard et al., 2018; Yan et al., 2016) or to find a representative partition of the nodes of a graph, which we refer to as graph partitioning (Chen et al., 2014; Nazi et al., 2019; Kawamoto et al., 2018; Bianchi et al., 2020).
Optimal Transport for structured data. Based on the theory of Optimal Transport (OT, Peyré & Cuturi, 2019), a novel approach to graph modeling has recently emerged from a series of works. Informally, the goal of OT is to match two probability distributions under the constraint of mass conservation and in order to minimize a given matching cost. OT originally tackles the problem
of comparing probability distributions whose supports lie on the same metric space, by means of the so-called Wasserstein distance. Extensions to graph data analysis were introduced by either embedding the graphs in a space endowed with Wasserstein geometry (Nikolentzos et al., 2017; Togninalli et al., 2019; Petric Maretic et al., 2019) or relying on the Gromov-Wasserstein (GW) distance (Mémoli, 2011; Sturm, 2012). The latter approach is a variant of the classical OT in which one aims at comparing probability distributions whose supports lie on different metric spaces by finding a matching of these distributions being as close as possible to an isometry. By expressing graphs as probability measures over spaces specific to their topology, where the structure of a graph is represented through a pairwise distance/similarity matrix between its nodes, the GW distance computes both a soft assignment matrix between the nodes of the two graphs and a notion of distance between them (see the left part of Figure 1). These properties have proven to be useful for a wide range of tasks such as graph matching and partitioning (Xu et al., 2019a; Chowdhury & Needham, 2021; Vayer et al., 2019), estimation of nonparametric graph models (graphons, Diaconis & Janson, 2007; Xu et al., 2021) or for graph dictionary learning (Vincent-Cuaz et al., 2021; Xu, 2020). GW was also extended to directed graphs in Chowdhury & Mémoli (2019) and to labeled graphs via the Fused Gromov-Wasserstein (FGW) distance in Vayer et al. (2020).
Despite those recent successes, applications of GW to graph modeling still have several limitations. First, finding the GW distance remains challenging, as it boils down to solving a difficult nonconvex quadratic program (Peyré et al., 2016; Solomon et al., 2016; Vayer et al., 2019), which, in practice, limits the size of the graphs that can be processed. A second limit naturally emerges when a probability mass function is introduced over the set of the graph nodes. The mass associated with a node refers to its relative importance and, without prior knowledge, each node is either assumed to share the same probability mass or to have one proportional to its degree. In this paper, we somehow argue that these choices can be suboptimal in several cases, and should be relaxed, with the additional benefit of lowering the computational complexity. As an illustration, consider Figure 1: on the left image, the GW matching is given between two graphs, with respectively two and three clusters, associated with uniform weights on the nodes. By relaxing the weight constraints over the second (middle image) or first graph (right image) we obtain different matchings, that can better preserve the structure of the source graph by reweighing the target nodes and thus selecting a meaningful subgraph.
Contributions. We introduce a new optimal transport based divergence between graphs derived from the GW distance. We call it the semi-relaxed Gromov-Wasserstein (srGW) divergence. After discussing its properties and motivating its use in ML applications, we propose an efficient solver for the corresponding optimization problem. Our solver better fits to modern parallel programming than exact solvers for GW do. We empirically demonstrate the relevance of our divergence for graph partitioning, Dictionary Learning (DL), clustering of graphs and graph completion tasks. With srGW, we recover SOTA performances at a significantly lower computational cost compared to methods based on pure GW.

heading: 2 MODELING GRAPHS WITH THE GROMOV-WASSERSTEIN DIVERGENCE
text: In this section we introduce more formally the GW distance and discuss two of its applications on graphs, namely graph partitioning and unsupervised graph representation learning. In the following the probability simplex with N-bins is denoted as ΣN := {h ∈ RN+ | ∑ i hi = 1}.
GW as graphs similarity. In the OT context, a graph G with n nodes can be modeled as a couple (C,h) where C ∈ Rn×n is a matrix encoding the connectivity between nodes and h ∈ Σn is a histogram, referred here as distribution, modeling the relative importance of the nodes within the graph. The matrix C can be arbitrarily chosen as the graph adjacency matrix, or any other matrix describing the relationships between nodes in the topology of the graph (e.g. adjacency, shortest-path, Laplacian). The distribution is often considered as uniform (h = 1n1n) but can also convey prior knowledge e.g. using the normalized degree distribution (Xu et al., 2019a). Consider now two graphs G = (C,h) and G = (C,h), respectively with n and m nodes, potentially different (n 6= m). The GW distance between G and G is defined as the result of the following optimization problem:
GW22(C,h,C,h) = min T1m=h T>1n=h ∑ ijkl ∣∣Cij − Ckl∣∣2 TikTjl with T ∈ Rn×m+ (1) The optimal coupling T ? acts as a probabilistic matching of nodes which tends to associate pairs of nodes that have similar pairwise relations in C and C respectively, while preserving masses h and h through its marginals constraints. GW defines a distance on the space of metric measure spaces (mm-spaces) invariant to measure preserving isometries (Mémoli, 2011; Sturm, 2012). In the case of graphs, which corresponds to discrete mm-spaces, such invariant is a permutation of the nodes that preserves the structures and the weights of the graphs (Vayer et al., 2020). First applications of GW to ML on graphics problems (Peyré et al., 2016; Solomon et al., 2016) motivated further connections with graph partitioning. In the GW sense, this paradigm is illustrated by finding an ideally partitioned graph G = (D,h) of m clustrers whose structure is a diagonal matrix D ∈ Rm×m, representing the cluster’s connections, and its distribution h estimates the proportion of the nodes in each cluster (Xu et al., 2019a). The OT plan between a graph G = (C,h) represented by its adjacency matrix C to this ideal graph can be used to recover G’s clusters, if their number and proportions match the number and the weights of G’s nodes. Xu et al. (2019a) suggested to empirically refine both graph distributions, by choosing h based on power-law transformations of the degree distribution of G and to deduce h from h by linear interpolation. Chowdhury & Needham (2021) proposed to use heat kernels from the Laplacian of G instead of its adjacency binary representation, and proved that the resulting GW partitioning is closely related to the well-known spectral clustering (Fiedler, 1973).
The GW distance has been extended to graphs with node attributes (typically Rd vectors) thanks to the Fused Gromov-Wasserstein distance (FGW) (Vayer et al., 2019). In this context, an attributed graph G with n nodes is a tuple G = (C,F ,h) where F ∈ Rn×d is its matrix of node features. FGW between two attributed graphs aims at finding an OT by minimizing a weighted sum of a GW cost on structures and a Wasserstein cost on features (balanced by a parameter α ∈ [0, 1]). Most of the applications of GW can be naturally extended with FGW to attributed graphs.
GW for Unsupervised Graph Representation Learning. More recently, GW has been used as a data fitting term for unsupervised graph representation learning by means of Dictionary Learning (DL) (Mairal et al., 2009; Schmitz et al., 2018). While DL methods mainly focus on vectorial data (Ng et al., 2002; Candès et al., 2011; Bobadilla et al., 2013), DL applied to graphs datasets consists in factorizing them as composition of graph primitives (or atoms) encoded as {(Ck,hk)}k∈[[K]]. The first approach proposed by (Xu, 2020) consists in a non-linear DL based on entropic GW barycenters. On the other hand, Vincent-Cuaz et al. (2021) proposed a linear DL method by modeling graphs as a linear combination of graph atoms thus reducing the computational cost. In both cases, the embedding problem that consists in the projection of any graph on the learned graph subspace requires solving a computationally intensive optimization problem.

heading: 3 SEMI-RELAXED GROMOV-WASSERSTEIN DIVERGENCE

heading: 3.1 DEFINITION AND PROPERTIES
text: Given two observed graphs G = (C,h) and G = (C,h) of n and m nodes, we propose to find a correspondence between them while relaxing the weights h on the second graph. To this end we introduce the semi-relaxed Gromov-Wasserstein divergence as :
srGW22(C,h,C) = min h∈Σm GW22(C,h,C,h) (2)
This means that we search for a reweighing of the nodes of G leading to a graph with structure C with minimal GW distance from G. While the optimization problem (2) above might seem complex to solve, it is actually equivalent to a GW problem where the mass constraints on the second marginal of T are relaxed, reducing the problem to:
srGW22(C,h,C) = min T1m=h ∑ ijkl |Cij − Ckl|2TikTjl with T ∈ Rn×m+ . (3)
From an optimal coupling T ? of problem (3), the optimal weights h ?
expressed in problem (2) can be recovered by computing T ?’s second marginal, i.e h ? = T ?>1n. This reformulation with relaxed marginal has been investigated in the context of the Wasserstein distance (Rabin et al., 2014; Flamary et al., 2016) and for relaxations of the GW problem in (Schmitzer & Schnörr, 2013) but was never investigated for the GW distance itself. To the best of our knowledge, the most similar related work is the Unbalanced GW Séjourné et al. (2020); Liu et al. (2020); Chapel et al. (2019) where one could recover srGW with different weighting over the marginal relaxations (∞ on the first marginal and 0 on the second) but this specific case was not discussed nor studied in these works.
A first interesting property of srGW is that since h is optimized in the the simplex Σm, its optimal value h ? can be sparse. As a consequence, parts of the graph G can be omitted in the comparison, similarly to a partial matching. This behavior is illustrated in the Figure 1, where two graphs with uniform distributions and structures C and C forming respectively 2 and 3 clusters are matched. The GW matching (left) between both graphs forces nodes of different clusters from C to be transported on one of the three clusters of C, leading to a high GW cost where clusters are not preserved. Whereas srGW can find a reasonable approximation of the structure of the left graph either though transporting on only two clusters (middle) or finding a structure with 3 clusters in a subgraph of the target graph with two clusters (right). A second natural observation resulting from the dependence of srGW2 of only one input distributions is its asymmetry, i.e. srGW2(C,h,C) 6= srGW2(C,h,C). Interestingly, srGW2 shares similar properties than GW as described in the next proposition:
Proposition 1 Let C ∈ Rn×n and C ∈ Rm×m be distance matrices and h ∈ Σn with supp(h) = [[n]]. Then srGW2(C,h,C) = 0 iff there exists h ∈ Σm with card(supp(h)) = n and a bijection σ : supp(h)→ [[n]] such that:
∀i ∈ supp(h), h(i) = h(σ(i)) (4) And: ∀k, l ∈ supp(h)2, Ckl = Cσ(k)σ(l). (5)
In other words srGW2(C,h,C) vanishes iff there exists a reweighing h ? ∈ Σm of the nodes of the second graph which cancels the GW distance. When it is the case, the induced graphs (C,h) and (C,h ? ) are isomorphic (Mémoli, 2011; Sturm, 2012). We refer the reader interested in the proofs of the equivalence and the Proposition 1 to the annex (section 7.2).

heading: 3.2 OPTIMIZATION AND ALGORITHMS
text: In this section we discuss the computational aspects of the srGW divergence and propose an algorithm to solve the related optimization problem (3). We also discuss variations resulting from entropic or/and sparse regularization of the initial quadratic problem.
Algorithm 1 CG solver for srGW 1: repeat 2: G(t) ← Compute gradient w.r.t T of (2). 3: X(t) ← minX1m=h
X≥0 〈X,G(t)〉
4: T (t+1) ← (1−γ?)T (t)+γ?X(t) with γ? ∈ [0, 1] from exact-line search. 5: until convergence.
Solving for the semi-relaxed GW. The optimization problem related to the calculation of srGW22(C,h,C) is a non-convex quadratic program similar to the one of GW with the important difference that the linear constraints are independent.
Consequently, we propose to solve (3) with a Conditional Gradient (CG) algorithm (Jaggi, 2013) that can benefit from those independent constraints and is known to converge to local stationary point on non-convex problems (Lacoste-Julien, 2016). This algorithm, provided in Alg.
1, consists in solving at each iteration (t) a linearization 〈X,G〉 of the problem (3) whereG is the gradient of the objective in (3). The solution of the linearized problem provides a descent direction X? − T , and a linesearch whose optimal step can be found in closed form to update the current solution T (Vayer et al., 2019). While both srGW and GW CG require at each iteration the computation of a gradient with complexity O(n2m+m2n), the main source of efficiency of our algorithm comes from the computation of the descent directions. In the GW case, one needs to solve an exact linear OT problem, while in our case, one just needs to solve several independent linear problems under a simplex constraint. This simply amounts to finding the minimum on the rows ofG as discussed in (Flamary et al., 2016, Equation (8)), within O(mn) operations, with potential parallelization with GPUs. Performance gains are illustrated in the experimental section.
Entropic regularization. Recent OT applications have shown the interest of adding an entropic regularization to the exact problem (3) ,e.g. (Cuturi, 2013; Peyré et al., 2016). Entropic regularization makes the optimization problem smooth and more robust while densifying the optimal transport plan. Similar to (Peyré et al., 2016; Xu et al., 2019b; Xie et al., 2020), we can use a mirror-descent scheme w.r.t. the Kullback-Leibler divergence (KL) to solve entropic regularized srGW. The problem boils down to find, at iteration t the coupling T (t+1) ← arg minT1m=h;T≥0〈T ,G
(t)〉 + KL(T |T (t)) where > 0, G(t) denotes the gradient of the GW loss at iteration t and KL(T |T (t)) is the KL divergence. These updates can be solved using the following closed-form Bregman Projections :
T (t+1) ← arg min T1m=h;T≥0
KL(T |K(t+1)) ⇔ T (t+1) ← diag ( h
K(t)1m
) K(t). (6)
whereK(t) = exp ( G(t) − log(T (t)) ) (exp and log are applied componentwise). Unlike existing solvers for GW based on entropic regularization (Peyré et al., 2016; Solomon et al., 2016) that rely on a Sinkhorn’s matrix scaling algorithm onK(t) at each iteration, our problem requires only one (left) scaling ofK(t) per iteration. We denote by srGWe(C,h,C) the result of this procedure.
Sparsity promoting regularization. As illustrated in Figure 1, srGW naturally leads to sparse solutions in h̄. To compress even more the localization over a few nodes of C, we can promote the sparsity of h through a penalization Ω(T ) = ∑ j |hj |1/2 which defines a concave function in the positive orthant R+. We adapt the Majorisation-Minimisation (MM) of Courty et al. (2014) that was introduced to solve classical OT with a similar regularizer. The resulting algorithm, which relies on a local linearization of Ω(T (t)), consists in iteratively solving the srGW or srGWe problems, regularized at iteration t + 1 by a linear OT cost of components R(t)i,j = λg 2 (h (t) j ) −1/2. Further detailed explanations on these algorithms can be found in section 7.3 of the annex.

heading: 4 LEARNING THE TARGET STRUCTURE
text: A dataset of K graphs D = {(Ck,hk)}k∈[[K]] is now considered, with heterogeneous structures and a variable number of nodes, denoted by (nk)k∈[[K]]. In the following, we introduce a novel graph dictionary learning (DL) whose peculiarity is to learn a unique dictionary element. Then we discuss how this dictionary can be used to perform graph completion, i.e. reconstruct the full structure of a graph from an observed subgraph.

heading: 4.1 A NEW GRAPH DICTIONARY LEARNING
text: Semi-relaxed Gromov-Wasserstein embedding. We first discuss how an observed graph can be represented in a dictionary with a unique element (or atom) C ∈ Rm×m, assumed to be known or designed through expert knowledge. First, one computes srGW22(Ck,hk,C) using the algorithmic solutions and regularization strategies detailed in Section 3.1. From the resulting optimal coupling T ?k , the optimal weights for the target graphC are recovered with h ? k = T ?> k 1nk . The graph (C,h ? k) can be seen as a projection of (Ck,hk) in the GW sense and the distribution on the nodes h ?
k is an embedding of the graph (Ck,hk). Representing a graph as a vector of weights h ?
k on a graph C is a new and elegant way to define a graph subspace that is orthogonal to other DL methods that either rely on GW barycenters (Xu, 2020) or linear representations (Vincent-Cuaz et al., 2021). One particularly
Algorithm 2 Stochastic update of the dictionary atom C
1: Sample a minibatch of graphs B := {(C(k),h(k))}k. 2: Get transports {T ?k }k∈B from srGW(Ck,hk,C) with Alg.1. 3: Compute the gradient ∇̃C of srGW with fixed {T ?k }k∈B and perform a projected gradient step
on symmetric non-negative matrices S: C ← ProjS(C − η∇̃C) (7)
interesting aspect of this modeling is that when h ?
k is sparse, only the subpart (or subgraph) of C corresponding to the nodes with non-zero weights in h ?
k is used.
srGW Dictionary Learning and online algorithm. Given a dataset of graphs D = {(Ck,hk)}k∈[K], we propose to learn the graph dictionary C ∈ Rm×m from the observed data, by optimizing:
min C∈Rm×m
1
K K∑ k=1 srGW22(Ck,hk,C). (8)
This problem is denoted as srGW Dictionary Learning. It can be seen as a srGW barycenter problem (Peyré et al., 2016) where we look for a graph structure C for which there exists node weights (h ?
k)k∈[[K]] leading to a minimal GW error. Interestingly this DL model requires only to solve the srGW problem to compute the embedding h ?
k since it can be recovered from the solution T ?
k of the problem.
We solve this non-convex optimization problem above with an online algorithm similar to the one first proposed in Mairal et al. (2009) for vectorial data and adapted by Vincent-Cuaz et al. (2021) for graph data. The core of the stochastic algorithm is depicted in Alg. 2. The main idea is to use batches of graphs to independently solve the embedding problems (via Alg.1), then compute estimates of the gradient ∇̃C with respect to C on each batch B. Finally a projected gradient on the set S of symmetric non-negative matrices is performed to update C. In practice we use Adam optimizer (Kingma & Ba, 2015) in all our experiments. The complexity of this stochastic algorithm is mostly bounded by computing the gradients, which can be done in O(n2km+m
2nk) (see Section 3.2). Hence, in a factorization context i.e. maxk (nk) >> m, the overall learning procedure has a quadratic complexity w.r.t. the maximum graphs size. Since the embedding h ?
k is a by-product of computing the different srGW, we do not need an iterative solver to estimate it. Consequently, it leads to a speed up on CPU of 2 to 3 orders of magnitude compared to our main competitors (see Section 5.3) whose DL methods, instead, require such iterative scheme.

heading: 4.2 DL-BASED MODEL FOR GRAPHS COMPLETION
text: The structureC estimated on the datasetD can be used to infer/complete a new graph from the dataset that is only partially observed. In this setting, we aim at recovering the full structure C ∈ Rn×n while only a subset of relations between nobs < n nodes is observed, denoted as Cobs ∈ Rnobs×nobs . This amounts to solving:
min Cimp
srGW22
( C̃,h,C ) , where C̃ = [ Cobs
... . . . Cimp
] , (9)
where only the n2 − n2obs coefficients collected into Cimp are optimized (and thus imputed). We solve the optimization problem above by a classical projected gradient descent. At each iteration we find an optimal coupling T ? of srGW that is used to calculate the gradient of srGW w.r.t Cimp. The latter is obtained as the gradient of the srGW cost function evaluated at the fixed optimal coupling T ? by using the Envelope Theorem (Bonnans & Shapiro, 2000). The projection step is here to enforce known properties of C, such as positivity and symmetry. In practice the estimated Cimp will have continuous values, so one has to apply a thresholding (with value 0.5) on Cimp to recover a binary adjacency matrix. The method can be easily extended to labeled graphs by also optimizing the node features of non-observed nodes.

heading: 5 NUMERICAL EXPERIMENTS
text: This section illustrates the behavior of srGW on graph partitioning (a.k.a. nodes clustering), clustering of graphs and graph completion. All the Python implementations in the experiments will be released on Github. For all experiments we provide e.g. validation grids, initializations and complementary metrics in the annex (see sections: partitioning 7.5.1,clustering 7.5.2, completion 7.5.3).1

heading: 5.1 GRAPH PARTITIONING
text: As discussed in Section 2, it is possible to achieve graph partitioning via OT by estimating a GW matching between the graph to partition G = (C,h) (n nodes) and a smaller graph G = (Iq,h), with q n nodes. The atom Iq is set as the identity matrix to enforce the emergence of densely connected groups (i.e. communities). The distribution h estimates the proportion of nodes in each cluster. We recall that h must be given to compute the GW distance, whereas it is estimated with srGW. All partitioning performances are measured by Adjusted Mutual Information (AMI, Vinh et al., 2010).
Simulated data. In Figure 2, we illustrate the behavior of GW and srGW partitioning on a toy graph, simulated according to SBM (Holland et al., 1983) with 3 clusters of different proportions. We see that miss-classification occurs either when the proportions h do not fit the true ones or when q 6= 3. On the other hand, clustering from srGW can simultaneously recover any cluster proportions (since it estimates them) and can select the actual number of clusters, using the sparse properties of h.
Real-world datasets. In order to benchmark the srGW partitioning on real (directed and undirected) graphs, we consider 4 datasets (details provided in Table 4 of the annex): a Wikipedia hyperlink network (Yang & Leskovec, 2015); a directed graph of email interactions between departments of a European research institute (Yin et al., 2017); an Amazon product network (Yang & Leskovec, 2015); a network of interactions between indian villages (Banerjee et al., 2013). For the directed graphs, we adopt the symmetrization
procedure described in Chowdhury & Needham (2021). Our main competitors are the two GW based partitioning methods proposed by Xu et al. (2019b) and Chowdhury & Needham (2021). The former (GWL) relies on adjacency matrices, the latter (SpecGWL) adopts heat kernels on the graph normalized laplacians (SpecGWL). The GW solver of Flamary et al. (2021) was used for these methods. For fairness, we also consider these two representations for srGW partitioning (namely srGW and srSpecGW). Finally, three competing methods specialized in graph partitioning are also considered: FastGreedy (Clauset et al., 2004), Louvain (with validation of its resolution parameter,
1code available at https://github.com/cedricvincentcuaz/srGW.
Blondel et al., 2008) and Infomap (Rosvall & Bergstrom, 2008). The graph partitioning results are reported in Table 1. Our method srGW always outperforms the GW based approaches and on this application the entropic regularization seems to improve the performance. We want to stress that our general purpose divergence srGW outperforms methods that have been specifically designed for nodes clustering tasks on 3 out of 6 datasets.

heading: 5.2 CLUSTERING OF GRAPHS DATASETS
text: Datasets and methods. We now show how the embeddings (h?k)k∈[[K]] provided by the srGW Dictionary Learning can be particularly useful for the task of graphs clustering. We considered here three types of datasets (details provided in Table 9 of the annex): i) social networks from IMDB-B and IMDB-M (Yanardag & Vishwanathan, 2015); ii) graphs with discrete features representing chemical compounds from MUTAG (Debnath et al., 1991) and cuneiform signs from PTC-MR (Krichene et al., 2015); iii) graphs with continuous features, namely BZR, COX2 (Sutherland et al., 2003) and PROTEINS, ENZYMES (Borgwardt & Kriegel, 2005). Our main competitors are the following OT-based SOTA models: i) GDL (Vincent-Cuaz et al., 2021) and its regularized version, namely GDLλ; ii) GWF (Xu, 2020), with both fixed (GWF-f) and random (GWF-r, default setting for the method) atom size; iii) GW kmeans (GW-k), a k-means equipped with GW distances and barycenters (Peyré et al., 2016).
Experimental settings. For all experiments we follow the benchmark proposed in Vincent-Cuaz et al. (2021). The clustering performances are measured by means of Rand Index (RI, Rand, 1971). The standard Euclidean distance is used to implement k-means over srGW and GWFs embeddings, but we use for GDL the dedicated Mahalanobis distance as described in Vincent-Cuaz et al. (2021). GW-k does not use any embedding since it directly computes (a GW) k-means over the input graphs. For each parameter configuration (number of atoms, number of nodes and regularization parameter, detailed in section 7.5.2) we run each experiment five times, independently. The mean RI over the five runs is computed and the dictionary configuration leading to the highest RI for each method is reported.
Results and discussion. Clustering performances and running times are reported in Tables 2 and 3, respectively. All variants of srGW DL are at least comparable with the SOTA GW based methods. Remarkably, the sparsity promoting variants always outperform other methods. Notably Table 3 shows embedding computation times of the order of the millisecond for srGW, two order of magnitude faster than the competitors.

heading: 5.3 GRAPHS COMPLETION
text: Finally, we present graph completion results on the real world datasets IMDB-B and MUTAG, using the approach proposed in 4.2. Since this completion problem has never been investigated by existing GW graph methods, we adapted the learning procedure used for srGW to GDL (Vincent-Cuaz et al., 2021).
Experimental setting. Since the datasets do not explicitly contain graphs with missing nodes, we proceed as follow: first we split the dataset into a training dataset (Dtrain) used to learn the dictionary and a test dataset (Dtest) reserved for the completion tasks. For each graph of C ∈ Dtest, we created incomplete graphs Cobs by independently removing 10% and 20% of their nodes, uniformly
at random. The partially observed graphs are then reconstructed using the procedure described in Section 4.2 and the average performance of each method is computed for 5 different dataset splits.
Results. The graph completion results are reported in Figure 3. Our srGW dictionary learning outperforms GDL consistently, when enough data is available to learn the atoms. When the proportion of train/test data varies, we can see that the linear GDL model that maintains the marginal constraints tends to be less sensitive to the scarcity of data. This can come form the fact that srGW is more flexible thanks to the optimization of h but can slightly overfit when few data is available. Sparsity promoting regularization can clearly compensate this overfitting and systematically leads to the best completion performances (high accuracy, low Means Square Error).

heading: 6 CONCLUSION
text: We introduce a new transport based divergence between structured data by relaxing the mass constraint on the second distribution of the Gromov-Wasserstein problem. After designing efficient solvers to estimate this divergence, called the semi-relaxed Gromov-Wasserstein (srGW), we suggest to learn a unique structure to describe a dataset of graphs in the srGW sense. This novel modeling can be seen as a Dictionary Learning approach where graphs are embedded as a subgraph of a single atom. Numerical experiments highlight the interest of our methods for graph partitioning, and unsupervised representation learning whose evaluation is conducted through clustering and completion of graphs.
We believe that this new divergence will unlock the potential of GW for graphs with unbalanced proportions of nodes. The associated fast numerical solvers allow to handle large size graph datasets, which was not possible with current GW solvers. One interesting future research direction includes an analysis of srGW to perform parameters estimation of stochastic block models. Also, as relaxing the second marginal constraint in the original optimization problem gives more degrees of freedom to the underlying problem, one can expect dedicated regularization schemes, over e.g. the level of sparsity of h, to address a variety of application needs. Finally, our method can be seen as a special relaxation of the subgraph isomorphism problem. It remains to be understood theoretically in which sense this relaxation holds.

heading: ACKNOWLEDGMENTS
text: This work is partially funded through the projects OATMIL ANR-17-CE23-0012, OTTOPIA ANR20-CHIA-0030 and 3IA Côte d’Azur Investments ANR-19-P3IA-0002 of the French National Research Agency (ANR). This research was produced within the framework of Energy4Climate Interdisciplinary Center (E4C) of IP Paris and Ecole des Ponts ParisTech. This research was supported by 3rd Programme d’Investissements d’Avenir ANR-18-EUR-0006-02. This action benefited from the support of the Chair ”Challenging Technology for Responsible Energy” led by l’X Ecole polytechnique and the Fondation de l’Ecole polytechnique, sponsored by TOTAL. This work is supported by the ACADEMICS grant of the IDEXLYON, project of the Universit de Lyon, PIA operated by ANR-16-IDEX-0005. The authors are grateful to the OPAL infrastructure from Université Côte d’Azur for providing resources and support.

heading: 7 APPENDIX

heading: 7.1 NOTATIONS & DEFINITIONS
text: In this section we recall the notations used in the rest of the appendix.
Notations. For a vector h ∈ Rm we define its support as supp(h) = {i ∈ [[m]]|hi 6= 0}. Note that if h ∈ Σm we have supp(h) = {i ∈ [[m]]|hi > 0}. The cardinal of a discrete set A is denoted |A|. For a 4-D tensor L = (Lijkl)ijkl we denote ⊗ the tensor-matrix multiplication such that for a given matrixM , L⊗M is the matrix with entries ( ∑ kl LijklMkl)ij .
Constraints on couplings. We introduce U(h,h) the set of all admissible couplings between h and h, i.e the set
U(h,h) := {T ∈ Rn×m+ |T1m = h,T>1n = h}.
We also introduce for any h ∈ Σn and m ∈ N∗, U(h,m) the set of all admissible couplings between h and any histogram of Σm, i.e the set
U(h,m) := {T ∈ Rn×m+ |T1m = h},
such that ∀T ∈ U(h,m), the second marginal h(= T>1n) belongs to Σm.
Gromov-Wasserstein distance. For any q ≥ 1, the Gromov-Wasserstein distance of order q between (C,h) and (C,h) is defined by:
GWqq(C,h,C,h) = min T∈U(h,h) Lq(C,C,T )
where Lq(C,C,T ) = ∑ i,j∈[[n]]2 ∑ k,l∈[[m]]2 |Cij − Ckl|qTikTjl.
This definition can be equivalently expressed with a tensor-matrix multiplication as
GW qq (C,h,C,h) = min T∈U(h,h) 〈Lq(C,C)⊗ T ,T 〉 (10)
where Lq(C,C) is the 4-D tensor such that Lq(C,C) = ( (Cij − Ckl)q ) ijkl .
Semi-relaxed Gromov-Wasserstein divergence. The semi-relaxed Gromov-Wasserstein divergence of order q satisfies the following equation:
srGWqq(C,h,C) = min T∈U(h,m) Lq(C,C,T )
or equivalently srGWqq(C,h,C) = min
T∈U(h,m) 〈Lq(C,C)⊗ T ,T 〉. (11)
Fused Gromov-Wasserstein distance. For attributed graphs i.e graphs with nodes features, we use the Fused Gromov-Wasserstein distance (Vayer et al., 2019). Consider two attributed graphs G = (C,F ,h) and G = ( C,F ,h ) where F = (Fi)i∈[[n]] ∈ Rn×d and F = (F j)j∈[[m]] ∈ Rm×d are their respective matrix of features whose rows are denoted by {Fi}i∈[[n]] and {F i}i∈[[m]]. Given a trade-off parameter between structures and features denoted α ∈ [0; 1], the FGW qq,α distance is defined as the result of the following optimization problem
min T∈U(h,h) (1− α) ∑ ij ‖Fi − Fj‖qqTij + α ∑ ijkl |Cij − Ckl|qTikTjl (12)
which is equivalent to
min T∈U(h,h)
〈(1− α)Mq(F ,F ) + αLq(C,C)⊗ T ,T 〉 (13)
whereMq(F ,F ) is a matrix with entriesMq(F ,F )ij = ‖Fi − Fj‖qq .
Semi-relaxed Gromov-Wasserstein distance. Then in a similar way than for the GromovWasserstein distance, the semi-relaxed Fused Gromov-Wasserstein of order q is defined for any trade-off parameter α ∈ [0; 1] as the result of the following optimization problem
srFGW qq,α(C,F ,h,C,F ) = min T∈U(h,m) 〈(1− α)Mq(F ,F ) + αLq(C,C)⊗ T ,T 〉 (14)
One can see these problems as regularized versions of the quadratic problem of GW/srGW where a linear term in T takes into account a similarity measure between features of attributed graphs.

heading: 7.2 SEMI-RELAXED (FUSED) GROMOV-WASSERSTEIN PROPERTIES

heading: 7.2.1 EQUIVALENCE OF THE OPTIMIZATION PROBLEMS
text: Let us begin with the proof of the equivalence between both optimization problems used to define our semi-relaxed Gromov-Wasserstein divergence. Consider two observed graph (C,h) and (C,h) of n and m nodes. Our first optimization problem reads as
h ? ∈ arg min
h∈Σm GWqq(C,h,C,h). (15)
The second optimization problem coming from the relaxation of the second marginal constraints (T>1n = h) of admissible couplings T reads as
min T∈U(h,m) ∑ ijkl |Cij − Ckl|2TikTjl. (16)
Note that in the following proofs we always assume that graphs correspond to finite discrete measures and that their geometry is well-defined in the sense that entries of C and C are always finite. This implies that graphs define compact spaces ensuring existence of optimal solutions for both problems.
Lemma 1 Problems 15 and 16 are equivalent.
Proof. Consider a solution of problem 15 denoted (h?1,T1) note that the definition implies that T1 ∈ U(h,h ? 1). Another observation is that given h ?
1, the transport plan T1 also belongs to arg minT∈U(h,h?1) Lq(C,C,T ) hence is an optimal solution of GWq(C,h,C,h ? 1).
Now consider a solution of problem 16 denoted T2 with second marginal h ?
2. By definition the couple (h ?
2,T2) is suboptimal for problem (15) i.e
Lq(C,C,T1) ≤ Lq(C,C,T2). (17)
And the symmetric also holds as T1 is a suboptimal admissible coupling for problem 16 i.e ,
Lq(C,C,T2) ≤ Lq(C,C,T1). (18)
These inequalities imply that Lq(C,C,T1) = Lq(C,C,T2). Therefore we necessarily have T1 ∈ arg minT∈U(h,m) Lq(C,C,T ), and (h ?
2,T2) ∈ arg minh∈Σm,T∈U(h1,h) Lq(C,C,T ). Hence the equality, GWq(C,h,C,h ? 1) = GWq(C,h,C,h ?
2), holds true. Therefore by double inclusion we have
arg min T∈U(h,m) Lq(C,C,T ) = arg min h∈Σm,T∈U(h,h) Lq(C,C,T ). (19)
Which is enough to prove that both problems are equivalent.
The equivalence between analog problems for the semi-relaxed Fused Gromov-Wasserstein divergence can be proved in the exact same way.

heading: 7.2.2 ZERO CONDITIONS OF (FUSED) SRGW
text: We prove the following result:
Proposition 1 Let C ∈ Rn×n and C ∈ Rm×m be distance matrices and h ∈ Σn with supp(h) = [[n]]. Then srGWq(C,h,C) = 0 iff there exists h ∈ Σm with card(supp(h)) = n and a bijection σ : supp(h)→ [[n]] such that:
∀i ∈ supp(h), h(i) = h(σ(i)) (20)
and: ∀k, l ∈ supp(h)2, Ckl = Cσ(k)σ(l). (21)
Algorithm 3 CG solver for srGW, optionally with a linear regularization termD ∈ Rn×n. Note that D = 0 for unregularized version of srGW of equation 25.
1: repeat 2: G(t) ← Compute gradient w.r.t T of (24) satisfying equation 26 applied in T (t). 3: Get direction X(t) problem: Solve independent subproblems on rows of G(t) thanks to the
equivalence stated here X(t) ← arg min
X∈U(h,m) 〈X,G(t) +D〉
⇔X(t) = (hix(t)i )i∈[[n]] ← arg min xi∈Σm
〈xi,G(t)i +Di〉. (22)
4: Get optimal step size γ? for the descent directionX(t) − T (t): for any γ ∈ [0, 1] let us denote Z(t)(γ) = T (t) + γ(X(t) − T (t)), then γ? is defined as
arg min γ∈[0,1]
〈L(C,C)⊗Z(t)(γ),Zt(γ)〉+ 〈D,Z(t)(γ)〉 (23)
factored as a second-order polynomial function of the form aγ2 + bγ + c, by using linearity of the tensor-multiplication ⊗ and of the scalar product. Then the closed form is obtained by simple analysis of coefficients a and b as in (Vayer et al., 2019).
5: T (t+1) ← Z(t)(γ?) = (1− γ?)T (t) + γ?X(t). 6: until convergence decided based on relative variation of the loss between step t and t+ 1.
Proof. The reasoning involved in this proof mostly relates on the definition of srGW as minh∈Σm GW q q(C,h,C,h).
(⇒) Assume that srGWq(C,h,C) = 0. Then we have GWq(C,h,C,h) = 0 for some h. By virtue to Gromov-Wasserstein properties (Sturm, 2012, Lemma 1.10) there exists a bijection σ between the support of the distributions which is distance preserving. In other words, there exists σ : supp(h)→ supp(h) = [[n]] such that h(j) = h(σ(j)) for all j ∈ supp(h) and Ckl = Cσ(k)σ(l) for all k, l ∈ supp(h)2.
(⇐) Consider T ∈ U(h,h) and the induced bijection σ satisfying equations 20 and 21. It is trivial to verify that L(C,C,T ) = 0 implying that GWq(C,h,C,h) = 0. Moreover as T ∈ U(h,m) since U(h,h) ⊂ U(h,m) and the same cost is involved in both transport problems, we have:
0 ≤ srGWq(C,h,C) ≤ GWq(C,h,C,h ? ) = 0 =⇒ srGWq(C,h,C) = 0.

heading: 7.3 ALGORITHMIC DETAILS
text: We provide in this section the algorithmic details completing our explanations given in the main paper (see subsection 3.2). Note that for all numerical experiments we considered q = 2, so the following algorithms are specific to this scenario.

heading: 7.3.1 CONDITIONAL GRADIENT SOLVER FOR SRGW AND SRFGW
text: The optimization problem related to computing srGW22(C,h,C) can be reformulated as
min T1m=h T≥0
vec(T )> ( C 2 ⊗K 1n1>n − 2C ⊗K C ) vec(T ). (24)
where ⊗K denotes the Kronecker product of two matrices, vec the column-stacking operator and the power operation on C is applied element-wise. For simplicity let us use the other equivalent formulation using the 4-D tensor notation i.e
min T∈U(h,h)
〈L(C,C)⊗ T ,T 〉 (25)
where L(C,C) = ( (Cij − Ckl)2 ) ijkl
. This is a non-convex optimization problem with the same objective and gradientG than GW which can be expressed in the general case as
G = L(C,C)⊗ T +L(C>,C>)⊗ T . (26)
Note that if C and C are symmetric matrices the gradient can be factored to
G = 2L(C,C)⊗ T . (27)
We propose to use a Conditional Gradient (CG) algorithm Jaggi (2013) to solve this problem, provided in Alg. 3.
For the sake of conciseness, we also introduce here a linear regularization term of the form 〈D,T 〉 which will be used while enforcing sparsity promoting regularization on srGW or for adapting this algorithm to srFGW whereD results from distances between features (up to proper scaling with the trade-off parameter α).
Then in this general case, our CG algorithm consists in solving at each iteration t the following linearization of the problem (24)
min X1m=h X≥0
〈X,G(t) +D〉 (28)
whereG(t) is the gradient at T (t) of (24). The optimization problem above can be very efficiently solved as discussed in (Flamary et al., 2016, Equation (8)). Indeed the problem above can clearly be reformulated as n independent linear problems under simplex constraints (each row ofX can be solved independently) of the form
min x>1m=hr,x≥0
x>gr (29)
where gr is the row r ofG. Optimizing a linear function over the simplex of dimensionality n can be done in O(m) because it consists in finding the smallest component i? = arg min gr in the linear cost, the solution being a scaled dirac vector hrδi? where all the mass is positioned on component i?. The solution of the linearized problem provides a descent directionX(t) − T (t), and a line-search is performed to get the optimal step size as described in equation 23. Note that the gradients at time step T (t) and T (t+1) have a simple linear relation, which we used in our implementation for conciseness but it is rather equivalent to the more classical implementation literally expressing terms involved in the line-search part.
Extension to srFGW. For two attributed graphs G = (C,F ,h) and (C,F ,h), of n and m nodes respectively, the semi-relaxed Fused Gromov-Wasserstein divergence of order 2 from G onto G is defined for any α ∈ [0, 1] as the result of the following optimization problem
srFGW2α(C,F ,h,C,F ) = min T∈U(h,m) 〈(1− α)M(F ,F ) + αL(C,C)⊗ T ,T 〉 (30)
where M(F ,F ) denotes the matrix of euclidean distances between features of F and F , i.e. Mij = ‖Fi − F j‖22. As a side note, one efficient way to compute this matrix in practice is by using the following factorizationM(F ,F ) = (F F )1d1>m+1n1>d (F F )−2FF > , where denotes the Hadamard product. The problem in equation 14 is a nonconvex quadratic problem which gradient w.r.t. T reads as : αG+ (1− α)M(F ,F ) (31) whereG corresponds to the gradient of the GW cost satisfying equation 26. Therefore we propose to tackle this problem by a straight-forward adaptation of the CG algorithm 3 detailed above by adding the multiplier α to the GW cost and settingD = (1− α)M(F ,F ).

heading: 7.3.2 ALGORITHMIC DETAILS AND GUARANTIES ON ENTROPIC SRGW.
text: Recent OT applications have shown the interest of adding an entropic regularization to the exact problem (24) Cuturi (2013); Peyré et al. (2016). We detail here how to design a mirror-descent scheme w.r.t. the Kullback-Leibler divergence (KL) to solve (24), in the same vein than Peyré et al. (2016); Xu et al. (2019b); Xie et al. (2020).
Algorithm 4 MD solver for entropic srGW, optionally with a linear regularization termD ∈ Rn×n. Note thatD = 0 for unregularized version of srGW of equation 25.
1: repeat 2: G(t) ← Compute gradient w.r.t T of (24) satisfying equation 26 applied in T (t). 3: Compute the matrixK(t)( ) following equations (36) and (38). 4: Get T (t+1) with the scaling ofK(t)( ) following equation (37). 5: until convergence decided based on relative variation of the loss between step t and t+ 1.
Algorithmic details. Indeed in order to solve the srGW problem stated in 16, we can use a mirrordescent scheme w.r.t. the KL geometry. At iteration t, the update of the current transport plan T (t) ∈ U(h,m) results from the following optimization problem
T (t+1) ← arg min T∈U(h,m) 〈G(T (t)),T 〉+ KL(T |T (t)) (32)
whereG satisfies 26 and KL(T |T (t)) = ∑ ij Ti,j log( Ti,j
T (t) i,j
)− Ti,j + T (t)i,j . Let us denote the entropy of any T ∈ Rn×m+ by H(T ) = − ∑ ij Tij(log Tij − 1), then the following relation can be proven
〈M ,T 〉 − H(T ) = KL(T | exp(−M ))⇔ KL(T |M) = 〈− logM ,T 〉 − H(T ) (33)
and leads to this equivalent formulation of 32:
T (t+1) ← arg min T∈U(h,m) 〈G(T (t))− logT (t),T 〉 − H(T ). (34)
DenotingM (t)( ) = G(T (t))− logT (t), overall we end up with the following iterations
T (t+1) ← arg min T∈U(h,m) 〈M (t)( ),T 〉 − H(T ). (35)
Which is equivalent thanks to the relation stated above to
T (t+1) ← arg min T∈U(h,m) KL(K(t)( )|T ) where K(t)( ) := exp{−M (t)( )/ }. (36)
Following the seminal work of (Benamou et al., 2015), the optimal T (t+1) is given by a simple scaling of the matrixK(t)( ) reading as
T (t+1) = diag
( h
K(t)( )1m
) K(t)( ). (37)
Note that an analog scheme is achievable while penalizing the srGW problem with a linear term of the form 〈D,T 〉. Which would simply result in a modification of the matrixM (t)( ) such that
M (t)( ) = G(T (t)) +D − logT (t). (38)
This more general setting is summarized in algorithm 4.
Similarly than for the CG algorithm 3, it is straight-forward to adapt the MD algorithm 4 to the semi-relaxed Fused Gromov-Wasserstein using FGW cost expressed in 31.

heading: 7.3.3 SPARSITY PROMOTING REGULARIZATION OF SRGW
text: In order to promote sparsity of the estimated target distribution while matching (C,h) to the structure C using srGW, we suggest to use Ω(T ) = ∑ j |hj |1/2 which defines a concave function in the positive orthant R+. This results in the following optimization problem:
min T∈U(h,m)
〈L(C,C)⊗ T ,T 〉+ λgΩ(T ) (39)
with λg ∈ R+ an hyperparameter. As mentioned in the main paper, equation 39 can be tackled with a Majorisation-Minimisation (MM) algorithm. MM consists in iteratively minimising an upper bound
Algorithm 5 MM solver for srGWg and srGWe+g
1: SetR(0) = 0. 2: repeat 3: Get optimal transport T (t) with second marginal h (t) from CG solver 3 (srGWg) or MD
solver 4 (srGWe+g) withD = R(t).
4: ComputeR(t+1) = λg2 (h (t) j ) −1/2 ij the new local linearization of Ω(T
(t)). 5: until convergence.
Algorithm 6 Stochastic update of the atom C
1: Sample a minibatch of graphs B := {(C(k),h(k))}k. 2: Get transports {T ?k }k∈B from srGW(Ck,hk,C) with Alg.1. 3: Compute the gradient ∇̃C of srGW with fixed {T ?k }k∈B and perform a projected gradient step
on symmetric non-negative matrices S: C ← ProjS(C − η∇̃C) (41)
of the objective function which is tight at the current iterate (Hunter & Lange, 2004). With this procedure, the objective function is guaranteed to decrease at every iteration. In our case, we only need to majorize the penalty term Ω(T ) to obtain a tractable function. Denoting h (t) = T (t)>1n, the estimate at iteration t, one can simply apply the tangent inequality∑ j √ hj ≤ ∑ j √ h (t) j + 1 2 ∑ j √ h (t) j (h− h(t))>1m. (40)
Using this inequality is equivalent to linearize the regularization term at h (t)
whose contribution can be absorbed into the inner product as 〈L(C,C)⊗ T +R(t),T 〉 where R(t) = λg2 (h (t) j ) −1/2 ij . The overall optimization procedure is summarized in 5. Note that the same procedure is used for promoting sparsity of srFGW using the adaptation of our CG solver 3 and MD solver 4 to this scenario.

heading: 7.4 LEARNING THE TARGET STRUCTURE
text: We detail in the following the algorithms for the srGW Dictionary Learning and its application to graphs completion.

heading: 7.4.1 SRGW DICTIONARY LEARNING.
text: We propose to learn the graph atom C ∈ Rm×m from the observed data D = {(Ck,hk)}k∈[K], by optimizing
min C∈Rm×m
1
K K∑ k=1 srGW22(Ck,hk,C), (42)
This is a nonconvex problem that we propose to tackle thanks to a stochastic gradient algorithm summarized in 6. At each iteration it consists in sampling a batch of graphs B = {(Ck,hk)}k∈[[B]] from the dataset and to embed each graph where hk = T>k 1n by solving independent srGW problems for the current state of the dictionary. Then we compute the estimate of the gradient over C reading as
∇̃C(.) = 2
B ∑ k∈[[B]] ( C hkh > k − T>k CkTk ) . (43)
Note that only the entries (i, j) of C such that i and j belongs to ∪k supp(hk) can have a non-null gradient. Finally, depending on the input structures we consider, we apply a projection on an adequate set S , for instance if input structures are adjacency matrices we consider S as the set of non-negative
Algorithm 7 Algorithms for graph completion using DL from srGW or GDL 1: Initialize randomly the entries Cimp by iid sampling from N (0.5, 0.01) (In a symmetric manner
if Cobs is symmetric). 2: repeat 3: Compute the optimal representations G̃(t) of the (C(t),h) onto the dictionary and optimal
transport T (t) between (C(t),h) and G̃(t): (srGW) : G̃(t) = (C,h (t) = T (t)>1n) where T (t) ← srGW(C(t),h,C)
(GDL) : G̃(t) = ( ∑ s∈[[s]] w(t)s Cs, h) where (w (t),T (t))← min w∈ΣS GW(C(t),h, ∑ s∈[[s]] wsCs,h)
(46) 4: Get C(t+1) from a projected gradient step w.r.t. the GW distance between (C(t),h) and G̃(t) with the optimal coupling T (t). 5: until convergence.
symmetric matrices (iterates can preserve symmetry but not necessarily non-negativity depending on the chosen learning rate).
Extension to attributed graphs. The stochastic algorithm described above can be adapted to a dataset of attributed graph D = {(Ck,Fk,hk)}k∈[[K]] with nodes features in Rd, by learning an attributed graph atom (C,F ). The optimization problem would now read as
min C∈Rm×m,F∈Rm×d
1
K ∑ k∈[[K]] srFGW22,α(Ck,Fk,hk,C,F ) (44)
for any α ∈ [0, 1]. The stochastic algorithm 6 is then adapted by first computing embeddings based on srFGW instead of srGW. Then simultaneously updating C following (44) up to the factor α and F thanks to the following equation
∇̃F (.) = 2
K ∑ k∈[[B]] ( diag(hk)F − T>k Fk ) . (45)

heading: 7.4.2 DL-BASED MODEL FOR GRAPHS COMPLETION
text: The estimated structure C, learned on the dataset D can be used to infer/complete a new graph form the dataset being only partially observed. Let us say that we want to recover the structure C ∈ Rn×n but we only observed the relations between nobs < n nodes denoted as Cobs ∈ Rnobs×nobs . We propose to solve the following optimization problems to recover the full matrix C while modeling graphs with our srGW Dictionary Learning
min Cimp
srGW22
( C̃,h,C ) , where C̃ = [ Cobs
... . . . Cimp
] , (47)
Graph completion with GDL For the linear dictionary GDL (Vincent-Cuaz et al., 2021) with several atoms denoted {Cs}s∈[[S]] sharing the same weights h we adapted our formulation to their model:
min Cimp min w∈Σs
GW 22 C,h, ∑ s∈[[S]] wsCs,h  (48) The matrix C has only the n2 − n2obs coefficients collected into Cimp are optimized (and thus imputed). This way Cimp expresses the connections between the imputed nodes, and between these nodes and the observed ones in Cobs. We tackle the non-convex problems above following the procedure described in Algorithm 7.
It consists in an alternating scheme where at each iteration t, we embed the current estimated graph C(t) on the dictionary depending on the chosen model (see (46)). Let us denote G̃(t) = (C̃(t), h̃(t))
this representation which respectively reads for srGW as G̃(t) = (C,h (t) ) and for GDL as G̃(t) = ( ∑ s∈[[S]] w (t) s Cs,h). Note that from this embedding step we also get the optimal transport T (t) from the GW distance between (C(t),h) and G̃(t). Then we can update the imputed part C(t)imp of C (t) while keeping fixed Cobs thanks to a projected gradient step with gradient w.r.t. C reading as
2 ( C(t) hh> − T (t)C̃(t)T (t)> ) (49)

heading: 7.5 DETAILS ON THE NUMERICAL EXPERIMENTS.

heading: 7.5.1 GRAPH PARTITIONING
text: Detailed partitioning benchmark. We detail here the benchmark between srGW and state-ofthe-art methods for graph partitioning on real (directed and undirected) graphs. We replicated the benchmark from (Chowdhury & Needham, 2021) using the 4 datasets whose preprocessing is detailed in their paper and resulting statistics are provided in Table 4. We considered the two GW based partitioning methods proposed by (Xu et al., 2019a) (GWL) and (Chowdhury & Needham, 2021) (SpecGWL). We benchmarked our methods denoted srGW and srSpecGW using respectively the adjacency and heat kernel on normalized Laplacian matrices as inputs. All these OT based methods depend on hyperparameters which can be tuned in an unsupervised way (i.e. without knowing the ground truth partition) based on modularity maximization (Chowdhury & Needham, 2021). Following their numerical experiments, we considered as input distribution h for the observed graph the parameterized power-law transformations of the form hi = pi∑
i pi where pi = (deg(i) +
a)b, with deg(i) the i-th node degree and real parameters a ∈ R and b ∈ [0, 1]. If the graph has no isolated nodes we chose a = 0 and a = 1 otherwise. b is validated within these 10 values {0, 0.0001, 0.005, ..., 0.1, 0.5, 1} which progressively transform the input distribution from the uniform to the normalized degree ones. An ablation study of this parameter is reported in Table 6. The heat parameter for SpecGWL and srSpecGW is tuned for each dataset within the range [1, 100] by recursively splitting this range into 5 values, find the parameter leading to maximum modularity, and repeat the same process on the induced new interval with this best parameter as center. This process is stopped based on the relative variation of maximum modularity between two successive iterates of precision 10−3. A similar scheme can be used to fine-tune the entropic parameter.
On our algorithm initialization. (Chowdhury & Needham, 2021) also discussed the sensitivity of GW solvers to the initialization and showed that GW matchings with heat kernels have considerably less spurious local minimum than GW applied on adjacency matrices. This way it is arguable to compare both methods by using the product hh > as default initialization. We observed for srGW that using an initialization based on the hard assignments of a kmeans on the rows of the input representations (up to the left scaling diag(h)), was a better trade-off for both kinds of representation. Hence we applied this scheme for our method in this benchmark. We illustrate in Table 4 how srGW refines these hard assignments by soft ones through OT.
Note that for these partitioning tasks we should/can not initialize the transport plan of our srGW solver using the product of h ∈ ΣN with a uniform target distribution h (0) = 1Q1Q, leading to
T (0) = 1Qh1 > Q. Indeed, for any symmetric input representation C of the graph, the partial derivative of our objective w.r.t the (p, q) ∈ [[N ]]× [[Q]] entries of T , satisfies
∂L ∂Tpq (C, IQ,T (0)) = 2 ∑ ij (Cip − δjq)2T (0)ij = 2 Q ∑ ij (C2ip + δjq − 2Cipδjq)hi
= 2 Q {Q ∑ i C2iphi + 1− 2 ∑ i Ciphi} (50)
This expression is independent of q ∈ [[Q]], so taking the minimum value over each row in the direction finding step of our CG algorithm (see algo 3) will lead toX = T (0). Then the line-search step involving, for any γ ∈ [0, 1], Z(0)(γ) = T (0) + γ(X − T (0)) will be independent of γ as Z(0)(γ) = T (0). This would imply that the algorithm will terminate with optimal solution T ? = T (0) being a non-informative coupling.
Parameterized input distributions: Ablation study. We report in Table 6 an ablation study of the parameter b ∈ [0, 1] introduced on the input graph distributions, first suggested by (Xu et al., 2019a). In most scenario and every OT based methods, a parameter b ∈]0, 1[ leads to best AMI performances (except for the dataset Village), while the common assumption of uniform distribution remains competitive. The use of raw normalized degree distributions consequently reduces partitioning performances of all methods. Hence these results first indicated that further research on the input distributions could be beneficial. They also suggest that the commonly used uniform input distributions provide a good compromise in terms of performances while being parameter-free.
Additional clustering metric. In order to complete our partitioning benchmark on real datasets, we report in 10 the Adjusted Rand Index (ARI). The comparison between ARI and AMI has been thoroughly investigated in (Romano et al., 2016) and led to the following conclusion: ARI should be used when the reference clustering has large equal sized clusters; AMI should be used when the reference clustering is unbalanced and there exist small clusters. Our method srGW always outperforms the GW based approaches for a given type of input structure representations (adjacency vs heat kernels) and on this application the entropic regularization seems to improve the performance. Note that this metric also leads to slight variations of rankings, which seem to
occur more often for methods explicitly based on modularity maximization. We want to stress that for this metric our general purpose divergence srGW outperforms methods that have been specifically designed for nodes clustering tasks on 4 out of 6 datasets (instead of 3 using the AMI).
srGW runtimes: CPU vs GPU for large graphs. Partitioning experiments with our methods were run on a GPU Tesla K80 as it brought a considerable speed up in terms of computation time compared to using CPUs as large graphs had to be processed. To illustrate this matter, we generated 10 graphs following Stochastic Block Models with 10 clusters, a varying number of nodes in {100, 200, ..., 2900, 3000} and the same connectivity matrix. We report in Table the averaged runtimes of one CG iteration depending on the size of the input graphs. For small graphs with a few hundreds of nodes, performances on CPUs or a single GPU are comparable. However, operating on GPU becomes very beneficial once graphs of a few thousand of nodes are processed.
Runtimes comparison on real datasets. We report in Table 8 the runtimes on CPUs for the partitioning benchmark on real datasets. For the sake of comparison, we ran the srGW partitioning experiments using CPUs instead of a GPU. We can observe that GW partitioning of raw adjacency matrices (GWL) is the most competitive in terms of computation time (while being the last in terms of clustering performances), on par with InfoMap. Our srGW partitioning methods lack behind the GW based ones in terms of speed. Note that we remain in the same order of complexity.
We want to stress that as illustrated in the previous paragraph, the use of CPUs for srGW multiplied by 10 to 20 times the computation times observed on a GPU. First the linear OT problems involved in each CG iteration of GW based methods are solved thanks to a C solver (POT’s implementation (Flamary et al., 2021)), whereas our srGW CG solver is fully implemented in Python. This can compensate the computational benefits of srGW solver regarding GW solver (see details in section 3.2 of the main paper). So the computation time became mostly affected by the number of computed gradients. We observed in these experiments that for the same precision level, srGW
needed considerably more iterations to converge than GW. Let us recall that our srGW partitioning method consists in matching an observed graph to the identity matrix C = I , whereas GW based partitioning methods use as target structure C = diag(h) as h is fixed beforehand. Our lack of heterogeneity in the chosen target structure could be the reason of these slow convergence patterns. hence advocates for further studies regarding srGW based graph partitioning using more informative structures.

heading: 7.5.2 DICTIONARY LEARNING: CLUSTERING EXPERIMENTS
text: We report in table 9 the statistics of the datasets used for our benchmark on clustering of many graphs.
Detailed settings for the clustering benchmark. We detail now the experimental setting used in the clustering benchmark derived from the benchmark conducted by Vincent-Cuaz et al. (2021). For datasets with attributes involving FGW, we validated 15 values of the trade-off parameter α
via a logspace search in (0, 0.5) and symmetrically (0.5, 1). For DL based approaches, a first step consists into learning the atoms. srGW dictionary sizes are tested in M ∈ {10, 20, 30, 40, 50}, the atom is initialized by randomly sampling its entries from N (0.5, 0.01) and made symmetric. The extension of srGW to attributed graphs, namely srFGW, is referred as srGW for conciseness in 2 of the main paper. One efficient way to initialize atoms features for minimizing our resulting reconstruction errors is to use a Kmeans algorithm seeking for M clusters on the nodes features observed in the dataset. For GDL and GWF, a variable number of S = βk atoms is validated, where k denotes the number of classes and β ∈ {2, 4, 6, 8}. The size of the atoms M is set to the median of observed graph sizes within the dataset for GDL and GWF-f. These methods initialize their atoms by sampling observed graphs within the dataset, with adequate sizes for GDL and GWF-f, while sizes of GWF-r are just determined by this random sampling procedure independently of the number of nodes distribution. For srGWg, srGWe+g andGDLλ, the coefficient of our respective sparsity promoting regularizers is validated within {0.001, 0.01, 0.1, 1.0}. Then for srGWe, srGWe+g, GWF-f and GWF-r, the entropic regularization coefficient is validated also within {0.001, 0.01, 0.1, 1.0}. Finally, we considered the same settings for the stochastic algorithm hyperparameters across all methods: learning rates are validated within {0.01, 0.001} while the batch size is validated within {16, 32}; We learn all models fixing a maximum number of epochs of 100 (over convergence requirements) and implemented an (unsupervised) early-stopping strategy which consists in computing the respective unmixings every 5 epochs and stop the learning process if the cumulated reconstruction error (Errt) does not improve anymore over 2 consecutive evaluations (i.e. Errt ≤ Errt+1 and Errt ≤ Errt+2). For the sake of consistency, we report in Table 7 the Averaged Mutual Information (AMI) performances on this benchmark, as we reported the RI in the main paper to be consistent with our main competitors.
Visualizations of srGW embeddings. We provide in Figure 5 some examples of graphs embeddings from the dataset IMDB-B learned on a srGW dictionary C of 10 nodes. We assigned different colors to each nodes of the graph atom (forth column) in order to visualize the correspondences recovered by the OT plan T ? resulting from the projection of the respective sample (C,h) ontoC in the srGW sense (third column). As the atom has continuous values we set the edges intensity of grey proportionally to the entries of C. By coloring nodes of the observed graphs based on the OT to its respective embedding (C)(second column) we clearly observe that key subgraphs information, such as clusters and hubs are captured within the embedding.
Impact of the graph atom size. Moreover if we increase the size of the dictionary, our embeddings are refined and can bring complementary structural information at a higher resolution e.g. finding substructures or variable connectivity between nodes in the same cluster. We illustrate these resolution patterns in Table 6 for embeddings learned on the IMDB-B dataset. We represent the embedded graph size distributions depending on the sizes of the graph atom learned thanks to srGW and its sparse variant srGWg. For any graph atom size, the mean of each embedded graph size distribution represented with a white dot is below the atom size, hence embeddings are sparse and subparts of the atom are indeed selected. Moreover, promoting sparsity of the embeddings (srGWg) lead to more concentrated embedded
graph size distributions with lower averaged sizes than its unregularized counterpart (srGW), as expected. Finally, these distributions seem to reach a stable configuration when the graph atom size is large enough. This argues in favor of the existence of a (rather small) threshold on the atom size where the heterogeneity of all the graphs contained in the dataset is well summarized in the dictionaries.

heading: 7.5.3 DICTIONARY LEARNING: COMPLETION EXPERIMENTS
text: Experiment details on graphs completion. For these completion experiments the exact same scheme, than in our clustering benchmark, is applied for learning srGW and GDL dictionaries on formed datasets Dtrain. On interesting specificity used on dataset IMDB-B was for the initialization. Indeed instead of initializing the structure atom iid from N (0.5, 0.01), we initialized the entries of Cimp denoting connections only between the new imputed nodes randomly. Then we initialized entries representating connections of these nodes to the ones of Cobs seeing the degree (scaled by maximum observed degrees) of each observed node as initial probability for a new node to be connected to this observed node. Then for each entries (i, j) where xi belongs to Cobs and xj to Cimp,Cij is initialized asN ( deg(xi)maxiobs deg(xi) , 0.01). This initialization led to better performances for both model so might be a first good practice to tackle completion of clustered graph using OT. For MUTAG we sticked with the initialization used for our and initialized imputed features randomly in the range of locally observed features.
Additional experiments on graphs completion. We complete the experiments on completion tasks of datasets IMDB-B and MUTAG reported in the subsection 5.3 of the main paper. Let us recall that for these experiments we fixed a percentage of imputed nodes (10 % and 20%) and looked at the evolution of completion performances over the proportion of train/test datasets. Here instead, we fix the proportion of the test dataset to 10% and make percentage of imputed nodes vary in {10, 15, 20, 25, 30}. A similar benchmark procedure than for experiments of the main paper is conducted.
These graph completion results are reported in Figure 7 for IMDB-B dataset and in Figure 8 for MUTAG dataset. Our srGW dictionary learning and its regularized variants outperform GDL and GDLλ consistently when the percentage of imputed nodes is not too high (< 20%), whereas this trend is reversed for high percentage of imputed nodes. Indeed, as srGW dictionaries capture subgraph patterns of variable resolutions from the input graphs, the scarcity of prior information in an observed graph leads to a too high number of valid possibilities to complete it. Whereas GDL dictionaries based on GW lead to more steady performances as they keep their focus on global stteructures. Interestingly, the sparsity promoting regularization can clearly compensate this kind of overfitting over subgraphs for higher levels of imputed nodes and systematically leads to better completion performances (high accuracy, low Means Square Error). Moreover, the entropic regularization of
srGW (srGWe and srGWe+g) can be favorably used to compensate this overfitting pattern for high percentages of imputed nodes (> 20%) and also pairs well with the sparse regularization (srGWe+g).

Review summaries:
The paper proposed a relaxed version of the Gromov-Wasserstein (GW) distance between two probability distributions on (graph) structured data. The main idea of new divergence is to relax the target distribution via an optimization setting. Based on the proposed divergence, the paper also introduced new formulations for existing problems of learning with graphs such as graph clustering (partitioning), graph dictionary learning, clustering of the graph, and graph completion. The experimental results for those applications on graph data have shown the efficiency of the proposed divergence.

In this paper the authors propose a modification of the Gromov-Wasserstein problem for matching given matrices $C$ and $\bar{C}$ with its respective histograms $h$ and $\bar{h}$, relaxing the restriction on $\bar{h}$. The problem is then given in an equivalent formulation, which is solved via Conditional Gradient.  The authors propose applications with graphs, such as graph partitioning, clustering, and completion via Graph Dictionary Learning. 

 The paper proposes to apply a version of Gromov-Wasserstein divergence to graphs. In particular, they relax the weight constraint of the second graph  and try to find a minimizer of the Gromov-Wasserstein distance. They named it Semi-relaxed Gromov Wasserstein distance (srGW) and apply it to solve some problems in graphs, namely, graph partition, graph clustering and graph completion. The experiments were carried out on synthetic data and real dataset such as: Wikipedia hyper link network, Amazon product network etc. The results show comparable and better performance over other methods such as  GW and spectral GW.

This paper proposes a semi-relaxed Gromov-Wasserstein (GW) dissimilarity for learning tasks on graphs. Here the "semi-relaxed" refers to removing one of the marginal constraints, with the practical effect that during GW graph matching, one is able to ignore some of the nodes in the source or target and thus obtain matchings that better respect the structure of the problem. What is quite interesting is that this simple relaxation immediately yields benefits in fundamental tasks such as graph partitioning, and the authors validate this observation with extensive experiments on graph datasets. The empirical results are strong, and show state-of-the-art performance compared to several recent baselines.