Globally solving the Gromov-Wasserstein problem for point clouds in low dimensional Euclidean spaces
Authors: Martin Ryner
Jan Kronqvist
Johan Karlsson
Venue: NeurIPS 2023 poster (2023)
Keywords: Gromov-Wasserstein problem
QAP
Global optimization
TL;DR: We develop a framework for globally solving the Gromov-Wasserstein problem between two sets of points in low dimensional spaces, where the discrepancy is the squared Euclidean norm.

Abstract: This paper presents a framework for computing the Gromov-Wasserstein problem between two sets of points in low dimensional spaces, where the discrepancy is the squared Euclidean norm.
The Gromov-Wasserstein problem is a generalization of the optimal transport problem that finds the assignment between two sets preserving pairwise distances as much as possible. This can be used to quantify the similarity between two formations or shapes, a common problem in AI and machine learning.
The problem can be formulated as a Quadratic Assignment Problem (QAP), which is in general computationally intractable even for small problems. Our framework addresses this challenge by reformulating the QAP as an optimization problem with a low-dimensional domain, leveraging the fact that the problem can be expressed as a concave quadratic optimization problem with low rank. The method scales well with the number of points, and it can be used to find the global solution for large-scale problems with thousands of points.
We compare the computational complexity of our approach with state-of-the-art methods on synthetic problems and apply it to a near-symmetrical problem which is of particular interest in computational biology.

Globally solving the Gromoy-Wasserster problem for point clouds in low dimensional Euclidean spaces
Martin Ryner
Vironova AB, Stockholm, Sweden
Division of Numerical Analysis, Optimization and Systems Theory,
Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden mardin.rymerréviroovvv.ccc com, martinrrkkth.s
Jan Kronqvist
Division of Numerical Analysis, Optimization and Systems Theory, Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden jankrākth. se
Johan Karlsson
Division of Numerical Analysis, Optimization and Systems Theory, Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden johan. karl harlsson@math kth. se
Abstract
This paper presents framework for computing the Gromov-Wasserstie problem between two sets of points in low dimensional spaces, where the discrepancy is the squared Euclidean norm. The Gromov- Wasserstein problem is generalization the optimal transport problem that finds the assignment between two preserving pairwise distances as much as possible. This can be used to quantify the similarity between formations shapes, common problem in Al and machine learning. The problem can be formulated as a Quadratic Assignment Problem (QAP). which is in general computationally intractable even for small problems. Our framework addresses this challenge by reformulating the QAP as an optimization problem with low-dimensional domain, leveraging the fact that the problem can be expressed concave quadratic optimization problem with low rank. The method scales well with the number of points, and can be used find the global solution for large-scale problems with thousands of points. We compare the computational complexity of our approach with state-of-the-art methods on synthetic problems and apply it to near-symmetrical problem which is of particular interest in computational biology.
Introduction
Many important applications in machine learning deal with comparing sequences, images, and higher dimensional data, where the data unstructured and not directly comparable. In physics, chemistry, biology, music, and linguistics, objects with greatly different properties often appear in symmetrical variations characterized by concepts such as isomerisms, chirality, harmonies, and alternations. Understanding, and being able analyze, these types of variations can be truly critical
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

as some variations in chemicals and biologicals may be toxic or even lethal. The Gromov- Wasserstein
framework [15] 16 177. has shown to be powerful approach for comparing and matching such data, as it is invariant to translations and rotation. The Saanoov-asssstteinn framework has, for example, been successfully applied to domain adaptation [29]. graph matching [28] metric alignment (9].
single-cell alignment [7], and word embedding
The task of evaluating the Typically, the computational burden grows exponentially with the number of points describing the compared objects. fact, a
problem is in general considered be intractable.
problem can be formulated quadratic assignment problem (QAP) [14, [ 4] which is known to be NP-Hard. Naturally, there has been plenty of research on local and approximate methods for solving Gromov- Wasserstein and QAP problems 20 23 22 27
25 However, objects containing symmetries or repeated patterns are particularly challenging for local optimization methods and may lead significant errors in the matching such objects with local optimization methods may accidentally find the sub-optimal reflections and rotations. The inability to detect such phenomena can have great impact on the discovery of isomerisms and subsequently attributes of crucial importance.
estimated discrepancy
In this paper, we develop rigorous method for globally optimizing Strmov-Wasserstee problems by calculating a sequence of iteratively improving upper- and lower bounds. We consider a general class of Gromov- Wasssesseis discrepancy problems where the points, representing the objects, belong to Euclidean space. We show that such Gromov- Wassserttsnttnnnn problems can be formulated exactly as low-rank QAPs. We build upon this low-rank QAP representation to develop an algorithm that scales well with the number of points. The proposed algorithm can be characterized
so-called cutting plane method 10) where we solve sequence of relaxed problems that are iteratively strengthened by generating and accumulating valid linear inequality constraints, i.e. cutting planes. The optimum of the relaxed problem provides a valid lower bound for the optimum of the
problem in each iteration. By solving computationally cheap optimal transportation problem 18 26 6). we obtain both an upper-bound and new cutting plane to strengthen the relaxation. We prove convergence for the proposed a algorithm, and present a computational study that clearly shows the algorithm's 's efficiency and that the performance scales well with the number of points.
The main contribution of the paper can be summarized as:
We identify a general class of Gromov- Wasserstein problems, for point clouds embedded in low dimensional Euclidean spaces, that be exactly represented as concave low-rank QAP. In particular, mappings of images fits well within our framework.
We develop method for solving this class of Cromov-Wasserstee problems by solving sequence of alternating problems, which are either low-dimensional or linear.
We prove that the proposed algorithm converges to global optimal solution. The algorithm produces an optimality certificate in each iteration, in the form of upper- and lower bounds, which informs us of the potential suboptimality if the algorithm is terminated early.
We present a numerical study, showing the efficiency of the proposed algorithm by comparing to other global optimization methods. We also illustrate the importance of globally solving Gromov- Wasserstein problems on a problem in computational biology.
In Section we introduce the Gromov Wasserstein problems and how it can be written as QAP. In Section we identify class of Gromov- Wasserrteinn discrepancy problems that can be written as a concave relaxed QAPs problem, and in Section we present the main methodology and an algorithm for solving this class of problems. Finally, in Section we present numerical results and an application in computational biology.
2 The Gromov-Wassersee discrepancy problem
Let 01 ... Zn E X' and 1/1 ... Un E ) be two sets of points and consider the problem of finding an assignment N between the point sets such that the pairwise distances dx Øz) and dy (B M*(r)) are as close as possible for i, 1,..., ....,, n, where dx and dy represents a notion of distance on the sets X and V. respectively. This can be formulated as the discrete Gromov-Wasserster Wasserstein discrepancy

---

