Linear-Time Gromov Wasserstein Distances using Low Rank Couplings and Costs
Authors: Meyer Scetbon
Gabriel Peyré
marco cuturi
Venue: NeurIPS 2021 Submitted (2021)
Keywords: Optimal Transport
Gromov-Wasserstein
Low-rank methods

Abstract: The ability to compare and align related datasets living in heterogeneous spaces plays an increasingly important role in machine learning. The Gromov-Wasserstein (GW) formalism can help tackle this problem. Its main goal is to seek an assignment (more generally a coupling matrix) that can register points across otherwise incomparable datasets. As a non-convex and quadratic generalization of optimal transport (OT), GW is NP-hard. Yet, heuristics are known to work reasonably well in practice, the state of the art approach being to solve a sequence of nested regularized OT problems. While popular, that heuristic remains too costly to scale, with cubic complexity in the number of samples $n$. We show in this paper how a recent variant of the Sinkhorn algorithm can substantially speed up the resolution of GW. That variant restricts the set of admissible couplings to those admitting a low rank factorization as the product of two sub-couplings. By updating alternatively each sub-coupling, our algorithm computes a stationary point of the problem in quadratic time with respect to the number of samples. When cost matrices have themselves low rank, our algorithm has time complexity $\mathcal{O}(n)$. We demonstrate the efficiency of our method on simulated and real data.

Linear-Time Gromov Wasserstein Distances using Low Rank Couplings and Costs
Anonymous Author(s)
Affiliation Address
email
Abstract
The ability to compare and align related datascts living in heterogeneous spaces plays an increasingly important role in machine learning. The Gromov- Wasserstein (GW) formalism can help tackle this problem. Its main goal is to seek assign- ment (more generally a coupling matrix) that can register points across otherwise incomparable datasets. As a non-convex and quadratic generalization of optimal transport (OT), GW is NP-hard. Yet, heuristics are known to work reasonably well in practice, the state of the approach being to solve sequence of nested regularized OT problems. While popular, that heuristic remains too costly to scale, with cubic complexity in the number of samples n.. We show in this paper how a recent variant of the Sinkhorn algorithm can substantially speed up the resolution of GW. That variant restricts the of admissible couplings those admitting low rank factorization as the product of two sub-coupe couplings. By updating alternatively each sub-coupling, our algorithm computes a stationary point of the problem in quadratic time with respect to the number of samples. When cost matrices have themselves low rank, our algorithm time complexity (n). We demonstrate the efficiency of our method on simulated and real data.
2
4
8
10
11
12
13
14
15
Introduction
18.
The ever increasing interest for Gromov-Wasse Wasserstein... Several problems in machine learning involve comparing families of points that live in heterogeneous spaces. This situation arises typically when realigning two distinct sets of feature representations obtained from the similar source. Recent applications single-cell genomics [15] and NLP [12, 1] provide two cases in point: Thousands of cells taken from the same tissue are split in two groups, each group is processed with different experimental protocol, resulting two distinct sets of heterogeneous feature vectors; Thousands of word embeddings for two languages are learned independently. In both cases, one expects to find
19
20
21
23
24
25 a meaningful way to register points across sets living in heteregencous spaces, since they contain similar overall information. That realignment is usually carried out using the -assoov-asssrsteiin (GW) machinery proposed by Mémoli [26] and Sturm [36], which seeks a relaxed assignment matrix that is "close" to isometry possible, using quadratic score to quantify that closeness. GW
27
29 has a lot practical appeal: It has been used in supervised learning [41], generative modeling [7], domain adaptation [9], structured prediction [37], quantum chemistry [27] and alignment layers [17].
31
despite its cubic cost. Because it is an NP-hard problem, these applications rely on approximating GW, typically by solving a sequence of OT problems using entropic regularization. This heuristic is efficient yet costly, since requires (n+3) operations to register two sets of n samples, a price that paid when re-instantiating each OT problem. Our goal is to reduce substantially that complexity by exploiting cow-factorization of both parameters (data) and variable (relaxed assignment) matrices in the GW problem, while maintaining state of the art performance in applications.
33
35
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.

---