problem
1 TEP 2
Σ (d)(21,0)) b) R T. (.
min
and where the assignment и is represented by permutation matrix d and P is the set of all n x n permutation matrices. The corresponding relaxed problem, where instead I is in the set of doubly stochastic matrices, denoted by P, is often referred as the Ktomov-Wasserster problem [20] In these formulations we note that
Σ
(dx( 2dx 24 (dy dy 15)
Ca) 2(C,I,IT,,) + (Co, C2)
where C2 = [dx( 2y. )[ Cy [dy(yj My ")]r"= and (..) denotes the standard (Frobenius) inner product. Since the first and third sums independent of r, solving the discrete Gromov- Wasserstein problem (1 is the same as solving quadratic assignment problem (QAP) on a simplified Koopmans-Beckung form [5] namely as
(C2I,ICC)) + (CC, Cx) + (Cy,
min PEP
This problem is in general NP-hard, and the number of variables scales with the number of data points, making (2) computationally intractable for problems of relevant size. Here, we focus on instances where the matrices C., C. are positive definite and low rank. By utilizing this structure, we develop an algorithm that guaranteed find globally optimal solution and scales well with the number of points.
3 The Gromov-Wasserster problem and low rank QAP
is when the point clouds belong to the Euclidean space and the squared Euclidean distance is used as discrepancy. That is, when the set of points are 21 .1....2. Zn E & and y/1 ..... Yn R'x, which we represent by the matrices
An important special case of the in problem, considered in 20 22
X (Z1, D2 ................................
Inn E &
Y= (V1, U2 ... Un) R
In this case it can be noted that the distance matrices C₂ and Cy has rank bounded by l, + 2 and ty + 2, respectively, which can be seen from the identities
C₂ = (1|æ; - 2 = 1m -x X + 1
(3a)
Cy = (lı - 2 Img T + my
(3b)
1,),, where m, (||=,11 || 11,,,.....1..... my (1 2 11 12.... 1 and 1 E Rex1 is a column vector of ones. This observation was also used in [22] for formulating the Gromov- Wasserstein problem quadratic problem of rank (l, + (), + 2) and developing fast algorithms for the problem. However, the rank can be even further reduced and the corresponding Gromov- Wasserstein problem can be formulated QAP problem of rank E,l, 22 Lemma 4.2.3] (cf. 20 Proposition 1]).
Proposition 1. (25) Lemma 4.2.3] Let I be doubly stochastic matrix and the matrices Cz and C2 given by 3 then it holds that
(C2T, rC,) (2XTYT (XXTYT,XXTT)) + (I,T) 21 my 1T T'mz:
where L 2m 4m, 1m Ty 4XTXxmm
Proof. The proposition follows by straightforward computations, where the expressions are simplified using 1 = I'1 1T See the appendix for the proof.
D

---

Hence, the low rank QAP formulation of the discrete Gromov- Waanoo--asssesteiin problem can be stated as
min TEP
-
-2XTYY, 2 (L, +
(4)
where L = 2mmm - 4m 1 - 4X X1m and cn = ((C,, Cz) + (Cy, C₁))
41 my m/)22 The Gromov- Wasserteens problem can also be rewritten similarly and formu- lated as
min TEP
-2XIY", 2XTY?) (L,1) co.
(5)
and since the objective function is concave, any optimal solution of (4) also an optimal solution of the relaxed problem.
Proposition 2. Any optimal solution of the discrete Cromov-Wassersterr problem (4). is also
optimal solution to the Gromov- Wasserstein problem (5). Conversely, problem (5) always has an
optimal solution in extreme point and any optimal extreme point to (5) is also optimal solution to
Proof. Since (5) is the minimization of concave objective function over convex sets P. it attains the optimal value in an extreme point of the feasible set. Since the permutation matrices are the extreme points to the doubly stochastic matrices, i.e., P ext (P), (5) attains its minimum on P. Further, the set of points in P for which (5) attains its minimum are the optimal solutions of (4) To show the converse statement, note that minimum exists since P is compact and the objective function is continuous. Further, since the objective function is concave, an optimum must be at an extreme point. Finally, since the extreme points of P is the permutation matrices P, any optimal extreme point of 5 is also feasible and optimal to
In the next section we will propose a methodology and an algorithm for solving this problem.
4 A cutting plane algorithm utilizing the low rank structure
By Proposition 2, we know that an optimal solution to the discrete Gromov- Wasserstein problem (4) can be obtained by solving the relaxed problem ( However, the relaxed problem (5) is still high-dimensional non-convex QP, which is NP-hard [1]]. The high dimensionality can, in particular, be a limiting factor in solving the problem. For example, it is known that the performance of spatial branch-and-bound, one of the main approaches for globally optimizing nonconvex problems [10), can scale poorly with the number of variables. Thus, directly optimizing either ( or (5) by spatial tranch-and-bounh is not computationally tractable for larger instances. Our idea is to use the low-rank formulation of the in problem and perform the optimization a projected subspace of dimension Ey + 1 by solving sequence of relaxed problems.
First, note that problem can be written as
min WE
-| - w w-+. Co
(6a)
subject to
W = 2XTY w= ()
(6b)
Equivalence of problems and 66) is shown by simply inserting the expressions for w and W into the objective function. Next, we project out the I' variables, and we define the feasible set in the
(w, W)-space as F = Projw. (W E xt, w E IR, r E P W = 2XTYT, w = (L,
Constructing H-representation of the polytope F, i... representing it by linear constraints of the form (Z., W) + a,w S B,, not trivial and the number constraints can grow exponentially with the number of data points. Therefore, we propose an algorithm based on a cutting plane scheme to optimize over F.
An extreme point of convex set is point in the set which does not lie in any open line segment joining points of the set.

---

Instead of directly optimizing the objective in (a over the feasible set F, which we don't 't have tractable representation for, we relax the problem as
min WVRR>xYy,weE
[ +
(7a)
subject to (Z, W) + a, S, for r r 1.,,,..,N
(7b)
The linear constraints (7b) are supporting hyperplanes of the feasible set F, which we will generate iteratively. The goal is to force the minimizer of problem 7 into the feasible set F by using relatively few linear constraints. Keep in mind, we don't 't need a full representation of set F, we only need to capture the shape of F in some areas of interest, e.g., the constraints defining the faces of F the optimal solution of problem (6) would suffice. The main advantage of the relaxation in problem (7) is that it contains far fewer variables than both problems ( and (6) and the dimensionality is independent of the number of data points, Problem 7 can, therefore, be solved much more efficiently, especially in early iterations when the number of constraints low. We will show that the constraints can be determined, as needed, by solving optimal transport problems. Based on this, we will develop an iterative approach that sequentially solves problem (7) and adds constraint until the the solutions is the same as (
To initialize the search we determine bounding box of F and use this to define set of constraints (7b). The bounding box is determined by the (elementwise) minimum and maximum of the variables and W given by
min 2(XTYT)) NJs S W. max 2(XTY") )rs for i 1, Eai j ... rep
ly
(8a)
rEP
min T, S S m,,, TEP
(8b)
rEP
which can each be computed efficiently by solving standard optimal transport problem. Initializing the set of constraints by the bounding box ensures that 7 is well-defined and bounded.
If the minimizer of problem 77 is within F, then we can stop as the solution is optimal for (6)
Otherwise, we improve the outer approximation of F by adding new a constraint defined by 4N+1+ aN+1 and +1+ Let (wy, Wy) be the current optimal solution (7 and assume that (WN, WN) é F, then we form new constraint, so-called cutting plane, that excludes (WN, Wx) from the feasible set of (7)
We form a new constraint based on the gradient of the objective function 7ap. which is given by -2
(e..ec(W))>( [[]] - w.) 1-1,-2 -vvec2)))
By letting Øx+1 = 1 and ZN+1 -2WN. N. the hyperplane defining the new constraint will have the (negative) gradient in the optimum (WN. Wx) as normal vector. Then we select 3y+1 such that the new constraint forms supporting hyperplane of F (7b). This can be found by solving the following optimal transport problem
|N+1 ==
max WERRENG,,WERRIE
(ZN+1, W) + 1N+
max (4X* WXY L.)).
(9)
TEP
subject
W = -- (L,1"
When solving this problem, we also obtain a solution Nx which is a doubly stochastic matrix (generically also a permutation matrix), which gives an upper bound for (6 and candidate for the optimal solution. In the following subsection, we prove that the that algorithm, described in globally optimal solution.
Algorithm converges
geometrical illustration of the algorithm is given in Figure (I). For illustrative purposes, we have used one-ooneeeete mie-imennional data resulting in two-dimensional problem in the (W, w)-space. The data sets consist of 6 points each where one of the data sets has a reflective symmetry. This results in two global optima and 6! projected permutations in P. The first solution of 7 is located at one of the corners of the bounding box and marked with a "1" (the subsequent solutions are marked "2" "5"). The infeasible point "1" is excluded from the search space by cutting plane (red line, marked with an "A A"". Following the same procedure we obtain point "2", and cutting plane "B". Adding further cutting planes excludes "3" and subsequently "4", resulting in the feasible and optimal point "5". Note that the cutting planes from iterations 3 and 4 almost overlap the cutting planes "A" and "B". since the gradients in the points I and 3 are very similar (the same for points 2 and 4).
Remember, minimizing the objective 6al an outer approximation of the feasible set.
5

---

Algorithm Gromov-Wasserster problem
Input E RR'x* Y R Xm < 0
(Define point clouds and give tolerance level)
Libouud -0, and bound <
(Set lower and upper bounds) (Set initial constraints)
(Zr, Orr, 3r) for = 1...., N from 8 where =2222 2
while Uhound Libundd
do
(wo WN) < Optimal solution to bouud IVW wl. W + Co
(Solve
(Update lower bound) (Solve 9
IN & Optimal solution to
9
Ubund bouund x2XINY''11 (L,I')) +
(Update upper bound)
(x+++111,++1) (2WN, (4X WNY + +1,N))
(Calculate constraints)
N + N
(Update iteration number)
end while
1.50
permudal
(L.F)
Sequentall bvert bounss
+1.25
125
-1.70
-1.28
-1.75
3.22
3.24
Figure 1: Left image: An illustrative example the method on one-dimensional data. Right image: The area around the two global optima highlighting the sequence of optimal extreme points in the approximate cover and generation cutting planes.
4.1 Proof of convergence of Algorithm
The main result considering convergence is presented in the following theorem.
Theorem The gap between the upper bound and lower bound in Algorithm converges (if the tolerance is e 0).
Proof. Consider the Nth iteration in Algorithm (T. let (wN, Wx)) an optimal solution to 7 and
Ty is an optimal solution to ( with corresponding points (wi, W) = ((I,, N), 2XTx, YT)), in the (w, W) Wy-space. Assume that the gap in the objective function between those points
E
||
| v.
(10)
The new constraint is then defined by by W + S BN+1 where ZN+1 2WN and 8x+1 2(W N, W/ w and thus for any point (w, W) that satisfy the constraint must hold that
0 5 2(WN. W W) の

---

By substituting ú from 10 we obtain
eN N,WN, N, W W) w + || + +NN 1W12
=W - w + 2N WN - W) - -II -WNNl
Wx
WN 2(W,, WN -W)
(JWN - w/ + |Wx (1- 4||W>(1?3>'
where we in the last step have used the Cauchy-Schwarz Schwwrz inequality.
For any iteration number & with > O, we have that N Wy+*) is feasible for (Z W) BN+1, and thus the Euclidean distance between (WN, Wx) and (WN+K: Wx+*) ) is at least eN/(1 [1 + 4||W' w II?)). If the gap in the algorithm does not converge to O, then there is an e >> 0 for which eN for all N and thus the distance between any two points in the sequence {(WN, Wx)} is bounded from below by e/(1 [1++ +4 max || W E (8))). However, since the infinite sequence of points {wN, WN)1 } N belong
bounded set defined by (8). there must be convergent subsequence, which contradicts that there is positive lower bound on the distance between any two points.
From Theorem I the gap between the upper and lower bound converges to zero, and thus Algorithm converges to globally optimal solution.
4.2 Considerations when solving the relaxed problem
Problem minimizes a concave function over a convex set. Thus, the solution located in the extreme points of the convex set, i.e., the outer approximation of F. The standard approach to solve such problems is by branch and bound methods. However, the low dimension and sequential generation of constraints make it viable search among the extreme points for optimal solution.
To simplify notation, we define x²" = (w vec(W)²) E IR" where r' := Ez + 1. Then we can write (7b) on the form A& S by. Note that, by construction, none of the constraints are strongly redundant as every constraint satisfied with equality for permutation. As the constraints are added sequentially, it is actually easy to compute the new extreme points by keeping track of previous extreme points as described in the following proposition.
Proposition 3. Assume that the extreme points {ok]ak of the convex set described by Ax S are known. When adding a constraint Ayyx < by, the additional extreme points are linear combinations pairs of existing extreme points Tk, and 2k, both satisfying the same r constraints with equality and A5 < b, and A > b, so that the combination satisfies A7, (AT,, + +11 >))x,.) = by-
Proof. This done by counting the number constraints satisfied with equality. See the appendix for the proof.
Especially, in lower dimensions, e.g., with two or three dimensional data, this approach of keeping track of all extreme points and calculating new extreme points after adding a constraint can very efficient for solving problem (7. More details of this method is provided in the appendix. In the numerical results, we present results where problem (7) is solved both by this extreme point search and the spatial branch and bound method in Gurobi.
5 Numerical results
5.1 Computational efficiency
In this section we compare the time solve the problem up to an accuracy measured in relative error with different methods: Algorithm when 7 is solved with the extreme point method as described in section (.2,, Algorithm when solved using Branch & bound using Gurobi 10.0 [33, MILP1 formulation in [1] implemented in Gurobi and finally when (6) is directly solved using Gurobi. The MILPI formulation can handle a larger class of problems, but is reported to handle very few dimensions. All computations were performed using Matlab on an Intel i5 2.9 GHz PCC The linear optimal mass problem ( was solved using the package which is based on the network simplex
[18]. The model problems tested are evenly distributed points in unit disc or ball which we denote
7

---

Table I: Computational efficiency. Computational time on the format [mean (low high)] from repeats for various problem geometries, dimensions and sizes, for the proposed extreme point method, the same method using branch and bound (B&B), the problem formulation MILPI [1]] and finally ( implemented in Gurobi via Matlab interface. An "_" indicates that the problem timed out, such way being incomparable with the proposed method. An "!" indicates that the problem reached 10*
iterations and stopped to the accuracy indicated.
Type
Las E
Rel. emror
Algorithm [s] Extreme point B&B
MILPI [s]
B&B [s]
u
10
2,2
10
0.14 (0.07-0.3) 21 (6-47)
(11-58)
0.15 (0.14-0.16)
100
2,2
10-s
0.48 (0.3-0.7) 86 (52-107)
25 (19-39)
U
500
2,2
10
11 (9-16) / 408 (269-511)
U
1000
2.2
10
69 (54-85) 576 (389-1059)
W
2000
2,2
10
460 (313-653)/-
U
10
2,3
10 10-8 8 1.8 (1.2-2.4) 133 (45-296)
105 (49-147)
2.4(L.8-3.4)
W
100
2,3
10 8
278 (99-813)/-
172 (133-221)
U
500
2,3
10
9568 /-
M
10
2,3
10 10-8
0.51 (0.39-0.65) 708 (233-1184)
146 (66-227)
3 (2.6-4.0)
M
100
2,3
10-8
86 (20-275)
95 (73-116)
NG
500
2,3
10-5
53101/
N2
10
3,3
10 10-2
1.8 (0.7-3.2) 142 (73-210)
117 (71-163) 0.2(0.1-0.3)
Nº
100
3,3
10-2 2
36 (22-55)/
45(36-65)
N₂ 500
3,3
10-2
436 (228-862)
NG
3,3
10-2
1.2 (0.5-2.3) / 22 (11-43)
(43-94)
0.2 1-0.3) 10 (9-12)
N3
100
3,3
10-2
(5-8) 91 (76-111)
NS
500
3.3
10-2
11 (9-16)/ 161 (104-226)
N3. 1000
3,3
10-2
25 (22-29)/ 176 (149-224)
N3
2000
3,3
10-2
93 (91-100) / 578 (429-691)
U, and normally distributed points o0, We denote M1 := 1), = di,,di 1, 1)) and N3 (O, diag(1, 1, )). See Table | for numerical results. Some notes the results
1. On 2-dimensional data (fz = ly 2), the extreme point method is particularly efficient.
2. For problems that need many extreme points (> 10%), which depends on the data itself, the handling of extreme points becomes the driver of computational cost.
3. Problems mainly containing reflections (e.g. A23) easier to solve than those with room for rotations.
4. Directly solving
(6) with Gurobi was not feasible for problems with n 2 500.
5.2 Comparison with local search method
We compare the results using the proposed method local search method provided the Git-hub repository for 200 The entropy regularization parameter is set to O and the method is run with random initializations (including the first lower bound [17] used with success in (22) until the relative eиror the global optimum is less than specific tolerance The problems are the same as the previous section. Note here that in the local search methods we need an oracle in order to determine when we have reached given performance level (which of course is not available in practice), whereas we the proposed method computes upper and lower bounds.
Results for two and three dimensional data (f, Ly = 2, 3) are presented in Table 2 The results show that the proposed method performs better than multi-starting the local method when Eg- = ly = For the matching of 2-dimensional data to 3-dimensional data, the local search method is surprisingly fast suggesting that the problems is of a completely different nature than when 2-d data matched to 2-d data or 3-d data is matched to 3-d data.

---

Note that we have chosen to compare with the method from [20] rather than (22] This is since we have not optimized the optimal transport computations using the low rank structures the problem. If we optimized the computations in this way, i.e., as in 21 we expect to get similar improvement as in 22 compared with 20
Table 2: Computational efficiency compared with local search [20]. Computational time on the format [mean (low high)] on two specific problems which contain near symmetries and the number of random initializations needed to achieve the required accuracy on the format [mean (low high)]. The number of initializations were limited to 1000.
Rel.
Algorithm 1 Exec. time [s]
Local search
Type
π
Eg Ey
error
Exec. time [s]
Initializations
Sucessful runs
u
100
2,2
10-6
0.5 (0.4-0.6)
4 (0.3-12.7)
64 (4-205)
5
U
200
2,2
10-6
1.6 (1.0-199
27 (13-42)
129 (59-197)
5
U
300
2,2
10-6
3.2 (2.7-3.9)
174 (33-458)
366 (69-959)
U
400
2,2
10-6
6 (5-8)
322 (97-685)
321 (97-781)
3
U 100
2,3
10-6 6
352 (99-669)
23 (8-60)
341 (110-853)
5
200
2,3
10-6
1238 (643-2612)
92 (11-168)
421(49-778)
3
U
300
2.3
10
3006 (601-4908)
131 (4-344)
270 (10-710)
5
U
400
2,3
10-6
4729 (4279-4868)
367 (5-343)
365 (1-859)
5
N1
100
2,2
10-6
0.5 (0.3-0.7)
0,6 (0.2-1.0)
10 (3-16)
5
N
200
2,2
10-6 6
0.9 (0.7-1.1)
13.4 (1-37)
61 (6-168)
5
N
300
2,2
10-6
2.4 (2.0-2.9)
37.5 (0.5-90)
75 (1-182)
5
N
400
2,2
10-6
4.4 (3.6-5.3)
149 (22-378)
142 (21-361)
4
5.3
Application to symmetrical data for morphological analysis
In this example investigate the impact of correctly evaluating the Gromov- Wasserstenn discrepancy compared to estimating it by local search. As a test case, we examine the ability to classify Adeno Associated Viral (AAV) particles based on the Ktomov--asserseee originating from transmission electron microscopy. hence & ly 2. AAV particles are nearly round viral particles with multiple near-rotational symmetries as illustrated in Figure By sampling 500 positions on each AAV particle proportional to the protein density, the point sets from pairs of particles X; and X, can subsequently be compared using the Gromov- Wasserstein discrepancy.
discrepancy on image data
Computing the tedious. Therefore, one may consider calculating the discrepancy to a subset of the objects that well distributed under the actually calculating all pairs of discrepancies, we use a greedy approach by defining the index subset Se [s:)+= by selecting the first object arbitrarily and then let the set grow by
discrepancy between all objects in large set X' {
discrepancy. To find such subset without
Sk+1 == {Sk, argmax min dew Xy).
(11)
VgS JES,
In this way all objects are closer than tolerance to object in the subset, and the way the subset is produced generates a monotonically decreasing tolerance. By using this procedure, every object obtains feature vector of distances to the objects indexed by Sg. Next the feature vectors used as input to a K-means clustering and classification quality in terms of purity, adjusted rand index and normalized mutual information compared an expert evaluation is presented in Figure
The example shows that the data contains symmetries, local search methods may get stuck on permutations that are locally optimal but that are far from globally optimal, Figure 2 Thus the discrepancy obtained from local search is inappropriate to use as a discriminating feature, as shown in the lower performance in the classification illustrated in Figure In the left subfigure of Figure 2 self-similarities are visited the proposed method, and these local optima are fact also almost optimal when used as initiation point using local search methods e.g., (20)
This example shows that when the Wasserennsse problem calculated accurately it provides valuable information and biological meaning it differentiates viral particles with different cargo and

---

variations in capsid structure, and, at the same time, finds the optimal orientation positively revealing possible chiralities or isomerisms. It also shows that when the distance calculated accurately, provides better decision support than using local search methods.
5.16% minror
8.77% mirror
0%, rotation
8-72%,
(Moud
Figure 2: Left: Each distance calculated with the proposed method and the same distance calculated with the local search method
[20]. For near symmetrical data, the confusion of the measurements of the local method is clear. Middle: The relative eтror of the local method compared to the global solution. Right: Matching of structure to itself to four different orientations visited by the algorithm where the color indicate the permutation r' The relative error of the Kromov-Wasserstein discrepancy and the isometry to the global optimum is written near each matching. The global optimum is able to correctly match itself, whereas the local search method may get stuck in local optima.
infomation
0.6
0.0
Adused Fand inoox
0.5
Putty
0.7
0.4
0.4
0.3
beal
0.8
[5.]
[SL]
The trajectory of the increased quality of classification compared expert evaluation when the distance is computed from all particles to the sequence of particles suggested in the text. The problem is unambiguous
Figure 3:
gain of quality using an exact evaluation (blue) of the over the local search (red).
6 Discussion
When using distances input for statistical analyses, the accuracy of the measurement set bound for the information resolution. If the measurement system introduces error a certain structure, this can produce artefacts in the result and affect decisions taken on the result. When using distances for such purposes, it is necessary to either know the measurement error, the artefacts being produced, or using accurate measurement system. In this paper we have provided method which computes the in problem accurately, which reduces the uncertainty of such considerations.
Acknowledgement
This work was funded by Vironova and the Swedish innovation agency Vinnova through the innovation milieu GeneNova (2021-02640).
10

---

References
[1] David Alvarez-Melis and Tommi Jaakkola. Gromov- Mamoov-asssesttein alignment of word embed- ding spaces. arXiv preprint arXiv:1809.00013 2018.
[2] Florian Beier, Robert Beinert, and Gabriele Steidl. On a linear Gromov-Wassersten distance. IEEE Transactions on Image Processing, 31:7292-7305, 2022.
[3] Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich. Displacement interpolation using lagrangian mass transport. In Proceedings of the 2011 SIGGRAPH Asia conference, pages 1-12, 2011.
[4] Rainer Burkard, Mauro Dell' Amico, and Silvano Martello. Assignment problems: revised reprint. SIAM, 2012.
[5] Rainer E Burkard, Eranda Cela, Panos M Pardalos, and Leonidas S Pitsoulis. The quadratic assignment problem. Handbook of combinatorial optimization, pages 1713-1809, Springer. 1998.
[6] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26:2292-2300, 2013.
[7] Pinar Demetci, Rebecca Santorella, Bjöm Sandstede, William Stafford Noble, and Ritambhara Singh. Gromov- Banov-Wassee Wasserstein optimal transport to align single-cell multi multi-omics data. BioRxiv, pages 2020-04, 2020.
[8] Li Deng. The MNIST database of handwritten digit images for machine learning research [best of the weh]. IEEE signal processing magazine, 29(6): 141-142, 2012.
[9] Danielle Ezuz, Justin Solomon, Vladimir G Kim, and Mirela Ben-Chen Gwenn: A metric alignment layer for deep shape analysis. Computer Graphics Forum, 36(5):99-55) 57, 2017.
[10] Christodoulos A Floudas. Deterministic global optimization: theory, methods and applications, volume 37. Springer Science & Business Media, 2013.
[11] Michael Friesen. Low rank quadratic assignment problem: Formulations and experimental analyvis. PhD) thesis, Science: Department of Mathematics, 2019.
[12] Ralph E Gomory. Solving linear programming problems in integers. Combinatorial Analysis, 10:211-215, 1960.
[13] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023.
[14] Tjalling C Koopmans and Martin Beckmann. Assignment problems and the location of economic activities. Econometrica: journal of the Econometric Society, pages 53-76, 1957.
[15] Facundo Memoli. On the use of Tromov-Hausdone distances for shape comparison. In M. Botsch, R. Pajarola, B. Chen, and M. Zwicker, editors, Eurographics Symposium on Point- Based Graphics. The Eurographics Association, 2007.
[16] Facundo Mémoli. Cromov-Hausdove distances in Euclidean spaces. In 2008 IEEE Computer Society Conference Computer Vision and Pattern Recognition Workshops, pages 1-8. IEEE, 2008.
[17] Facundo Mémoli. Cromov-Wassersterr distances and the metric approach to object matching. Foundations of computational mathematics, 11(4):417-487, 2011.
[18] J.B. Orlin. A polynomial time primal network simplex algorithm for minimum cost flows. Mathematical Programming, 78(2):109-129, 1997.
[19] Panos M Pardalos and Stephen A Vavasis. Quadratic programming with one negative eigenvalue NP-hard. Journal of Global optimization, 1(1):15-22, 1991.
[20] Gabriel Peyré, Marco Cuturi, and Justin Solomon. Gromov- Krrmov-Wassereee Wasserstein averaging of kernel and distance matrices. In International Conference on Machine Learning, pages 2664-2672. PMLR, 2016.

---

[21] Meyer Scetbon, Marco Cuturi, and Gabriel Peyré. Low-rank Sinkhorn factorization. In International Conference on Machine Learning, pages 9344-9354 PMLR, 2021.
[22] Meyer Scetbon, Gabriel Peyré, and Marco Cuturi. Linear-time Gromov-Wassersiee Wasserstei distances using low rank couplings and costs. In International Conference on Machine Learning, pages 19347-19365. PMLR, 2022.
[23] Justin Solomon, Gabriel Peyré, Vladimir G Kim, and Suvrit Sra. Entropic metric alignment for correspondence problems. ACM Transactions Graphics (TOG), 35,))11-33, 2016.
[24] Cole Stiegler. Efficient local optimization for low-rank large-scale instances of the quadratic assignment problem. The University of Iowa, 2018.
[25] Titouan Vayer. A contribution to optimal transport on incomparable spaces. arXiv preprint arXiv:201 1.04447, 2020.
[26] C. Villani. Topics in Optimal Transportation, volume 58. Graduate studies in Mathematics, AMS, 2003.
[27] Hongteng Xu, Dixin Luo, and Lawrence Carin. Scalable Gromov-Wasserseee Wasserstein learning for graph partitioning and matching. Advances in neural information processing systems, 32:3052-300 20022, 2019.
[28] Hongteng Xu. Dixin Luo, Hongyuan Zha, and Lawrence Carin Duke. Gromov-Wasserva learning for graph matching and node embedding. In International conference on machine learning, pages 6932 6941. PMLR, 2019.
[29] Yuguang Yan, Wen Li, Hanrui Wu, Huaqing Min, Mingkui Tan, and Qingyao Wu. Semi- supervised optimal transport for heterogeneous domain adaptation. In ISCAI, volume 7, pages 2969-2975, 2018.
12

---

A Appendix
A.1 Proof of Proposition
Proposition 4. Let D be doubly stochastic matrix and the matrices C2 and Cy given by
C2 = (1124 - - -= 1m - -2X X + m
Cy = (lly- 8z11) 1m. Ty + my 17
then it holds that
(C,I, = (2XTX"XXXTY") (L,T) + 21 my 1" 1 where L m 4m, 4x"X
m2
Proof. The proposition follows by the following straightforward computations, where just expand the expressions and use T1 = 1 1, thus
tr(CuTC,r") tr((1mg - 2X X m. 1 T my my r
J( =mfr1mm - T + m2 1
-2tr(X"XT(1mm) - " + 1
+1"r(1m, 2Y Y + m,1 1,11) " m.
-m "m" 2m 2mmTmmy TT + nm TT,,
-2m7r X X1 + 4 4 XIY
-1 XI'm
+um² m. - 21 yr m, + +1"m,, 1",,,mmm mz
XTX 4m ry 2mm=Imm
x XTm, + 21 my m,,1' m2
=(2XTY 2XXTT")
+ m2 4m, 1r 4x x2mm,
+ +1 my,"" m.. 1m
A.2 Proof of Proposition 3
Proposition 6. Assume that the extreme points {RR]a of the conver set described by As S b are known. When adding a constraint ATya < by, the additional extreme points are linear combinations of pairs of existing extreme points Ik, and œk, both satisfying the same T- constraints with equality and A 2k,, I bn and A Z bn so that the combination satisfies AT ()TK, + +(1 +11 (1=A)xk,) = by-
Proof. Let W be matrix whose columns consist of the extreme points defined by the N 1 constraints Ax < b. Also, let a7 < by be an additional constraint, be unit vector with on position l; and let a parametrize the convex cone on W, i.e. 1 = 1, 2 0 so that AWa S describes all points the convex set. Suppose that B₁, describes the indices of the constraints that define the k:th extreme point by letting AB, be the sub matrix of A including the rows denoted by the indices in BR. Then AB, W ex 6B.- It then follows that AW AWee,, +1 A)ex2); = b, ifand only if E Bt, and j E BR, and 0 < X < 1. The construction in the proposition, AzE, + (1 - >) so by then satisfies r constraints with equality and all other constraints with inequality, i.e. the point is an extreme point to the set. We also note that every multiple combination (three or more) of extreme points sharing r constraints are not extreme points as they are linear combination of the pairs given by a(Arx + (1 A1)=2)) + (1 - a)222 + (1 - )2a)))) Every linear combination constraints, i.c they are not extreme points.
that AT (AIR,
(1 -1)r))))))
pairs of points sharing less than r constraints will not satisfy r
13

---

A.3 A description of the implementation of the extreme point method
The handling of the extreme points in the paper is done by keeping track of the extreme points, their connection to the boundary constraints, and lookup tables for the adjacent extreme points, i... extreme points that satisfies the same 1 constraints with equality, where is the rank of the problem. Let the extreme points be described in the matrix E where each column describes an extreme point. Thus AE Σ b1, if the constraints described by the matrix A and vector b such that Ax Σ is the constraint equations.
Let the adjacency be described in (sparse) matrix with binary elements D where D.,- 1 if extreme point and j are adjacent. Let us also keep track of the constraints that are satisfied by an extreme point with equality. For this purpose let B be matrix with binary elements in which the element Bu = 1 if extreme point satisfies constraint i with equality. Thus, Dij = 1 (B"B),, = This is one of the computational drivers for the proposed extreme point method.
When new constraint (An bn.) added, the new extreme points are generated by linear combina- tion of the infeasible extreme points I := {i: A, E, > b, and their adjacent feasible extreme points f: = {j j) 1} where E Let # indicate the number of elements of a finite set, then E LIEZ # (f:) is the number of new extreme points. We place the new extreme points in the matrix E
by adding them the end as a matrix The new matrix containing the extreme points
En = (F P)
The matrix keeping track of which extreme points satisfies which constraints with equality is extended with
B 0
Bn
where ==
= B., & B. B.,efle preferably implemented using bitwise operators. Here, the last row describes the newly added constraint and k re-enumeration of the new extreme points.
The new adjacency matrix D,, be concatenated with the old D and two additional matrices
D Dn = OT N
O
where O,, = the old extreme point i adjacent to the new extreme point j. This information is already available for us, since the new extreme points adjacent fy. Finally, Ni. if (C* C)i 2. This is by far the most computationally expensive operation in the proposed algo- rithm, which can implemented with std:popcount the standard c++ library. For 3-dimensional problems, around 120 bits needs to be compared between all new extreme points.
A.4
Additional tests
A.5
Convergence rate
The convergence proof in Proposition 3 does not include rate of convergence. In Figure 4 show the convergence trajectory for the problems 24 and N₁ for 2-dimensional data. The tests show that the convergence rate is linear to its nature up to number of iterations where the gap closes completely.
In comparison the result of the proposed method we present the same convergence evaluation with the local search method in [20]. Figurd]sshws the trajectory of convergence. The tests show that rate of convergence sublinear to its nature. The number of initializations needed to achieve pre-determined accuracy increase with the number of points for the local search. Note here that in order to determine when to stop one needs to know the optimal value, which is not available for the local search method.
14

---

10P
10a
10
10%
10/
Borations The 100
Ioraaioos
m=200
n-400
n=400
10
10
10F
Reaatee era
10
10
100
100
350
Ilorations
Boraaions
Figure 4: Top The trajectory of the relative error for the proposed method for the problem 2/ 2-dimensional data. Bottom row: The trajectory of the relative error for the proposed method for the problem Wi 2-dimensional data. The solid blue line indicates the mean convergence rate for 20 runs and the dashed red line indicates standard deviation from the mean.
107
10"
Retalive
10
10/
10"1
1019
10
intistizations
Intlaiization
106
n=100
10P
n--200
n-300
109
n=400
Racaho anar
12
10
50
10
200
200
500
400 Intollnntione
infodlonoons
Intlls Lealiame
inaiskratione
Figure 5: Top row: The trajectory of the relative eπror for the local search method for the problem 24 on 2-dimensional data. Bottom row: The trajectory of the relative error for the local search method for the problem Mi on 2-dimensional data. The solid blue line indicates the mean convergence rate for 20 runs and the dashed red line indicates standard deviation from the mean.
A.6 Evaluation the MNIST dataset
In this section we present the performance of the proposed method on shapes originating from the MNIST dataset [8]. Even though the arabic numerals not entirely unique to reflections (e.g. "2" and "5") and rotations (e.g. "6" and "9"), the ability to compute the exact Gromov- Wasserstein discepancy is of interest. When the images doesn't contain the same number of points, evenly distributed points are taken from the image with more points. The maximum number points was 400. Number of iterations are presented in Table 3 and the computational time is presented in Table 4 showing that the proposed method works as anticipated on the data set. Examples of correspondances of numerals are shown in Figure 6
15

---

Table 3: Number of iterations to compute the distance between two numeral shapes to relative error gap 10 Numbers presented are median in the upper table and [min max] the lower table. The number of random tests were 10 for each combination.
0
2
5
8
9
0
92.5
55.5
98.5
75.5
115.5
106.5
80
81
90
90
21
46.5
36
70.5
50.5
47
50
49
2
52
53.5
64
58.5
60.5
57.5
58
43
40.5
52.5
72.5
46.5
56
56.5
4
84.5
57
63.5
65
73.5
61.5
53.5
64.5
49.5
6
44
45.5
67
56.5
48.5
64
47.5
8
62.5
55.5
52.5
0
2
3
4
6
7
8
9
67-134
39-80
67-223
51-139
86-184
69-131
62-97
63-143
73-102
53-163
10-39
32-78
30-64
49-98
36-61
35-88
27-69
43-59
27-93
29-77
36-78
37-75
44-88
47-95
38-72
43-86
36-83
3
29-67
31-116
47-60
42-89
27-77
48-85
27-87
50-107
40-111
40-111
52-92
47-108
50-80
5
31-74
45-88
43-87
54-94
36-97
26-80
36-88
54-123
32-94
24-76
52-75
30-87
40-82
42-71
9
29-65
Table Computational time for the distance between two numeral shapes relative error gap of 10-$. Numbers presented are [median]s in the upper table and [min max]s in the lower table. The number random tests were for each combination.
1
2
3
4
5
6
7
8
0.85
0.48
1.03
0.76
1.09
1.33
0.83
0.81
0.98
0.95
0.05
0.18
0.15
0.32
0.22
0.18
0.15
0.3
0.29
0.34
0.33
0.41
0.39
0.46
0.37
0.47
0.34
3
0.31
0.7
0.41
0.64
0.39
0.52
0.51
0.85
0.61
0.6
0.67
0.86
5
0.53
0.62
0.49
0.75
0.5
0.54
0.49
0.96
0.66
7
0.61
0.84
0.59
0.77
0.74
9
0.92
0
2
5
8
9
0
0.5-1.4
0.3-1
0.6-2.4
0.5-1.5
0.9-1.9
0.6-1.6
0.6-1.2
0.5-1.7
0.6-1.4
0.5-2.1
1
0-0.1
0.1-0.4
0.1-0.4
0.2-0.4
0.2-0.3
0.1-0.4
0.1-0.4
0.2-0.4
0.1-0.5
2
0.2-0.5
0.2-0.6
0.2-0.5
0.3-0.6
0.3-0.8
0.3-0.5
0.3-0.7
0.2-0.7
3
0.2-0.7
0.2-1.2
0.3-0.5
0.3-0.9
0.2-0.7
0.4-0.7
0.2-0.8
4
0.4-1
0.3-1
0.3-1.3
0.4-0.8
0.4-1.3
0.4-0.9
5
0.2-0.9
0.4-1
0.3-0.9
0.5-1.1
0.4-1.1
6
0.2-1
0.3-1
0.6-1.8
0.3-1.2
7
0.3-0.9
0.7-1.2
0.1-1.4
8
0.6-1.4
0.5-1.1
9
0.5-1.1
16

---

Figure 6: Left image: Matching of the numeral "6" the numeral "9". Right image: Matching of the numeral "2" to the numeral "5". The colors represent which point in one numeral corresponds to another point the other numeral.
17


Reviews:
[{"review_id": "72QKhbNnMO", "reviewer": "Reviewer_pM3y", "summary": "This paper proposes a new algorithm to solve the Gromov-Wasserstein problem between two sets of points in Euclidean spaces when the ground cost is the squared Euclidean norm. This is done by first reformulating the Gromov-Wasserstein problem as a low-rank QAP problem, \nthen relaxing the set of admissible couplings in a way that the optimal solutions of the obtained relaxed problem stays the same, and then using a cutting-plane method to solve the relaxed problem. The obtained algorithm provides at each iteration a lower bound and an upper bound on the value of the optimal solution, and it can be proven that the gap between those bounds theoretically converges to zero at infinity, implying that the proposed algorithm theoretically converges to a global optimum. The experiments offer a comparison in term of computational efficiency between the proposed algorithm and a more \"traditional\" algorithm, and emphasize the importance of converging to a global optimum when using Gromov for applied problems. ", "questions": "The following questions are related to the steps of reasoning that I find could be clearer:\n- 1/ What is the definition of the rank of a QAP problem? \n- 2/ Why do we need to relax Problem (4) to doubly stochastic matrices instead of permutation matrices before applying \nthe cutting-plane method? \n- 3/ Can you explain how to obtain the initial $ (Z_r,\\alpha_r,\\beta_r)_r $ from solving the bounding box problem? \n- 4/ I'm a bit confused about Equation (8a): are we looking for a different $ \\Gamma $ for each $ i,j$? The way it is written suggests we're looking for the same for all $ i $ and $ j $. \n- 5/ $ N $ stands for the current iteration number,  but also for the number of initial constraints. I find that confusing. \n\nExperiments:\n- 1/ In Figure 2 left: are the differents structure plotted are local optima obtained using a local search method? I'm not sure \nto understand Figure 2 right either: the blue is the result of your algorithm? can you explain why it is better than the result with the local search method and where do we see the comparison with the expert evaluation? \n- 2/ How the methods compare itself (in terms of computational efficiency) to approached methods (as for instance Entropic Gromov, sliced Gromov,etc) with several random initializations? \n", "limitations": "The authors has not clearly discussed the limations of the algorithm. \nIn my opinion, the principal limitation is that the method remains computationally costly, and so, as this is the case for almost all methods that compute Gromov, it is still unusable in higher dimension or when the number of points is large. \n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "I enjoyed reading this paper. To the best of my knowledge, the proposed algorithm is new. \nBeside the proposed algorithm in itself, the paper highlights that most of algorithms which compute Gromov only converge to a local minimum, which is a problem that is often ignored in the literature. Moreover, when solving the\nGromov-Wasserstein problem with \"traditional\" algorithms, it is not even possible to know in practice whether the solution\nwe've converged to is a local or global optimum. Thus, proposing an algorithm that is guaranteed to converge towards the global optimum is a real strength in my opinion. ", "weaknesses": "- As someone familiar with the Gromov-Wasserstein problem but not expert in optimization, I find this paper a bit difficult \nto read because some important steps in the mathematical reasoning are lacking of details and so, in my opinion, too much of the reasoning is left to reader.\n- I think the experiment section could be improved: the experiment of Figure 2 could be clearer, and it could be nice \nto have a synthetic experiment that exhibit a case where the proposed algorithm converges to the global optimum whereas the local search method remains stuck in a local optimum.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you sincerely for the comments, feedback, and questions. See below for our answers.\n\nWeaknesess:\n1. We will go over the mathematical derivations again and try to improve the readability for the final version of the paper. \n2. We are planning the expand the numerical results in the final version. We will improve the description of the experimental setup in Figure 2. In the supplementary material 1.4.3, we presented the statistical performance of the local search method on this example, which in mean has a relative error of above 10%. We will move this into the paper in the final version. We also plan to change Figure 1 so that it illustrates that one could get stuck in a second local optimum using a local method. Here we can identify the whole region of attraction for the local optimum that is not a global optimum.\n\nQuestions:\n1. A QAP problem may be written on the standard form $\\min_{\\Gamma \\in P} trace(W\\Gamma D\\Gamma^T)$. Now $ trace(W\\Gamma D\\Gamma^T) = vec(W\\Gamma D)^T vec(\\Gamma) = vec(\\Gamma) ^T (D^T\\otimes W)^T vec(\\Gamma)$. Then we denote the rank of $D^T\\otimes W$ as the rank of the QAP. In the GW case using the squared Euclidean norm we have the quadratic expression $trace (X^TX \\Gamma Y^TY \\Gamma^T)$. If we identify $W = X^TX$ and $D = Y^TY$, then the rank is the product of the dimensionality $l_x$ and $l_y$ of $X$ and $Y$ respectively, because of the eigenvalue structure of the kronecker product.\n2. We are not aware of any convenient way to represent the set of points that corresponds to permutations in the restricted low dimensional space. The main purpose of the relaxation is to obtain easier subproblems, and to the best of our knowledge, optimizing over permutation matrices is in general a difficult combinatorial problem.\n3. Equation 8 describes the bounding box problem, as they describe the min and max values of the elements in the low rank representation. Equation 8a are $\\ell_x\\ell_y$ equations indexed by $i = 1... \\ell_x$ and $j = 1...\\ell_y$ giving $Z_{i,j} = \\delta_{i,j}$, $\\alpha_{i,j} = 0$ and $\\beta_{i,j} = max_\\Gamma (2X\\Gamma Y^T)_{i,j}$ for the maximum expression for example. \nNote also that these inequalities can be written more plainly as $W^{\\min}\\le W\\le W^{\\max}$\nwhere $W^{\\min}$ and $W^{\\max}$ are specified elementwise by the minima and maxima in Equation 8a.   \n4. Equation 8a  is $\\ell_x\\ell_y$ equations with one optimal $\\Gamma$ for each $(i,j)$.\n5. N stands for the number of constraints. We will clarify this and how it relates the iteration number in the final version.\n\nExperiments:\n1. The proposed algorithm provides a sequence of permutations that may be close to locally optimal as visualized in figure 2 left. The objective value in these permutations may be very close to the global optimum (written in relative error) but suggest a permutation which is very far away from the global optimum permutation geometrically, here being rotations and mirror symmetries. This is only to exemplify the complexity of the problem to be solved. Figure 2 right shows classification quality metrics compared to the expert classification (ground truth), i.e., 1 would be a perfect result. The figure shows that the proposed method results in a better solution with regard to all the quality metrics. We will try to clarify this better in the final version. \n2. A key advantage with this method is that one knows when optimum is reached. When considering local search (e.g., Entropic Gromov) it is difficult to know how close the best solution is to the optimal solution. In particular, this is the case when considering problems with symmetries, which are common in for example computational biology. In the final version, we will more clearly highlight the advantages of our global optimization approach compared to local search with multiple initializations. The sliced Gromov-Wasserstein version handles problems when the spaces are of different nature. However having a distance matrix, under certain conditions it is possible to find a representation in Euclid coordinates.\n\n\nLimitations:\n\nIt is true that the method is computationally costly if the number of dimensions increases. However, for small dimensional problems (e.g., 2d and 3d objects which are common in practice), the method scales well in the number of points (in contrast to most other methods for solving GW problems)."}, {"review_id": "F6khLciCuY", "reviewer": "Reviewer_G7H7", "summary": "This paper solves the Gromov Wasserstein (GW) distance problem for squared Euclidean norm by considering a low-dimensional space on which the computation is performed, using a cutting-plane method. To this end, they write the GW problem as a low-rank optimization problem. Their algorithm is supported by theoretical convergence guarantees.\n", "questions": "The matrices $C_x$ and $C_y$ are assumed to be positive definite and of low rank. How can the low-rank assumption be ensured? What does this mean in terms of the point clouds $X$ and $Y$?\n\nCan the cutting plane algorithm trick be extended to the case where both point clouds have weights (i.e. a non-uniform discrete measure), and for which the probability measures do not have the same number of points (therefore the mass in the transport plan $\\Gamma$ could split)?", "limitations": "A comparison with more recent methods than [16] for computing the Gromov Wasserstein distance should be included.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This article makes it possible to compute the Gromov-Wasserstein distance efficiently. In a way, this paper extends the results of [17] which reformulated the Gromov-Wasserstein problem as a low-rank quadratic problem. The idea of using the cutting plane algorithm to efficiently solve this optimization problem in a projected low-dimensional subspace is new, and the proposed algorithm converges to a global optimal solution.", "weaknesses": "Solving problem (9) at each iteration of the algorithms remains computationally heavy, as it corresponds to solving an optimal transport problem with $\\Gamma$ a $n\\times n$ matrix.\n\nFor a more relevant contribution, and as the use of the low-rank quadratic structure of the Gromov-Wasserstein optimization problem has already been exploited in [17], the authors should compare their results with [17] in experiments, and not only with the local search method of [16].", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you sincerely for your comments and questions. See below for our answers. \n\nFirst we would like to clarify that the paper is not just an extension of the work [16] and [17]. Even though we also consider the GW-problem, our method is guaranteed to reach the global optimum for the problems we consider. This is to the best of our knowledge the first results that guarantees a global optimal solution for a class of GW-problems where the domain is multidimensional (for 1-d such results exists). Also, note that even though the work [17] uses similar low rank structures, the structures are used completely different. In [17], the low rank structures are used to speed up the method from [16]. In our paper, the low rank structure is used to ensure convergence to a global optimum (and not just a local).\n \nWeaknesses:\n\nThe proposed method considers the topic of a global optimum of the GW distance. The main point of the paper is to show that one can find a global optimal solution within a reasonable time, and we do not claim to be faster than [16] or [17] on any GW-problems. The local search algorithms converge faster to a local optimum, but the advantage of our approach is that we are guaranteed to reach a global optimum. For many applications, this is critical.\n\nAlso note that we have not focused on optimizing the computational complexity in each of the subparts of the algorithm. For example, we do not utilize Sinkhorn’s method for solving the optimal transport subproblems (equation 9). It is true that [17] utilizes low rank structures similar to the ones we consider for improving the computational complexity of [16]. Similarly, we could use the methods from [Scetbon, Cuturi, Peyré - ICML, 2021] for improving the computational speed of the subproblems in our method. But since we have not optimized these computations, we think that it makes more sense to compare with the method in [16]. We will add this comment to the final version of the paper to motivate why we compare with [16] instead of [17]\n\nQuestions:\n1. With Squared Euclidean distance, this is ensured trivially as $C_x$ and $C_y$ will have rank at most $\\ell_x+2$ and $\\ell_y+2$. For more general settings, it gets a bit more complicated but this is out of scope of this paper. Work on low rank approximations of point clouds to their distance is a whole separate field of investigation, which we can add references to.\n2. Yes, then $L = (1^T\\mu + 1^T\\nu) m_xm_y^T -4m_x\\nu ^T Y^T Y-4X^T X\\mu m_y^T$, $c_0 =(\\langle C_x,C_x\\rangle+ \\langle C_y,C_y\\rangle -4\\nu^Tm_y \\mu^Tm_x)/2$ where $\\mu$ and $\\nu$ are the marginals.\n\nLimitations:\n\nPlease see the answer in the weaknesses section.\n"}, {"review_id": "UmUO3cJLaV", "reviewer": "Reviewer_FaYe", "summary": "- The paper introduces a novel algorithm for efficiently computing the Gromov-Wasserstein (GW) distance for low-dimensional point clouds. The proposed method offers a transformative approach by transforming the GW distance problem into a sequence of concave optimization problems over convex sets, thereby improving computational efficiency.\n\n- By leveraging this new algorithm, the authors conducted comparative evaluations against existing techniques for computing the GW distance. The results demonstrate significant improvements in terms of computational speed. ", "questions": "- line 190: in the last inequality, am I right that a square root is missing?\n- line 196: Is there a typo in the lower bound?", "limitations": "N/A", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- One of the key strengths of this paper is the introduction of a novel and interesting algorithm. The proposed method offers a fresh perspective on computing the Gromov-Wasserstein (GW) distance for low-dimensional point clouds. \n- The research holds significant importance in the field of machine learning and data science, given the growing popularity of the GW distance. By proposing a method that accelerates the computation of the GW distance, the paper addresses a practical need and offers a valuable contribution to the field.\n- The clarity of presentation is a notable strength of this paper. The authors effectively communicate their ideas, methodologies, and findings, making the paper accessible to readers. \n- The experiments well support the acceleration effect of the proposed algorithms.", "weaknesses": "A potential weakness of the paper is the lack of a complexity analysis for Algorithm 1, which could be addressed to further enhance the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you sincerely for your comments, feedback, and questions. See below for our answers.\n\nWeaknesses:\n\nWe agree that this is a weakness and we are currently working on quantifying the convergence rate. For a given dimensionality and number of points, it should be straightforward to show that our algorithm finds an $\\epsilon$ accurate solution in $O((1/\\epsilon)^{\\ell_x\\times \\ell_y+1})$ iterations. However, this trivial complexity bound is not very satisfying and it doesn't take in consideration the increasing complexity of the QP problem. We believe that it should be possible to prove a faster convergence rate. In practice, we have experienced that the method converges much faster as exemplified in section 1.4.1 in the supplementary material.\n\nQuestions:\n1. Line 190: Yes, thank you for finding this error.\n2. Line 196: The notation might be confusing, what we mean is that the lower bound is limited to the largest possible $\\|W\\|_F$ obtained in the bounding box. We will improve the readability of this.\n"}, {"review_id": "bpUTHfwAH7", "reviewer": "Reviewer_yHuE", "summary": "The Gromov-Wasserstein optimal transport problem is a non-convex problem known to be hard, closely related to QAP. The authors consder its special case when the two involved metric spaces are Euclidean with small numbers of dimensions. We want to permute one set of points such that the sum of squares of differences of Euclidean distances between all point pairs is minimized. This simplification of this problem (compared to the general GW OT problem) is possible because the matrix of Euclidean pairwise distances for any set of points in $R^k$ has rank at most $k+2$. This allows the authors to formulate the problem as minimization of a simple concave quadratic function of a small ($kl+1$ where $k,l$ are the dimensions of the spaces) number of parameters over a convex polyhedral feasible set, which is a projection of the set of doubly-stochastic matrices. This problem is solved to global optimality by alternating two steps: (1) globally solve the concave minimization problem over an outer approximation of the feasible set, (2) improve this outer approximation by generating a cutting plane (a linear inequality valid for the above projection of the set of doubly-stochastic matrices). Step 1 can be done either by an off-the-shelf concave minimization algorithm (such as branch&bound) or by a reasoning over the set of extremal points of the feasible set (which is fast for low-dimensional spaces). Step 2 leads to solving the ordinary linear (Kantorovich) OT problem. The algorithm is proved to converge to a global optimum.\n\nThe proposed algorithm is first tested on synthetic data. Here, we match two randomly generated (either uniformly on a disc or normally distributed around the origin) point sets from $R^2$ or $R^3$. This is compared to other globally optimal methods solved by Gurobi and to the local search method [17] with multiple initializations. Second, the method is tested on a real application from biology. These experiments show that the method is almost always much faster than the other methods (the global methods are in fact usable only for small instances, due to the large number of parameters).", "questions": "Why were experiments on more real instance types not included, such as with shape/object matching? Is there a substantial obstacle or you just considered this as out of the scope of this paper?\n\nCan the method be extended to be resilient to outliers? One option for this would be to use 1-norm rather than 2-norm in pairwise distances - but this would violate the low-rank assumption and make the approach inapplicable.\n\nCan the trick be applied to a wider class of problems, such as a wider subclass of QAP?", "limitations": "The experiments are limited, not convincingly showing efficiency of the method on a wide enough class of real instances.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "A simple algorithm for a difficult problem, without many tuning parameters, which performs well in experiments.\nClearly written.", "weaknesses": "The most important issue is limited experimental testing. The main experiment is done on synthetic, randomly generated data. However, it is known that optimization algorithms often behave very differently on random data and on real data.\n\nMoreover, I do not quite understand the setup in the biological experiment, especially I get lost in the 1st paragraph in section 5.1. Please explain better. What is $n, l_x, l_y$ in this case?\n\nTo show more clearly the strengths and weaknesses of the approach, it should be tested on as many instance types as possible. Currently I am not convinced that the method would not last very (unacceptably) long for some instances from practice.\nOne option is to make more extensive synthetic experiments, with more complex data than just uniformly/normally distributed. E.g., one can take a non-random set of points in 2D or 3D (draw a shape or use a shape from a public shape database) and to synthetically generate the second set by permuting these points and adding noise to them (gaussian, uniform) or replace some of them with outliers.\n\nAnother obvious suggestion is shape/object matching on real data, where the point distances are measured by Euclidean metric (rather than geodesic). Possible inspirations for such experiments is [13, 17] (but, I believe, also other works).\n\nAs for novelty: Although the method just combines tools well-known in optimization, the main idea (decoupling the concave minimization part and the cutting plane part) is clever. Unfortunately, I do not know the relevant literature well enough to be sure that a similar approach has not been proposed before. I am surprised that no globally optimal algorithm for globally solving this problem has been proposed before by others - such algorithms have been apparently proposed only for general GW problem and QAP.\n\nThe initial formulation (1) of the GW problem seems to be wrong because, to my knowledge, the GW optimal-transport problem optimizes over doubly-stochastic matrices rather than permutation matrices (i.e., it allows \"soft matching\"). Therefore, the relaxation from (4) to (5) is unnecessary and in fact misleading. If (5) minimizes a concave function over a convex set, the two problems are indeed equivalent, but the explanation should go the other way.\n\nMinor issues / suggestions:\n\n- 77: \"the first two sums\" should be \"the 1st and 3rd sum\"\n\n- In (2), it would be logical to omit the 2nd and 3rd terms because they do not depend on Gamma (as noted above).\n\n- 92: You may wish to justify/cite why a Euclid distance matrix has rank at most $l_x+2$.\n\n- Proposition 1 should not be really a proposition, it is just a algebraic manipulation.\n\n- 133: The part defining the set $\\mathcal{F}$ is over-complicated and confusing. Why don't you write just \n$\\mathcal{F} = \\\\{ (2X\\Gamma Y, \\langle L,\\Gamma\\rangle) \\mid \\Gamma \\in P \\\\}$, which shows that $\\mathcal{F}$ is a linear map of $P$.\n\n- In the final version, Table 1 should report also the number $N$ of iterations (= added cutting planes) needed to achieve the prescribed accuracy. Currently, this is only in the supplement.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you sincerely for your comments and questions. See below, for our answers.\n\nWeaknesses:\n\nLimited experimental testing: \n\nThe main focus of this paper is on difficult problems that appear in, e.g., computational biology and where there are symmetries in the data and thus there is a large set of local optima. For many problems with image matching on simple figures (cats, dogs, etc.), both local search methods methods and the proposed global optimization method would be very quick. The local search would typically give the correct result and be quicker, but without any guarantees. Our method would be slower, but still fast, and come with guarantee. However, it makes sense to add some additional such examples, and illustrate that our method works also for this case. Thus we will include additional numerical simulations on known simple geometrical objects, such as MNIST. In brief, comparing two figures in MNIST takes in mean 0.7 seconds (standard deviation 0.39 seconds) running in mean 64 iterations (standard deviation 24 iterations) on problems of size 169 points in mean (standard deviation 30 points). In comparison [16] has a mean relative error of about 10% in its measurements of the GW discrepancy, which for some figures are more than the mean gap between the different classes\n\nSetup in the biological experiment: \n\nIn the biological experiment $n$ is the number of points, $\\ell_x$ is the dimensionality of the space which the points in $X$ is located in and $\\ell_y$ accordingly. In the biological experiment in section 5.3: $n=500$ and $\\ell_x=\\ell_y = 2$ (since we compare 2-d images). This shall be included in the text. Thank you for pointing this out.\n\n\n\nEuclidean metric: By using the Euclidean metric rending a high rank on the distance matrices with mixed positive and negative eigenvalues, such problems falls outside the class for which our paper considers. We are currently working on extending the results in this direction.\n\nGW problem: The original Gromov-Wasserstein description [13] indeed relax the problem to general marginals, and we shall add this in the introduction before Equation (1) to correct the terminology. Thank you for pointing this out. When the marginals are uniform over n points the result is still valid with permutations as the resulting matching, and closely related to the pre-model using couplings and correspondences in the Gromov-Hausdorff distance. One may consider discrete marginals $\\mu$ and $\\nu$ in this model. Then, the only difference is that $L = (1^T\\mu + 1^T\\nu) m_xm_y^T -4m_x\\nu ^T Y^T Y-4X^T X\\mu m_y^T$ and $c_0 =(\\langle C_x,C_x\\rangle+ \\langle C_y,C_y\\rangle -4\\nu^Tm_y \\mu^Tm_x)/2$.\n\nThe minor issues:\n\nline 77: Indeed, thank you for pointing this out!\nEqn (2): We decided to keep this as it is relevant for the discrepancy, even though not relevant for the optimization problem.\nLine 92: We can add this close to (3) where this is clear by a sum of matrices with rank 1, $\\ell_x$ and 1.\nProposition 1. Yes, but we added to increase readability and clarity.\nLine 133. Thank you for this suggestion!\nTable 1. Yes, we shall add the number of added cutting planes to the tables.\n\nQuestions:\n\nWe realize that there is a gap between the random data and the biologic example. We are planning to add a low dimensional geometrical data examples e.g. MNIST to show the performance and when symmetry becomes the key issue for local search methods. Such data is in scope for the method indeed.\n\nIn order to handle ouliers and 1-norm distance matrices we are violating the low rank assumption indeed. Currently such problems are not within the problem class considered in the paper, but we are looking to generalize the framework for handling such problems.\n\nAs the QAP problem is restricted to permutations, it would be a natural extension to consider a wider class of matrices. When the matrices contain both positive and negative eigenvalues, this becomes a much harder problem, which is currently out of scope of the method, but we are looking into handling such problems.\n\n"}, {"review_id": "WnNoYJTuGi", "reviewer": "Reviewer_xsEn", "summary": "This paper consider the Gromov-Wasserstien (GW) problem with quadratic cost, a (non-convex) quadratic optimization problem over the space of probability measures (in this work, restricted to uniform discrete measures supported on $n$ points) which takes the form (in discrete setting):\n$$\\min_{\\Gamma} \\sum_{1 \\leq i,j \\leq n} ||x_i - x_{i'}|^2 - |y_j - y_{j'}|^2|^2 \\Gamma_{ij} \\Gamma_{i'j'},$$\nwhere $\\Gamma$ is a doubly stochastic matrix of size $n \\times n$. \n\nThis problem is known to be hard to solve globally. This paper proposes an approach to do so that is tractable when the point clouds $X= (x_1,\\dots,x_n)$ and $Y = (y_1,\\dots,y_n)$ are in low dimension---say $d$ [edit: fix rendering]. The core idea is that the GW problem can be reparametrized as \n\n\\begin{equation}\\tag{1}\n\\min_{(W,w), \\Gamma} - |W|^2 - w + c_0,\n\\end{equation}\nwhere $\\Gamma$ is still a bistochastic matrix, $W,w$ must satisfy the relation $W = 2 X \\Gamma Y^T$ and $w = \\braket{L,\\Gamma}$, where $L$ and $c_0$ are constant term (wrt $w,W,\\Gamma$). Note that $W$ can be understood as the correlation between $X$ and $Y$ wrt the joint law $\\Gamma$, which is of size $d \\times d$, so losely speaking, the quadratic part of the GW problem does not depend on $\\Gamma$ (which is of size $n \\times n$) but only on the correlation it induces (which is a much smaller object when $d$ is small). \n\nFrom this reformulation, the idea of the paper to globally solve GW is the following: \n- Let $\\mathcal{F}$ denote the constraint polytope which links $W,w$ and $\\Gamma$. This space is too complex to be represented explicitly (by linear constraints). \n- Build a bounding box of $\\mathcal{F}$. This is doable because we know that the correlation matrix between $X$ and $Y$ has (explicit) bounds on its entries---for instance when $d=1$ this is simply the cost of the increasing (resp. decreasing) matching. \n- Build iteratively a sequence of \"lower approximation\" $(H_N)_N$ of $\\mathcal{F}$ (made of supporting hyperplanes of $\\mathcal{F}$), in sense that minimizing the functional in (1) over $H_N$ yields a lower bound $L_N$ for the GW problem. Note: this is a concave minimization problem, which is tractable in low dimension ($d \\times d + 1$ here). \n- Let $(W_N, w_N)$ be the solution on $H_N$. Typically $W_N,w_N \\not\\in \\mathcal{F}$ (otherwise, we have found a global minimizer of GW), but from this we can build a new constraint $H_{N+1}$ through solving a standard OT problem, providing a doubly stochastic matrix $\\Gamma_N$ which is by definition sub-optimal, hence giving an upper-bound $U_N$ for GW. \n- Eventually, the authors prove that $U_N - L_N \\to 0$ as $N \\to \\infty$, so the proposed approach yields a practical (in low-dimension) algorithm that converges to a **global** solution of the GW problem. ", "questions": "I think that the work may be improved by adding: \n- A discussion / numerical illustration on the convergence (rate) of the algorithm (Section 1. in Weaknesses), or an explanation of why this is not relevant. \n- A more precise comparison with concurrent works, in particular [17] (which is the closest to the current work as far as I can tell), from both a theoretical and numerical perspective ; where scalability but also quality of the result (which, I guess, may favor the proposed approach). If the comparison is not meaningful, please explain why. ", "limitations": "The authors clearly state that the work is dedicated to low-dimensional point cloud, which is its main limitation. \n\nAside from that, I do not see any potential societal impact specific to this work. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper considers the difficult (and important, in my opinion) question of globally solving the GW problem. \nThe proposed approach is based on existing ideas (parametrization of GW by low-dimensional matrices) but pushes them further to get an original and interesting approach. Having globally converging algorithms for GW, even if restricted to low dimensional setups with \"not-so-many points\", can be useful to assess the quality of other algorithms that may only converge locally (maybe some of them still converge globally \"most of the time\", etc.). \n\nThe paper is clear and does not sacrifice mathematical technicality. Proofs of Prop 2 and Thm 1 have been checked and no major flaw was identified (aside from small details). \n\nI also appreciate that the paper immediately acknowledges that its approach is limited to low-dimensional problems (which is not a major issue; using GW in low dimension is completely natural). ", "weaknesses": "# 1. On the convergence of the algorithm\n\nFrom my understanding, the proof of global convergence of the algorithm (Theorem 1) is purely asymptotic: the gap $\\epsilon_N = |U_N - L_N|$ is controled by $|\\theta_N - \\theta_{N+k}|,\\ (\\forall k)$ where $\\theta_N = (W_N, w_N)$ is compactly supported, which implies by contradiction that $0$ must be the single accumulation point of $\\epsilon_N$. \n\nWhile this makes sense (up to few technical considerations, see below), this can be considered as a weakness (at least from a theoretical viewpoint): the convergence of $\\epsilon_N \\to 0$ is controlled by \"how fast $\\theta_N$ accumulates\", and even in low dimension, without further investigation (i.e. thinking that $\\theta_N$ moves arbitrarily in a compact set), it may take a lot of time to reach low gap $\\epsilon$ (and this seems to be suggested by the experiments, where the criterion is set to $10^{-8}$ when $d=2$, but to $10^{-2}$ when $d=3$). \n\nI am not saying that this is what happen in practice, nor that this invalidates the contribution of this theoretical result, but a discussion on this may be welcome. \n\n\n# 2. On numerical experiments\n\nWhile the proposed experiments are conducted in a reasonable way, they remain somewhat limited in my opinion. In particular, \n- [related to the point above] I would have appreciated to have more illustrative experiments on the algorithm behavior, its convergence, etc. \n- As far as I can tell, the paper only compares with the work \"local search\" [16] (Peyré et al., 2016). Why not comparing with more modern works, as [17] (Scetbon et al., 2022), [D] by Sejourné et al. (2021) or [E] (Li et al., 2023)? Note that [17] seems to handle larger instances ($n=10^5$ points), but may fail to globally converge as far as I can tell. Showcasing the strengths of the current work vs [17] (probably the closest work in its spirit), even in illustrative scenarios, would be of interest. \n\n\n\n# Minor comments (rather suggestions than actual weaknesses)\n\n1. In Algorithm 1, the upper bound is updated as $U_{N+1} \\leftarrow \\min(U_N, \\text{OT cost} )$. From my understanding, this makes $U_N$ non-increasing, while $L_N$ is increasing (as a minimization problem with more and more constraints), so that $\\epsilon_N$ is decreasing. This is never mentioned as far as I can tell; but it seems to be used in the proof when saying that \"if $\\epsilon_N \\not\\to 0$, then it must be lower bounded\". Am I correct? Also, the proof says \"assume that $\\epsilon_N = \\dots$\" (which to me means \"assume that $U_{N+1} = \\text{OT cost}$ rather than $U_{N+1} = U_N$\"), but never discusses the possibility that $U_{N+1} = U_N$. In any case, this does not invalidate the proof as, from my understanding, the $\\min$ in the algorithm is not required to prove convergence (even if we cannot ensure that $\\epsilon_N$ is decreasing). But in any ways, I think that this is worth some clarification. \n\n2. The discussion related work and context can be improved. For instance, I do not think that [16] is a suited reference in the introduction with mentioning the GW problem for the first time. It would rather be credited to either Mémoli (2011), or Sturm (2012 - \"The space of spaces\"). What is new vs known (from [16], or for instance from [A, Sec 2.2.3] and related works) in sections 2 and 3 should also be highlighted. Similarly, using that the quadratic term in GW only depends on $X \\Gamma Y^T$ is also used (among other) in [B, C] (note : [C] is a preprint put on arxiv after neurIPS submission, this is a suggestion for the revised version, not a criticism). \n\n3. [typo] line 90 : \"has\" should be \"have\" I think?\n\n4. [typo] line 114, I think that \"maximum\" should be \"minimum\" ? \n\n5. I wonder how useful are the variables $Z_N, \\alpha_N$, given that $Z_N$ is simply $2 W_N$ and $\\alpha_N = 1$. I understand that this is a convenient way to write a \"general\" hyperplane equation $\\braket{Z,W} + \\alpha w \\leq \\beta$, but to me it turned out to hinder the reading a bit. \n\n6. [typo] line 168, \"that the that\". \n\n\n# References\n\n- [A] A contribution to Optimal Transport on incomparable spaces, T. Vayer, 2020\n- [B] On the existence of Monge maps for the Gromov-Wasserstien problem, Dumont et al., 2023.\n- [C] The Gromov-Wasserstein distance between spheres, Arya et al., 2023. \n- [D] The Unbalanced Gromov Wasserstein Distance: Conic Formulation and Relaxation, Sejourné et al., 2021\n- [E] A Convergent Single-Loop Algorithm for Relaxation of Gromov-Wasserstein in Graph Data, Li et al., 2023", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you sincerely for your comments and questions. See below for our answers.\n\nWeaknesses:\nOn the convergence of the algorithm:\n\nWe are working on quantifying the convergence rate. The number of iterations are bounded by  ($O((1/\\epsilon)^{\\ell_x\\ell_y+1})$), to cover the compact set, but the low dimensional QP problem is increasing in complexity for every iteration. The experiments in the supplementary material 1.4.1 suggest much faster rates, but this is work in progress.\n\nIf the problem is extremly symmetric, e.g. equidistant points on the surface of a sphere are being matched to equidistant points on another sphere, then the convergence would be very slow. We have not introduced symmetry breaking concepts into the scope of the method yet.\n\nOn numerical experiments:\n\n1. We are planning to add another experiment matching geometrical data e.g. MNIST objects, where we can highlight issues with symmetry using local methods.  In brief, comparing two figures in MNIST takes in mean 0.7 seconds (standard deviation 0.39 seconds) running in mean 64 iterations (standard deviation 24 iterations) on problems of size 169 points in mean (standard deviation 30 points). In comparison [16] has a mean relative error of about 10% in its measurements of the GW discrepancy which for some figures are more than the mean gap between the different classes.\n\n2. Indeed, the result in [17] and [D] increase the performance solving the linear OT, which can be incorporated in our algorithm as well. However the main purpose of this paper is to show a simple method to globally solve the GW problem with accuracy certificates. Since we have not optimized the subparts of our algorithm using the computational methods from, e.g., [17], we thing that is makes more sense to compare with [16]. Regarding scalability, note that our method fully decomposes the computations in terms of non-convex low dimensional quadratic problems and linear OT problem, thus one can apply any of the modern tricks on the OT problems in order to improve speed and problem sizes. This has not been a focus of this work and we have just used of the shelf solvers for the OT problem in order to simplify the description and implementation. \n\nMinor comments:\n\n1. Indeed, the upper bound $U_N$ is the least upper boud, which decreases, this the gap is monotonically decreasing. This can be stated more clearly.\n\n2. Indeed, this is an error that occured during redisposition of the manuscript. To our understanding, the Gromov-Wasserstein distance should be credited to Mémoli [13] or the 2007 proceedings article. Thank you for pointing out the newly published work!\n\n3-6. Indeed, thank you for finding the typos. \n\nQuestions:\n\n1. We have answered this under section weaknessess.\n\n2. The proposed method considers the topic of a global optimum of the GW distance. The main point of the paper is to show that one can find a global optimal solution within a reasonable time, and we do not claim to be faster than [16] or [17] on any GW-problems. The local search algorithms converge faster to a local optimum, but the advantage of our approach is that we are guaranteed to reach a global optimum. For many applications, this can be critical.\nAlso note that we have not focused on optimizing the computational complexity in each of the subparts of the algorithm. For example, we do not utilize Sinkhorn’s method for solving the optimal transport subproblems (Equation 9). It is true that [17] utilizes low rank structures similar to the ones we consider for improving the computational complexity of [16]. Similarly, we could use the methods from [Scetbon, Cuturi, Peyré - ICML, 2021] for improving the computational speed of the subproblems in our method. But since we have not optimized these computations, we think that it makes more sense to compare with the method in [16]. We will add this comment to the final version of the paper to motivate why we compare with [16] instead of [17].\n\n"}]

Author rebuttal:
Dear reviewers,

Thank you sincerely for the work you have put into reviewing the manuscript paper. Your reviews have been thorough, pointing at strengths and weaknesses in the paper and suggested clarifications that benefit the presentation and clarity of the paper. We will do our very best to incorporate the suggestions and remarks under the page limit constraint. In the rebuttal section in each review, we have answered questions and comments in order to make it more clear what answer goes to what question.

Again, thank you for your reviews and feedback.


Meta-review: All the reviews are positive. The paper introduces a novel algorithm for an important problem (GW in Euclidean spaces), leveraging the specific structure of the problem. The rebuttal have helped improved the presentation of the method towards a wider audience. The final version should include a better account for previous works on the same problem, for instance the PhD work of Vayer which already leveraged the concavity of GW between Euclidean spaces. I recommend acceptance.
Decision: Accept (poster)