37
Wasserstein: from cubic to linear complexity. A comparatively simpler problem is the registration of two populations embedded in the same space. This corresponds the classic optimal transport
38
39 (OT) problem, which has received considerable attention in ML [28]. OT has found applications
40 in computer vision [29], NL.P [24], single cell tracking [33] or multi-task regression in neuro- imaging [22]. While the OT problem is originally cast as a linear program, with a O( log(n)) cost, many of these works rely on solving instead a penalized OT problem using Sinkhorn's 's algorithm [34, 13]. In its most naive implementation, the Sinkhorn has quadratic complexity [2]. Recent works achieve (nn) complexity by targeting the matrix-vector updates Sinkhorn's 's algorithm using low-rank approximations of the data kernel matrix [4, 3, 31]. This idea be further improved by imposing the low-rank constraint on the optimization variables of the original OT problem [19]. to modify Sinkorn's 's steps by enforcing low rank factorization of the coupling variable [32].
41
42
43
44
46
Gromov- Wassosstsss::::: from NP-hard linear approximations. The GW problem replaces the linear objective function in OT by non-comvex quadratic objective. Much like OT is re- laxation of the optimal assignment problem, GW can be seen as relaxation of the quadratic assignment problem (QAP), Both GW and QAP are NP-hard to solve [8]. In practice, iteratively minimizing linearization of that quadratic objective using Sinkhorn works surprisingly well [20, 35]. This method corresponds to mirror-descent scheme [27]. 54 and in the special case of Euclidean distance matrices, the loss is concave and it can be also interpreted as bi-linear
48
49
50
51
52
M
56
relaxation [23]. In the most general case, this results in an
Onn*) algorithm (the objective is a quadratic function of n x n. relaxed assignment matrix), that is reduced to ((n)) when using separable losses [27]. price that remains too
58
60
high for several ML applications. It is possible to replace the GW distance by cheaper yet only distantly related prox- ies, such as lower bounds based on OT [26] (see also [30]) or sliced projections [38]. Whether GW can be efficiently sped up remains an open question. We propose in this work novel approach that leverages, as done recently for
61
62
63
64
65
86
OT, low-rank methods. A very recent line of works attacks this problem by quantizing first the two input spaces to solve GW problem of reduced size, thus effectively pro- ducing an ad-hoc low-rank coupling [11]. A nice feature of this approach is that it maintains the triangular inequal- ity and provides a valid upper-bound on the GW distance, Related approaches which also approximate GW distance
87
GW loss: 0.35357
GW loss: 0,35389
68
Figure 1: Top row: we compute the GW coupling between two curves in 2D and 3D, with n m = 10000 points. These points are endowed with the squared L2 distance. Bottom row: coupling obtained with the SoTA entropic approach [20, 271, compared with our linear method with rank r = 10. See Appendix D. for more details.
69
70
71
73 using clustering methods (possibly a recursive way) [6] and [40]. We take in this paper direct approach: 75 instead of separating clustering and GW resolution in 2 independent steps, propose do address them simulta- neously: our method seeks the least-costly (in GW sense) coupling with low rank constraint, illustrated in Fig.
77
79
Contributions We introduce the low-rank-GNV problem, by imposing a low rank constraint on feasible couplings. This method works hand-in-hand with entropic regularization and leads to Sinkhorn- like algorithm. Because of its exclusive reliance on matrix-vector products, the method streams well on GPUs. This method can also leverage low-rank factorizations of the input data matrices to further reduce the complexity of each iteration to reach linear time. Numerical evaluations simulated and real datasets show that this low-rank approximation maintains the favorable property of antropi-rregllaiie GW (namely its ability to compute "good" local minima) for a linear computational price, thus paving the way for larger scale uses of GW in ML.
81
83
85
2 Background on the Gromoy- Wassersttnss Framework
87
88 Comparing measured metric spaces. Let (X,dx) and (), dy: be two metric spaces, and l and two discrete probability measures on X' and V, respectively. We write μ := Σ a,d,, and
90 V: :: Em where n, m 2 1, a, b are two histograms in the probability simplicies An, Am of
respective size and m, and (I1.... In), (V1...., .... 9m) are two families in X' and y. For g 2 1.

---

92 let us also denote == (F)) E R=x" and B (d). Zy JIss, Sm E R R
93 two pairwise cost matrices between the points in the respective supports of p and v. The Gromov-
94 Wasserstein (GW) discrepancy between two discrete metric measure spaces (n, dx) and (v, dy) is
95 the solution of the following non-convex quadratic problem, instantiated here for simplicity as
96 function of (a, A) and (b, B), which contain all the information that is needed:
GW((a, A), (b, B)) min EA.B(P), where Ia, = {P R+ Pell..a
Pim P, b). (1)
87
and the energy AA is a quadratic function parameterized by a loss L: RR R
(P. = E (as)) 3.
98 A typical choice of the loss is the TP distance L(a, = b/2 with p 2 1. In that case, [26]
99 proves that G/ defines a distance on the space of metric measure spaces quotiented by measure-
100
preserving isometries. When p 2, as we consider from now on, the GW objective can be evaluated efficiently using the marginal constraints imposed on P. as follows [27]:
(P) = (A a, a} + (B b) 2(APB.P)
(3)
102
Indeed, (3) can be computed efficiently in O(n2m + nm²) operations, using only matrix/matrix
multiplications, instead of the O(n²m² complexity of the naive evaluation of (2).
103
104
Entropic Gromoy-Wasserve
The original GW problem (1) can be regularized using an entropic
105
term [20, 35, 27], leading to the following problem:
GW (la, (b, B)) = min (P(())HH((P PElla.,
-
(4)
where H(P) := : -i, (log((...)) 1) is the entropy of P. By applying a Mirror Descent (MD) scheme with respect the KL divergence and by choosing the step-size to be y I/E, Peyré et al. [27] provide simple algorithm which consists in solving a se-
106
107
108
109
quence of regularized OT problem as presented in Algorithm 1. Indeed, each KI. pro- jection in Algorithm can be computed efficiently thanks to the Sinkhorn algorithm [13].
110
111
112
Computational complexity. Given cost matrix C, the KL projection of K, onto the polytope II b). where KL(P, Q) (P, log(P/Q) 11), is carried out in the inner loop of Algo. using the Sinkhorn algorithm, through matrix-vector products This quadratic complexity (in red) is dominated by the cost of updating matrix C' each iteration in Algorithm 1, which requires O(n + nm?) algebraic operations (cubic, in violet). As noted above, evaluating the objective EA.a(P) (P) has the same order. In the following we show that by considering low rank exact decomposition (or approximation) of the distance matrices, the cubic cost of reupdating C and subsequently evaluating εA,B can be brought down to quadratic.
113
Algorithm 1 Entropic-GW
114
Inputs: A, B, b,E P = ab?
115
116
for l = 0,... do
117
C t- -4APB nm (n+m)
Kg f- exp(-C/ca nm P+ argmin (L Ks) ((m))
119
120
PELI(a.b)
121
end Result: (A,)(P) nm (n+m)
123
3 Exploiting Low-Rank Factorization for Cost Matrices
125
126
Exact factorization of cost matrices. In this section we consider the case where the cost matrices A and B admit low-rank factorization. More precisely, we make the following assumption.
128
Assumption 1. Assume that A and B admit low-rank factorization, that is there exists A1, Au E xdd and B B2 E Rm such that A = A AT and B = B, BI, where d & n,d' d' & m.
129
130
case in point is when both and are squared Euclidean distance matrices, with a sample size that is larger than ambient dimension. This case is highly relevant, covering many applications of OT to MI. The d & n assumption is also likely to hold for most applications, since cases where d > n
132
are known to pose challenges to the estimation of OT [16, 39J. Writing X [21...., zn] dx if
3

---

134
i.j" then one has, writing 2 (Xx2) R" that A = =1,, + 1== - -x Therefore by denoting A1 = |2, In. - v2x*] E 1nx (d-2) and A [1, v2x*] E (2x(d+2)
A = |
135
136
we obtain the factorization above.
137
Under Assumption 1, the complexity of Algo. is downgraded to quadratic in sample size: the two operations that make Algo. cubic lie in the updates of the cost and the computation of the objective. Observe that for any given P E IRnxm, one can compute at each iteration
138
139
C 4= B
140
in nm(d= + d²) + dn + m) algebraic operations. Moreover thanks to the reformulation of (P) given in (3), one can compute it in quadratic time as well. Indeed writing G1 :=
141
PB, and G2 == A P both in Rdxd", one has (APB, P) = (G, @ G2) Com- puting GI, G2 given P requires only 2(rmd + mdd"), and computing their dot product adds dd' algebraic operations. The overall complexity to compute EA,B(P) is O(nmd + mdd").
142
143
144
145
146
General distance matrices. When the original cost matrices A, are not low-rank but describe distances, we propose to use a recent body of work that output their low-rank approximation in linear time [5, 21J. These algorithms produce, for any distance matrix D Rxxm and > O, matrices D1 E 12axd D2 E Ruxdd in (( +
Algorithm 2 Quadratic Entropic-GW
Inputs: A2, B1, B2, a, b,s P = ab² nm
148
149
for l = 0.... do
150
G2 x Ab PB, nmd mdd' nmd ndd K₂ t exp(- - -C/E) nm
C t- -4A,
152
poly( ( )) algebraic operations such that with probability at least 0.99 one has
153
P argmin (L(P,Kg) (nm)
154
PeTT(a,))
DD - D2 17-- < I Calli- + *1
end
4 f (A02 a. a) + ? b) Onm)
where Ca denotes the best rank-d approximation
G2 A PBi nmd mdd?
156
to D. We fall back on this approach to obtain a
Gi + A PB2 amd mdd
low-rank factorization of a distance matrix in lin-
02 5- 221 @ G2)La O(dd')
158
ear time whenever needed, aware that this incurs
EA,B(P) + +
an additional approximation. See Appendix B
Return: (AP)(
180
for more details.
Imposing a Low Nonnegative Low-Rank for the Coupling
161
In this section, we shift our attention to a different opportunity for speed-ups, without assuming that Assumption holds: we regularize the GW problem problem by decomposing the coupling as a product of two low-rank couplings, in the footsteps of [18, 32), using the following definition:
163
165
Definition Given M E R=xm, the nonnegative (NN) rank of M is the smallest number of nonnegative rank-one matrices into which the matrix can be decomposed additively:
4
(k+ min g|M
R4 Vi, = R4 2
187 Following [18, 32], we propose to constrain GW, enforcing a rank r on the coupling:
(* A), (b, B)) := min EA,B(P), where Ia,b()) = (P E Na rk4 S (5) Pell,,(r)
188 Note that the minimum is always attained as (r,,(r) is compact and the objective is continuous.
169 In [32], the authors show that one can parameterize any coupling in (r,( as a product of two
170
low-rank couplings linked by a common marginal. For any g E A,, the interior of A,, writing
go p E Raxm P Q diag(1/9) Q E Ha.g and R E ll
171 one has that UgEA: Hg II. (r). Therefore GW-LR introduced in (5) can be reformulated as the following optimization problem
172
min GW-LR( ((a, A), (b. B)) (Q.R.))eC(a.brr)
εA,A(Q diag(1/g)
(6)

---

173
where C(a.b,r) == == Ci(a,b,r) (r with
Ci(a,b,r) = { (Q. R. g) E R4x x RTxT x (RT) s.t. Q1, = =a, R1, = -).
C2(r) = ( (12, R, E R R x R s.t. 1² = R² 1m = 99
174
Stabilization of the Method. [32] propose to stabilize the objective defined in (6) by adding to the constraints lower bound a on the weight vector g such that g 2 a coordinate-wise Indeed, as a solution of (6) must satisfies > 0 coordinate-wise, then for a sufficiently small, the solution of the same problem where one adds the constraint g 2 a will remain the same. Therefore let us introduce our new set of constraints C(a, b,r,a) a) == Ca, b, r, a) (2 where (a b, b, a)) :=
175
178
177
178
179
C, (a, b) {{Q, R. g) g 2 a}. Another way to stabilize the method is by considering a double regularization scheme as proposed in [32] where in addition of constraining the nonnegative rank of the coupling, we regularize the objective by adding an entropic term (Q. R, g). which is to be understood as that of the values of the three respective entropies evaluated for each term.
180
181
182
(W-LR!ya ((a, A), (b, B)) =
min (Q.R.g)EC(abbrra)
E diag(1//)))) (1)))(())()))))))))))))))))))))) - e (HQ, R,g))
(7)
Mirror Descent Scheme. As in [27], propose to use MD scheme with respect to the KL
divergence to approximate GW-LRR()) in (7). More precisely, for any e 2 O, the MD scheme leads for all & 2 0 to the following updates which require solving a convex barycenter problem per step:
184
185
(Qx+1, 9k+1) == argmin KL(Ç, K&)
(8)
(EC(a,b.r,))
where (Qu, 90) E b(a,b,r) is initial point such that Qu > 0 and Ro > 0. Pk = Qn diag /1 1/9x) K₁ == (K k( (()), K = exp(4y BRg diag(1/gk)
187
(7e 1) log((u)), )) K( (2) == BX DQx diag(1 1/gk) (ne - 1) log((R RR)), K(
188
exp(- / (7E 1) log(sx)) with wale =- 1Q RR)s, for all i (I...., r) and
189
190
is a positive step size. Solving (8) can be done efficiently thanks to the Dykstra's Algorithm as
191
showed in [32]. See Appendix C for more details.
192
Initialization. To initialize our algorithm, we adapt the First Lower Bound of [26] to our case of interest. More preciscly, we show the following Proposition. See appendix A for the proof.
193
Proposition 1. Let us denote ā A R., ÿ = B02b R. and C ((ē; ;))))as E R
194
195
Then for all & 2 0 and r 2 1 we have,
GW-LR (r) ((a, (b, B)) 2
(č, Q diag(1/9) -H( min (Q.R..)C((.bbrra))
(9)
196
Note that the RHS of the inequality (9) is exactly the problem studied in [32] for which an algorithm was proposed. Therefore to initialize our algorithm, we propose to use their approach. Note that here
the cost Č' is the squared Euclidean distance between two families (#1,.. (.. f.) } and (ūI...., Dm } in I-D which admits a low-rank factorization. Therefore we can apply the linear-time version of the algorithm presented in [32] to compute the solution. Algorithm 3 summarizes our approach, where
198
199
200
201
(-)) denotes the operator extracting the diagonal of a square matrix.
Computational Cost. Computing the initialization goes through the computations of ä and ÿ which requires On2² + =²) algebraic operations. Moreover, applying the algorithm proposed in [32] when the underlying cost is the squared Euclidean distances between two families in 1-D needs only O(n + m)r) algebraic operations. Solving the barycenter problem as defined in (8) can be done efficiently thanks to Dykstra's 's Algorithm. Indeed in [32, Algorithm 2] the authors show that
203
205
given (K(2) k k((3)) ), each iteration of their algorithm requires only O(nn + m)r) algebraic operations since it involves only matrix/vector multiplications. However computing the kernel
207
matrices (K(") [2(2)) x(²)) at each iteration of Algorithm requires quadratic complexity with
209
210
respect to the number of samples. Overall the proposed algorithm, while faster than the cubic
implementation proposed [27], still needs ((([2++²])] operations per iteration. the following we will see that by combining both nonnegative low-rank constraints on the coupling and low-rank approximations of the distance matrices, we can obtain a linear time algorithm with respect to the number of samples which computes an approximation of the GW distance.
211
213
5

---

Algorithm 3 Low-Rank GW, GW-LR (3)
((a, A), (b, B))
Eid
Inputs: A, B, a, b, E, I f A-4%,
B
O(m² +
+ Step (+)
21 f 502, 22 - 702 (m++n)
G + [2,, -v2]. C2 + [m 22 227"
O(n+m)
(Q, R, +
argmin (Q. (E...)))((...rrra))
(C, Q (iag(L/g)
= R)))
O(nn m))
for k L. ... do
C - - AQ diag(1/g) C2 + RBB O((n² + m))r)
+ Step (*x)
K()) + exp(4xC0 CyRdiag
(7s 1) log(())
((m +
K(2) + exp(4>02 Qiadiag(L/g) - (nE 1) log( (t))
((
w L-D(Q""COO C2C,), K( expp-A/w/gp" (re (1e--1)og(00 l))g)))) O(rr?)
Q, argmin KL(S, (k(l), K(2), ((())) n)r
(EC(a,brr,)
end
C c - (ī, a) (g.b) n+m
f- + diag(1/g), C2 + RTB ((n2+m²)z) + Step (**)
G- C2R, G +- CLG. G Q" diag(1/9) O((m+
02 - GT) I
- 0 +
Return: 3
Convergence of the mirror descent. Even if the objective (7) is not convex in (Q, R. g), we obtain the non-asymptotic stationary convergence of the MD algorithm in this setting. For that purpose 217 we consider the same convergence criterion as the one proposed in [32] to obtain non-asymptotic
215
218
stationary convergence of the MD scheme defined as
A (5, 7) = 72 (KL GE,a (L(G.
where (e,a(E. 7) argmincec(o.b.rr)) (VEA, + KL(().)) For any 1/r 2 a > 0, we
show in the following proposition the non-asymptotic stationary convergence of the MD scheme applied to the problem (7). See Appendix A for the proof.
221
Proposition 2. Let E 2 0, a > and 2 1. By denoting Le,a := 27([|]]a2|BBlaaaa + +.)
222
223
and by considering a constant stepsize in the MD scheme (8) 2Le,n'
we obtain that
min A (1Q RR, g2),)) Do I
N
(7) ((a, A), (b, B)) is the distance of the initial value where D= EA,B E (Qo diag GW-LR
224
225
the optimal one.
c.a ((a, A), (b, B)) GW-LR!")((e A), B)). Thus Proposition 2 show that our algorithm reach stationary point of (7). In particular, if 0, the proposed algorithm converges towards stationary point of (5).
Recall that for a sufficiently small, we have GW-LR (r)
226
228
Double Low-rank Approach for Linear Time GW
229
Almost all operations in Algorithm 3 are linear time, except for the three updates highlighted in red, involving C1 and C2, and the computations of ã = A a and ÿ = B as they still require a quadratic number of algebraic operations. When adding Assumption i from $3 to the rank constrained approach from $4, we notice that the strengths of both approaches can work hand in hand, both in casier initial evaluations of ã, ÿ, but, most importantly, each new recomputation of a factorized linearization of the quadratic objective:
231
233
234
235
Linear time outer norms. Because A admits a low-rank factorization. one can obtain a low-rank factorization for _ Indeed, remark that for 2, E Rd y) Eije1 Tı,,W.y. Therefore

---

by studying the rows of matrices A1 == [af) :_.;;aa) and A2 := ; ".)),, if one writes
(x)) == Vect(ax")) E RRd where Vect(-) the vectorization operation, we obtain that
A = A T where A = [<((aaa")
(a!")),, Ae [ (a
In Algorithm 3, the line "Step (+)" can thus be replaced by I + A
T a and ā + B3 B2 b Note that computing A₁ given A1 requires only (nd?) operations, so that this alternate code only takes (nd))) + (m()))² operations.
236
237
238
239
Linear time linearization of the GW objective. The linearization step, the critical step in Algo.
that consists in updating C' at each iteration, consumes substantial portion of the computational budget of GW. Introducing the low-rank Sinkhorn approach makes this step quadratic in Algo.3: 3; the complexity of that step also quadratic using the low-rank assumption on costs A and B, in Algo.2. 2. There is therefore an opportunity to marry both to speed-up that important step. We argue that this is indeed what happens, in the sense that combining the two yields indeed linear time complexities in sample sizes, by replacing in Algorithm 3, the lines "Step (**)" by
241
242
243
244
245
C1 + - A1 Q diag(1/g) and C2 + BBB
Note that this speed-up would not be achieved using other approaches that output low rank approximation of the transport plan [4, 3, 31J. The crucial obstacle to using these methods here is that the cost matrix C' in GW is "synthetic". in the sense that it is the output of a matrix product APB involving the very last transport P. This stands in stark contrast with the requirements in [4, 3. 31] that the kernel matrix corresponding to K2 C- admits favorable properties, such as being p.s.d or admitting an explicit (random or not) finite dimensional feature approximation. Since C' changes
247
249
250
251
252
at each iteration in Algo. I, they are not directly applicable.
Combining the results in $4 with those from $B results in updates for C₁ and C2 that only require 253
(nrdd) and (mrrd' operations.
254
255
Linear time GW. Finally all the quadratic operations appearing in Algorithm (3) can be replaced by linear counterparts. The iterations that have not been modified had an overall complexity of O (mr(r + d') + ar(r + d)) at each iteration. The initialization and linearization steps can now be performed in linear time, with respective complexity of respectively O(n(r + d²) + m((d'² + r)) and (((r((r + d) + mr(r + d')).
257
258
259
Experiments
280
Our goal in this section is to demonstrate that, for far smaller computational budget, the GW-LR approach is competitive with the direct entropic approach on datasets that are either synthesized to exhibit local clusters, or directly validated on a real high-dimensional dataset well. Because both approaches have different hyperparameters, our goal is to stick to a realistic evaluation that stresses both optimality of solutions as a function of computational effort, well performance in real life applications. We start by investigating the sensitivity of hyperparamaters and on our method. Since GW not convex, these may interact in unexpected ways. Experiments were run on a personal MacBook Pro 2019 laptop. We reused code from github. com//neyerscatbon//D0 and downloaded genomics data from github. cowccoo/////o wwwoo//singgaaa/SSac..
262
283
254
265
266
267
268
269
270
Benchmarks. We consider three synthetic problems and one real world problem to evaluate time- accuracy trade-offs, and also compare the couplings obtained by method and that of the entropic version [27]. More precisely, we compare the quadratic approach in GW-LR computed with algorithm (3) (and its linear time counterpat, Lin GW-LR as presented in $5), with Entropic-
272
273
274
GW, the cubic implementation of [27] (as well as its quadratic counterpart, Quad Entropic-GW presented in Algo. 2). For GW-LR and Lin GW-LR, and in all experiments, we the lower bound on entries of g to a = 10
275
278
Initialization To initialize all algorithms with common strategy, we adapted the first lower bound of [26, Def. 6.1] to the entropic case. In all experiments showing time-accuracy tradeoffs, we choose provide platform independent quantities. Accuracy is measured by evaluating the ground-truth energy E, (even in scenarios when the method uses a low rank approximation for A, B at optimization time).
278
279
to number of operations
280
7

---

100
V:50
V: v:33
10
0.100
0.300
0.100
0.300
GM loes
0.075
0.075
0.075
0.050
0.050
0.050
0.050
0.025
0.025
0.025
0.000
10
10
10
10
00000
10
30
30
16
0,000
10'
10
10"
10"
0.000
10
10
10
10
Operations
Operations
Operations
Operetions
CW-LR, (- S00
GW-LR, lle
GW-UR.
Entropic-GW
Quad Entropic-GW
500
GWAR.
56
The number of cluster in each distribution is 10 and the number of samples is n = m = 5000. The ground cost is the Euclidean distance. As can evaluate the distance between two arbitrary points, we can obtain in linear-time an efficient approximation presented in 3. The rank of their factorizations is fixed to be d d' 100. GW-LR and Entropic- GW corresponds to the case where the full matrices A and B are considered while Lin GW-LR and Quad Entropic-GW take inputs the low-rank approximations of the distance matrices. We plot the time- accurrcy tradeoff for multiple choices y and rank r defined as a fraction of n. For Entropic-GW and Quad Entropic-GW, set 1/7 proposed in [27]. Recall that for low-rank methods, we set E =0. 0.
Figure 3:
the distance matrices A and as
282
Sensitivity to y and Here we aim showing the dependence in both y and e of our proposed method. In Figure 2, we compare the GW loss obtained by our algorithm when varying E and y on two mixtures. We
283
284
285
-121
286
show that when E 0, the proposed method manage
287
to consistently obtain small GW loss whatever is. By allowing e > 0, the algorithm able reach even smaller GW loss, however, the choice of e depends highly on y. Therefore in the following experiments, 0 for dependence in y and e of our method in other settings and observe similar behaviors. See Appendix D.2 for more details
GW loss
286
022
290
fix
method. We also show the
292
log Y
294
Figure 2: In this experiment, consider two mixtures of (2 and 3) Gaussians in respec- tively 5-D and 10-D, sampled as discrete mea- sures with n m = 5000 points, see more details setup in Appendix D.2. The ground cost is the squared Euclidean distance. which provides an exact low-rank factorization of the cost as presented in $ 3. Results on speed (in Appendix) are therefore obtained using Lin GW-LR. The nonnegative rank of the coupling is set to r 50 n//100 We plot the GW loss obtained by Lin GW-LR when varying e for multiple choices of 7. Both size and color have been used to quantify visually the value of the loss at that parameter pair. Oc- casional inversions are due to the nonconvex nature of the GW problem.
295
Remark 1. As shown Figure 8 in Appendix D.2, allowing >> O may also increase the speed of con- vergence of the algorithm. However choosing well e for given y must be done carefully and we prefer in the following experiments to present the performance of our method in the simplest setting where = O.
296
297
299
300
301
Synthetic low-rank problem In this experiment aim at comparing the time-accuracy tradeoff of the different methods when the underlying distribu- tions has low-rank structure. For that purpose, we consider two distributions in respectively 10-D and
303
304
305
306
15-D, where the support of each distributions is the concatenation of clusters of points, and where the eu- clidean distance between the centroids of the clusters is bigger than threshold 3. Here we set 3 = 10. Both distributions are uniform, have the same number
307
308
309
310
of clusters and the same number of points in each cluster. Some illustrations of the simulated data provided in Appendix D.3. In Figure 3, when the underlying cost the (not squared) Euclidean distance, our methods manage to consistently obtain similar accuracy that the ones obtained by
312
313
314
entropic methods, with very low rank r = n n/500, while being orders of magnitude faster. Figure 4,
we also compare the time-accuracy tradeoff's in the more favorable case where the underlying cost the squared Fuclidean distance and obtain similar results. We also show more experiments for
316
different number of clusters in Appendix D.3, leading to similar conclusions.
8

---

250
100
50
V: 30
oss
30
10
10
10 20
10 Operations
Operations
Operations
Operations
1003
GW-LR, 280
100
Entropic-GW
Figure The number of clusters in each distribution 5 and the number of samples considered here is n m 10000. The ground cost is the squared Euclidean distance. We compare Lin GW-LR and Quad Entropic-GW as we have an exact factorization the matrices and B. We plot the time-accuracy tradeoff when varying 7 for multiple choices of For Quad Entropic- we set 8 1/7 and Lin GW-LR we set
250
100
y: 5:
y: 10
0.4
0.1
0/4
0-4
08
20
800
0.
202
60
82
00
30
00 23
Entropic-GW
GWALR, = 10
GW-LR, ==
GW-LR, = 100
Figure 5: We plot, for each cells of the SNAREseq dataset, the FOSCTTM ranked in the increasing order for both GW-LR and Entropic-GW.
y:
y: 100
0.10 y:
y:
0.10
0.10
0.10
0.04
10a
10a Operations
10
10a
10t
Operations
Operotions
Operations
GW-LR,
GW-LR, r 50
GWVLR, 100
Entropic-GW
Figure 6: Plot of the time-accuracy tradeoff when varying for multiple choices of rank on the SNAREseq dataset. For Entropic-GW we set 1/7 for GW-LR, we set
0.
Experiments Single Cell Genomics Data. We reproduce the single-cell alignment experiment introduced [14]. The dataset consists in single-cell multi-omics data generated by co-assays. In that setup, the ground truth one-to-one correspondence information between cells known, and can therefore be used to benchmark GW strategies. The dataset considered is the SNAREseq [10], with m 1047. We apply the exact same pre-processing steps as proposed in [14] by computing intra-domain distance matrices A and B with a k-NN graphs based on correlations, to compute shortest path distance matrices. Note that that case, one cannot obtain directly linear time low-rank factorization of A and B using [5, 21], since the shortest path distances need to be computed Figure we compare the alignment performance through the "fraction of samples closer than the true match" (FOSCTTM) introduced in [25]. We see that both algorithm obtain similar performance. However, in Figure 5, we show that whatever the y chosen, GW-LR reaches better accuracy while being order of magnitude faster than Entropic-GW for a very small rank r 10.
319
320
321
322
323
324
325
first. Therefore we only consider the quadratic GW-LR and the cubic Entropic-GW.
326
328
329
330
331
Conclusion. While the factorization introduced in [32] held the promise to speed up classic OT, we have shown in this work that it delivers an even larger impact when applied the GW problem: Indeed, the combination of low-rank Sinkhorn factorization with-low rank cost matrices is the only
332
334
one, to our knowledge, that ensures that the linearization step the GW objective can be carried
335
with linear complexity, throughout outer iterations. This linear complexity is comparable to that of
396
the most recent OT solvers, yet still retains the appealing properties of the Entropic approach, such as
stability and convergence to meaningful solutions.

---

References
338
339
[1] Jean Alaux, Edouard Grave, Marco Cuturi, and Armand Joulin. Unsupervised hyperalignment for multilingual word embeddings. arXiv preprint 111.0 2018.
340
[2] Jason Altschuler, Jonathan Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration. arXiv preprint arXiv: 705449944, 2017.
342
343
[3] Jason Altschuler, Francis Bach, Alessandro Rudi, and Jonathan Niles-Weed Massively scalable sinkhorn distances via the nyström method, 2018.
345
[4] Jason Altschuler, Francis Bach, Alessandro Rudi, and Jonathan Weed. Approximating the quadratic transportation metric in near-linear time. arXiv preprint arXiv: 1810,1004, 2018.
346
347
[5] Ainesh Bakshi and David Woodruff. Sublinear time low-rank approximation of distance matrices, 2018
348
349
[6] Andrew J Blumberg, Mathieu Carriere, Michael A Mandell, Raul Rabadan, and Soledad Villar. Mrec: a fast and versatile framework for aligning and matching point clouds with applications to single cell molecular data, arXiv preprint arXiv:2001 01666,, 2020.
350
351
352
[7] Charlotte Bunne, David Alvarez-Melis, Andreas Krause, and Stefanie Jegelka, Learning generative models across incomparable spaces. arXiv preprint 1eXivv1905.00 05461, 2019.
353
354
[8] Rainer E Burkard, Eranda Cela, Panos M Pardalos, and Leonidas S Pitsoulis. The quadratic assignment problem. In Handbook of combinatorial optimization, pages 1713-1809. Springer, 1998.
355
356
[9]
Lactitia Chapel, Mokhtar Alaya, and Gilles Gasso. Partial optimal transport with applications on positive-unlabelld learning. In Advances in Neural Information Processing Systems 33 (NeurlPS 2020), 2020.
358
359
[10] Song Chen, Blue B Lake, and Kun Zhang. High-throughhut sequencing of the transcriptome and chromatin accessibility in the same cell. Nature biotechnology. 37(12): 1452-11774, 2019.
381
382
Samir Chowdhury, David Miller, and Tom Needham. Quantized coaoovwwasserstein arXiv preprint arXivv2104.0201. 2021.
364
12]
Alexis Conneau, Guillaume Lample, Marc' Aurelio Ranzato, Ludovic Denoyer, and Hervé Jégou. Word translation without parallel data. arXiv preprint 17Xxiv:771.044000 2017.
365
366
[13] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in neural information processing systems, pages 2292-2300, 2013.
367
388
[14] Pinar Demetci, Rebecca Santorella, Bjöm Sandstede, William Stafford Noble, and Ritambhara
389
Singh. Gromov-wasserstin optimal transport to align single-cell multi-omies data. bioRxiv, 2020. doi: 10. /2..1//.....22.....
370
371
[15] Pinar Demetci, Rebecca Santorella, Bjorn Sandstede, William Stafford Noble, and Ritambhara Singh. wrssev-wasserstein optimal transport align single-cell multi-omics data. BioRxiv, 2020.
372
373
374
[16] Richard Mansfield Dudley et al. Weak convergence of probabilities on nonseparable metric spaces and empirical measures on euclidean spaces. Illinois Journal of Mathematics, 10(1): 109-126, 1966.
375
376
377
[17] Danielle Ezuz, Justin Solomon, Vladimir G Kim, and Mirela Ben-Chen. Gwenn: A metric alignment layer for deep shape analysis. In Computer Graphics Forum, volume 36, pages 49-57. Wiley Online Library, 2017.
378
379
380
[18] Aden Forrow, Jan- Christian Hütter, Mor Nitzan, Philippe Rigollet, Geoffrey Schiebinger, and Jonathan Weed. Statistical optimal transport via factored couplings, 2018.
381
[19] Aden Forrow, Jan Aan-Chrittian Hütter, Mor Nitzan, Philippe Rigollet, Geoffrey Schiebinger, and Jonathan Weed. Statistical optimal transport via factored couplings. In The 22nd International Conference Artificial Intelligence and Statistics, pages
383
PMLR, 2019.
[20] Steven Gold and Anand Rangarajan. Softassign versus softmax: Benchmarks in combinatorial optimization. Advances in neural information processing systems, pages 626-632, 1996.
386
387
[21] Piotr Indyk, Ali Vakilian, Tal Wagner, and David Woodruff. Sample-optimal low-rank approxi- mation of distance matrices, 2019.
10

---

389
[22]
Hicham Janati, Thomas Bazeille, Bertrand Thirion, Marco Cuturi, and Alexandre Gramfort. Multi-subject meg/ceg source imaging with sparse multi-task regression. Neurolmage, page 1 16887, 2020.
390
[23]
Hiroshi Konno. Maximization of convex quadratic function under linear constraints. Mathe- matical programming, 11(1):11-1127, 1976.
393
394
[24]
Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Q Weinberger. From word embeddings to document distances. In Proc. of the 32nd Intern. Conf. on Machine Learning, pages 957-966, 2015.
395
387
[5 Jie Liu, Yuanhao Huang, Ritambhara Singh, Jean-Philippe Vert, and William Stafford Noble. Jointly embedding multiple single-cell omics measurements. BioRxiv, page 644310, 2019.
398
[26]
Facundo Mémoli. Gromov- distances and the metric approach to object matching. Foundations of computational mathematics, 11(4):417-487, 2011.
400
401
[27]
Gabriel Peyré, Marco Cuturi, and Justin Solomon. Gromov-wasserstein averaging of kernel and distance matrices. In International Conference on Machine Learning, pages 2664-2672, 2016.
402
403
[28]
Gabriel Peyré and Marco Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 11(5-6), 2019. ISSN 1935-8245.
405
29
Yossi Rubner, Carlo Tomasi, and Leonidas J. Guibas. The earth mover's distance metric for image retrieval. International Journal of Computer Vision, 40(2)999-111, November 2000.
408
[30] Ryoma Sato, Marco Cuturi, Makoto Yamada, and Hisashi Kashima. Fast and robust comparison of probability measures in heterogeneous spaces. arXiv preprint arXiv:2002.0161 2020.
408
409
[3]] Meyer Scetbon and Marco Cuturi. Linear time sinkhorn divergences using positive features, 2020.
410
411
[32] Meyer Scetbon, Marco Cuturi, and Gabriel Peyré. Low-rank sinkhorn factorization, 2021.
412
[33] Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al. of single-cell gene expression identifies developmental trajectories in reprogramming. Cell, 176 (4):928 01188994,.211 2019.
413
analysis
414
415
416
[34] Richard Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matrices. Ann. Math. Statist., 35:876-879, 1964.
418
[35] Justin Solomon, Gabriel Peyré, Vladimir G Kim, and Suvrit Sra. Entropic metric alignment for correspondence problems. ACM Transactions on Graphics (TOG), 35(4):1-13, 2016.
419
[66]
KhrloTh Sturm. The space of spaces: curvature bounds and gradient flows on the space of metric measure spaces. arXiv preprint arXiv:12008 0434, 2012.
421
[37] Titouan Vayer, Laetita Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty. Fused gromov-wasserstein distance for structured objects: theoretical foundations and mathematical properties. arXiv preprint arXiv: /. 002834, 2018 8.
423
[38] Titouan Vayer, Rémi Flamary, Romain Tavenard, Laetitia Chapel, and Nicolas Courty. Sliced gramov-wasserstenn. arXiv preprint
428
2019.
[39] Jonathan Weed and Francis Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in wasserstein distance. Bernoulli, :2200-2262022488 2019.
428
429
[40] Hongteng Xu, Dixin Luo, and Lawrence Carin. Scalable grosov-wasserstein learning for graph partitioning and matching. arXiv preprini 1rXiv:1905.07660 2019.
430
431
[41] Hongteng Xu, Dixin Luo, Hongyuan Zha, and Lawrence Carin Duke. Grosov-wassersteinn learning for graph matching and node embedding. In International conference machine learning, pages 6932-6941. PMLR, 2019.
433
Checklist
1. For all authors...
11

---

438
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's s contributions and scope? [Yes] The main claim of the paper is a linear time (w.r.t. sample size, as commonly understood in OT) computation for GW. Sections $3,4 and 5 build up that answer. This is experimentally validated across several experiments, both synthetic, to help the reader form intuitions, and on real data where GW was deemed useful.
437
438
439
440
441
442
(b) Did you describe the limitations of your York?[Yes] Because the GW problem is non- convex, these limitations are naturally discussed in the experimental section, Section $6. We discuss the effects of 7 and E on the method, which are not casy to parse due to the non-convexity of the method.
443
444
445
446
(c) Did you discuss any potential negative societal impacts of your Work?[N/A| As a purely methodological paper, we do not envision potentially negative impact of this work on its own.
448
449
(d) Have you read the ethics review guidelines and ensured that your paper conforms to [hem?[Yes] We have read these guidelines and confirm our paper does conform to them.
450
451
452
2. If you are including theoretical results...
453
(a) Did you state the full set of assumptions of all theoretical results?[ Yes] Our theoretical results of two nature: a bound in Proposition 1, and guarantee convergence in Proposition 2, which has no direct consequence on our empirical findings, yet remain useful. Theory is not our main contribution, but rather algorithms.
455
456
457
(b) Did you include complete proofs of all theoretical results? Yes] all proofs are in the appendix. Space in the main body of the paper was prioritized to include experimental validation, which, for this non-convex problem, we believe to be equally important.
458
459
460
3. If you ran experiments...
461
(a) Did you include the code, data, and instructions needed reproduce the main ex- perimental results (either in the supplemental material or as a URL)?[Yes] We have included portions of the code that we have used. We pledge to make the entire code available later in the reviewing process.
462
463
484
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?[Yes Yes] The main contribution is compute efficiency, that we have considered across several parameter choices.
466
468
(c) Did you report error bars (c.g., with respect the random seed after running ex- periments multiple [imes)?[N/A] Our experiments are deterministic, since we use a predefined initialization.
489
470
471
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)?[Yes) The experiments were run with basic computational means, a macbook pro.
473
474
4. If you are using existing assets (c.g., code, data, models) or curating/relcusing new assets...
475
(a) If your work uses existing assets, did you cite the creators?[[es)] We have reused a single toolbox, and accessed data available publicly.
476
(b) Did you mention the license of the assets? [Yes] All licenses open source, see supplementary.
478
479
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] Yes, code shared in the supplementary.
480
481
(d) Did you discuss whether and how consent was obtained from people whose data you' 're using/curating? [Yes]
482
483
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
484
485
5. If you used crowdsourcing or conducted research with human subjects...
486
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
487
488
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
12

---

490
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
491
13


Reviews:
[{"review_id": "joWTOcnypN", "reviewer": "Reviewer_FxAF", "summary": "The paper proposes a new approximating algorithm for Gromov-Wasserstein distance with $L_2$-norm loss that runs in a linear time regarding the number of samples. This algorithm is rooted in the combination of two low-rank techniques: one is the speed-up on cost matrix computation using low-rank factorization, and the other comes from the additional low-rank constraints on the set of feasible couplings. They also provide a guarantee on the convergence to a stationary point of the algorithm. Furthermore, they experimentally show that the proposed algorithm can accurately compute GW distance compared to previous techniques while being significantly faster.", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "--- Strength:\n\nThe paper is well-written and easy to follow. The complexity notes on each algorithm make it easy to see which parts to improve next.\n\n--- Weakness:\n\n(1) The contribution is quite limited. To achieve linear computation, the author only combines existing techniques to improve different parts of an existing algorithm. While \\textbf{Quad Entropic-GW} is just a result of matrix computation speed-up under low-rank assumption, \\textbf{Lin GW-LR} is a straightforward extension of \\cite{scetbon2021lowrank}.\n\n(2) In experiments, $\\varepsilon = 0$ for \\textbf{Lin GW-LR} and \\textbf{GW-LR}, and $\\varepsilon = 1 / \\gamma$ for others. It means that the objectives they compute are different: while \\textbf{Lin GW-LR} and \\textbf{GW-LR} compute the optimal unregularized GW objective with an additional constraint, \\textbf{Quad Entropic-GW} and \\textbf{Entropic-GW} compute the optimal regularized GW objective. Is it fair to put them in comparison, especially when $\\gamma$ is small (i.e, $\\varepsilon = 1 / \\gamma$ is large)?\n\n(3) Do the authors have insights into the optimal solution of $(5)$? What effects does the low-rank constraint have on the optimal GW solution (and consequently, the optimal correspondence between two marginals)?\n\n--- Typos:\n\n(1) The notation $A^{\\odot2}$, which I believe to be $A \\circ A$, i.e, Hadamard product, is not defined.\n\n(2) Line 176 : satisfies $\\to$ satisfy, $\\alpha$ sufficient small $\\to$ $\\alpha$ being sufficient small.\n    \n(3) Double citation ([14], [15]); ([18], [19]).\n\n--- References:\n\nThe references with the complexity of approximating optimal transport are very lacking. The authors need to include several more related references. Here are a few examples: [1] is for complexity of Greenkhorn algorithm, [2] is for the complexity of the acceleration of Sinkhorn and Greenkhorn algorithms, [3] and [4] are for the optimal complexity of approximating OT.\n\n[1] T. Lin, N. Ho, M. I. Jordan. On efficient optimal transport: an analysis of greedy and accelerated mirror descent algorithms. ICML, 2019.\n\n[2] T. Lin, N. Ho, M. I. Jordan. On the efficiency of the Sinkhorn and Greenkhorn algorithms and their acceleration for optimal transport. Arxiv preprint Arxiv: 1906.01437, 2020.\n\n[3] J. Blanchet, A. Jambulapati, C. Kent, and A. Sidford. Towards optimal running times for optimal transport. ArXiv preprint Arxiv: 1810.07717, 2018\n\n[4] K. Quanrud. Approximating optimal transport with linear programs. In SOSA, pages 61–69,2019.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I do not see any foreseeable limitations and negative societal impact of this work.", "review_text": "", "extra_scores": {"ethical_concerns": "I do not detect any ethical issues with the paper.", "needs_ethics_review": "No", "ethics_review_area": ["I don’t know"], "time_spent_reviewing": "4 hours"}, "extra_text": {}}, {"review_id": "iiDx5_o1nf4", "reviewer": "Reviewer_Pn9k", "summary": "This paper presents mainly two contributions around the computability/approximation of the Gromov-Wasserstein distance (NP-hard in general). The first contribution is to consider cost matrices that admit a low-rank approximation. This allows roughly to go from a cubic problem in complexity to a quadratic problem. The second contribution is to propose to optimize the couplings of the Gromov problem under low-rank constraints. Putting end to end these contributions allow to compute an approximation of GW in $O(n)$ (with respect to the number of samples).\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "- Pros:\n\t- The article is didactic and very well written\n\t- The problem addressed is important and far from trivial: the proposed method is interesting\n- Cons:\n\t- The experiments are restricted to small toys which is a shame given that the overall message of the article is that it allows to scale GW. \n\t- The theoretical and algorithmic contributions of this article are quite incremental: there is, in my opinion, a significant overlap between this article and [1] (ref [32] in the original article) \n\t- Some experimental points seem to me rather curious.\n\nOverall the article is well written, pleasant to follow and didactic (for example Figure 1  well illustrates the low-rank approach). I thank the authors for that. The proposed approach is interesting and the fact that it allows to compute a GW in $O(n)$ is of a certain interest. Below are my various comments:\n\n- On the factorization of the costs by low-rank matrices\n\nThe idea of exploiting the low rank properties of cost matrices to go from a cubic to a quadratic cost evaluation is interesting. The first contribution (\"Exact factorization of cost matrices\") is more of an \"observation\": it proves that when the pairwise distance matrices are low-rank (as for a squared euclidean distance) then one can implement the Gromov cost in $O(n^2)$ instead of the \"naive\" cost $O(n^{3})$. However this is not true in the general case and requires a strong prior on the cost matrices. Consequently the authors propose an algorithm to approximate these distances in product of low-rank matrices, based on [5,21] (refs of the article).\n\nThis point, according to me, is one of the most interesting points of the article, in terms of applicability. Because it allows to define a general \"low-rank\" solver for Gromov with a quadratic evaluation cost. I think that this part would really deserve more details. Indeed, as it is, it is directly borrowed from [1, paragraph General Case: Distance Matrix] for the Wasserstein context and is almost identically rewritten. It would have been very interesting, for example, to understand whether, in the Gromov setting, one can quantify the approximation of GW with the low rank cost matrices to that of GW with the initial cost matrices. \n\nCan we prove, for example, that the approximated GW is approximately the true GW plus an error term (which depends, for example, on the $d$ arrival dimension and the singular values of the distance matrices)? In practice, is this approximation reasonable? For example, can we have an idea of the trade-off between the $d$ dimension of arrival and the GW approximation?\n\n- On the low-rank approach on couplings: \n\nI find this part interesting in itself but rather incremental compared to [1]. The analysis of the convergence of the Mirror Descent is almost identical to the analysis done in [1] and does not bring any additional theoretical insight: the only change is the smoothness constant of the loss, for the rest the proof is very similar. It even seems to suggest that it is as simple to approximate the problem when it is quadratic as when it is linear, which is quite surprising. I think this part requires more discussions regarding the the differences between [1] and GW-LR w.r.t. the algorithmic solution. Overall it is difficult to consider this part as a strong contribution in itlself because most of the arguments of [1] are extended in a direct way to GW.\n\n- About the experiments:\n\nThe experiments presented in the article focus on answering the following question: for a given operation budget, does the (double)low-rank approach achieve a reasonable distance for Gromov compared to the entropy solution. In this sense, I find the different experiments proposed to be relevant: they show the interest of GW-LR. However, there are some points that I find a bit unclear and that seem to me rather curious. \n\nFirst of all, the authors choose $\\varepsilon= 0$ in their experiment which is quite confusing because it implies that there is, in the end, no entropic regularization. They justify the fact that for any choice of $\\varepsilon$ the loss obtained is \"small\" and that \"the choice of $\\varepsilon$ depends highly on $\\gamma$\". So this suggests that the entropy regularization is not really useful, which is surprising compared to the global message of the article, or that it is difficult to tune the hyperparameters. Moreover, it seems to me in this case rather unfair to compare with the entropic GW with a choice $\\varepsilon = 1/\\gamma$ and not to show also this setting for GW-LR. In practice, is it interesting to consider for example the $\\gamma$ defined by Proposition 2 and to leave tune only the $\\varepsilon$ ? \n\nI think it would also have been interesting to see how the method performs against non-regularized GW solvers such as a conditional gradient [2] (since the setting chosen in the article is $\\varepsilon=0$)\n\nMore generally, the proposed method introduces many hyperparameters: $\\varepsilon,\\gamma,r,\\alpha$ and, given that it already seems difficult to tune $\\varepsilon,\\gamma$, in view of the previous discussion, I am quite skeptical about the applicability for concrete and big data sets. Especially since the experiments focus on toy datasets and do not show any application of (Lin)-GW-LR to concrete application cases. It would have been interesting to use this approach for analysis on large graphs or shape-matching where GW is not applicable because of the size of the data. As it is, I find that the experiments are insufficient to show the usefulness of (Lin)-GW-LR on real data sets. \n\nIn conclusion, the article tackles an important problem and offers interesting solutions. However, considering that the experiments are not really convincing and that the article is very similar to [1] I would tend to recommend its rejection.\n\n- Small remark:\n\t- It is not really clear to me why the low-rank procedure on the cost matrices cannot be used for the Single Cell experiment. After calculating the shortest-path why can't we precalculate the low-rank matrices and compute GW-LR this way?\n\t- Typo: (row 64) \"sped\" -> \"speed\"\n\n[1] Meyer Scetbon, Marco Cuturi, and Gabriel Peyré. Low-rank Sinkhorn factorization, ICML 2021.\n\n[2] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, Titouan Vayer. POT: Python Optimal Transport. JMLR 2021. \n\n\n---- After Rebuttal ----\n\nI would like to thank the authors for their response. I am increasing my score to 6. I think this article proposes an interesting contribution to a difficult problem. I acknowledge, in terms of novelty, that the regularization scheme based on [32] should be interpreted similarly as adding the entropic reg to the GW problem. Yet, I think it would be valuable to specialize a little bit more the low rank approach to the GW setting. Even if the authors think that this is a minor result I believe that comparing GW-LR to GW in term of objective function could be interesting.\n\nMoreover, I still think that the experiments should be improved and that, as it is, they are quite limited to toys: the low-rank regularization is proved to be useful on real data only for the single-cell genomics data (where $n,m$ are quite small) where it is only compared to entropic GW. I believe that more ambitious experiments should be done on real datasets in order to show the advantages of low-rank regularization for GW (large graphs matching/classification, 3D mesh comparison). Finally, I think that the double regularization should be more discussed and applied (it is not really used in the experiments since $\\epsilon = 0$ is the default setting). As it is I think the paper does not really show the interest of adding the entropic regularization over the low-rank constraints. \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Authors did not provide any potential negative societal impact of their work however, it is difficult for me to foresee these impacts as the work is mostly an optimization paper.\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "bS-HsZB4sRg", "reviewer": "Reviewer_nPX5", "summary": "In general, computing GW distance exactly is NP-hard. To compute nonconvex approximations of it generally take time $O(n^3)$-$O(n^4)$, depending on if the loss is separable. This work follows a recent trend in works that use low-rank approximations to cost matrices and couplings that have been used in the optimal transport literature to develop a method that has linear scaling in n, the number of points. Experiments demonstrate the efficiency and accuracy of the proposed method on some baseline experiments. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Upsides: \n\n+ The method successfully transfers low-rank ideas to the Gromov-Wasserstein setting. \n+ The method shows how low-rank approximations make all operations efficient, including the outer norms and linearization of the GW objective in Section 5.\n+ All of the steps and approximations of the methods are clearly laid out and discussed.\n+ This work achieves the fastest approximation of GW distance to date in terms of scaling with $n$. \n+ Despite the nonconvexity, the method performs well on a suite of examples, illustrating some sort of benign nonconvexity in practice.\n\nDownsides:\n\n- A lot of the method seems to be a combination of past low-rank approximations for each part of the pipeline, and the novelty lies in applying them all as a whole. \n- The method is built on top of nonconvex approximations that are still not guaranteed in general. It is hard to get a sense of the failure regimes of GW from this work in terms of the nonconvexity, as well as the additional nonconvexity added by the low-rank approximations. While it is a bit easier to see how the distance/kernel matrix can be well-approximate by a low-rank matrix, it is not as easy to see how the low-rank coupling constraint performs in practice.\n- It is easy to get lost in the whole methodology, and it would be helpful for the reader to have a roadmap of the method as well as the low-rank approximations in place.\n- In the experiments of Figure 1, why doesn’t error decrease with increasing rank?\n- What practical guidance is there for choosing the ranks in practice?\n\nAdditional questions:\n\n*Could one choose different ranks for each low-rank approximation, i.e., different ranks for the distance matrices and the couplings?\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Societal impacts are adequately discussed in the sense that it develops a general purpose tool, and GW distance has the potential for applications in places like NLP and single-cell genomics. Limitations due to the nonconvexity are discussed. However, failure regimes due to this do not appear in the work.", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "ZHVBKUYRAlu", "reviewer": "Reviewer_38E2", "summary": "In this paper, the authors proposed a new variant of GW distance with linear computational complexity. The proposed LinearGW-LR leverages the low-rank property of the distance matrices and the low-rank assumption imposed on the target optimal transport matrix, which combines the advantage of the quadratic entropic GW and that of the low-rank Sinkhorn factorization in [32]. In the experiments, the authors show the feasibility of the proposed method in various tasks, especially in the cases of low-dimensional data.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Pros:\n1. The paper is well-written and organized clearly. The notations and symbols are defined well. \n\n2. I like the detailed computational complexity analysis provided by the authors, which seldom appears in NeurIPS and is valuable for the community.\n\n3. Detailed analytic experiments on synthetic data are shown, which verifies the feasibility of the proposed method.\n\nCons:\n1. The overall algorithm looks like a straightforward modification of the method in [32], and the feasibility of the proposed method is based on three conditions: (1) low-dimensional data; (2) Euclidean distance matrix; (3) the optimal transport has a good low-rank approximation. Requiring all the three conditions jointly may limit the application of the method in practice.\n\n2. The experiments on real-world data are insufficient. Is it possible to apply LinearGW-LR to more practical tasks, e.g., 2D/3D point cloud alignment? (which I think is a potential application suitable for the proposed method)\n\n3. In my opinion, the main limitation of the proposed method is whether the target optimal transport can be approximated well by a low-rank matrix. It would be nice if the authors can study and visualize more optimal transports learned by the LinearGW-LR and quantitatively analyze the degradation caused by the low-rank approximation.\n\n4. As mentioned in Line 74, some methods [6, 40] apply recursive mechanisms to accelerate the computation of GW. For example, the S-GWL in [40] actually learns low-rank OT explicitly as the authors did, while it is from the viewpoint of GW barycenter. I think it is necessary to consider [6, 40] as baselines.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I don't think the proposed work has any negative societal impact.", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The focus of the submission is speeding up the computation of the Gromov-Wasserstein (GW) distance (capable of capturing the discrepancy of distributions defined on different spaces) using a variant of the Sinkhorn technique. The approach improves the original cubic complexity (w.r.t. the sample size) of GW computation to quadratic or to linear depending on the low-rank structure imposed.

Estimating discrepancy measures is an important problem of the area. Unfortunately, as assessed by the reviewers the work and the overall methodology closely follow the ideas and implementation pushed forward in [32] which significantly limits the novelty of the current paper. In addition, the submission could be further improved by a more thorough review of approximation schemes for the base optimal transport problem.

More work is required before publication.
Decision: Reject