domain: math stat

text:
Ernesto Araya

Hemant Tyagi
email: hemant.tyagi@inria.fr

Graph Matching via convex relaxation to the simplex

This paper addresses the Graph Matching problem, which consists of finding the best possible alignment between two input graphs, and has many applications in computer vision, network deanonymization and protein alignment. A common approach to tackle this problem is through convex relaxations of the NP-hard Quadratic Assignment Problem (QAP).

 in the noiseless setting. Our method is evaluated on both synthetic and real data, demonstrating superior statistical performance compared to existing convex relaxation methods with similar computational costs.

Introduction

In the context of the Graph Matching (GM) problem, also referred to as network alignment, we consider two input undirected graphs, G = ([n], E) and G ′ = ([n], E ′ ), which are assumed to share the same set of vertices, denoted as [n] = {1, • • • , n}. The objective is to find a bijective map, represented by a permutation x : [n] → [n], between the vertices of G and G ′ , which is intended to align their edges as effectively as possible. The alignment is considered optimal when the largest possible number of edges in G and G ′ satisfy the relation {i, j} ∈ E ⇐⇒ {x(i), x(j)} ∈ E ′ . This problem has numerous applications, such as computer vision [START_REF] Zhou | A survey on graph matching in computer vision[END_REF], network de-anonymization [START_REF] Narayanan | De-anonymizing social networks[END_REF], pattern recognition [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF][START_REF] Emmert-Streib | Fifty years of graph matching, network alignment and network comparison[END_REF], protein-protein interactions and computational biology [START_REF] Zaslavskiy | Global alignment of protein-protein interaction networks by graph matching methods[END_REF][START_REF] Singh | Global alignment of multiple protein interaction networks with application to functional orthology detection[END_REF], to name a few. For instance, it has been used in computer vision to track a moving object by identifying the correspondence between the different portions in two different frames.

Mathematically, the Graph Matching problem can be formulated as an instance of the NP-hard quadratic assignment problem, as follows max x∈Sn n i,j=1

A ij B x(i)x(j) ≡ max X∈Pn ⟨AX, XB⟩ F , (1.1)
where A, B ∈ R n×n are the (possibly weighted) adjacency matrices of the pairs of input graphs to be matched and P n (resp. S n ) is the set of n×n permutation matrices (resp. the set of permutation maps). In words, a solution of (1.1) will maximize the total edge alignment between A and B. For our purposes, it will be more useful to consider the following equivalent formulation, min

X∈Pn ∥AX -XB∥ 2 F , (1.2)
which has the advantage of possessing a convex objective function that we define by

E(X) := ∥AX -XB∥ 2 F , (1.3)
where the input graphs A and B will be clear from the context. It will be useful to consider the vector form of (1.2), and for that we write with a slight abuse of notation1 E(vec(X)) = vec(X) T H vec(X), (1.4) where H := (Id ⊗A -B ⊗ Id) 2 and vec(X) represents the vector form of X obtained by stacking its columns (notations are described in Section 1.2). Despite the convexity of E(•) the problem (1.2) is NP-Hard in the worst case [START_REF] Burkard | The quadratic assignment problem[END_REF], which can be mainly attributed to its combinatorial constraints, and even finding a constant factor approximation remains a hard problem [START_REF] Sahni | P-complete approximation problems[END_REF]. On the other hand, notice that when B = X * T AX * for some X * ∈ P n , (1.2) is equivalent to the well-known Graph Isomorphism Problem [START_REF] Babai | Group, graphs, algorithms: the graph isomorphism problem[END_REF].

A general algorithmic strategy In spite of its worst-case complexity, several efficient algorithms have been proven to solve the problem (1.5) within a specific class of instances. Most of the time, these instances are sampled by a statistical model, which generates matrices A, B according to a ground truth permutation X * ∈ P n (we will formally introduce statistical models in Section 2.3). The meta-strategy used by most graph matching algorithms consists of two steps that can be summarized as follows.

1. Similarity finding stage. We construct a data-dependent similarity matrix X(A, B) ∈ R n×n , where we expect that Xij (A, B) captures the likelihood of X * ij = 1, where X * is the ground truth permutation.

2. Rounding stage. The similarity matrix X is then prejected onto the set of permutations, which, when minimizing the Frobenius norm, is equivalent to solving a maximum linear assignment problem. Since analyzing linear assignment can be quite challenging, many algorithms opt for using a greedy rounding strategy. This approach is simpler to analyze, computationally more efficient, and often exhibits comparable error rates to linear assignment in experiments.

Convex relaxations of graph matching. One common approach for the similarity finding stage consists of defining X(A, B) as the solution of a relaxation of (1.2), where the set of permutation is relaxed to a continuous constraint set. Of special importance is the family of convex relaxations which can be expressed in general as min X∈K ∥AX -XB∥ 2 F .

(1.5)

where K ⊆ R n×n is a convex set such that P n ⊆ K. Some of the choices for the set K that have been investigated in the literature include the hyperplane defined by 1 T X1 = n [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF] and the Birkhoff polytope of doubly stochastic matrices [START_REF] Aflalo | On convex relaxation of graph isomorphism[END_REF]. The latter can be considered the gold standard, in terms of tightness of convex relaxations, but its analysis has proven to be challenging, and obtaining theoretical guarantees remains largely open.

Graph matching on the simplex and mirror descent. We consider the set of matrices

∆ n 2 := {X ∈ R n×n : 1 T X1 = 1, X ij ≥ 0, ∀i, j ∈ [n]},
which can be identified with the usual unit simplex in R n 2 , in the sense that X ∈ ∆ n 2 if and only if vec(X) belongs to the unit simplex (we will refer to ∆ n 2 as the unit simplex in the sequel). By choosing K = ∆ n 2 in (1.5), we arrive to the following convex relaxation min

X∈∆ n 2 ∥AX -XB∥ 2 F . (1.6)
This formulation is referred to as graph matching on the simplex. Despite (1.6) representing a natural problem and being an advancement over other convex relaxations (e.g., it is tighter than the 1 T X1 = n constraint in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]), it has been largely overlooked in the literature. We argue that (1.6) presents at least two advantages over existing work.

• Regularization through positivity. In the works [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF][START_REF] Zhou Fan | Spectral graph matching and regularized quadratic relaxations II: Erdős-Rényi graphs and universality[END_REF], the hyperplane constraint K = {X ∈ R n×n : 1 T X1 = n} was considered and an additional regularizing term of the form η∥X∥ 2 F was added to the objective E(•). By directly incorporating the positivity constraints, we achieve a similar regularization effect without the need for a cumbersome explicit regularizer. This approach helps us avoid the task of parameter selection for such a parameter in practice. We show in experiments that it outperforms [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF][START_REF] Zhou Fan | Spectral graph matching and regularized quadratic relaxations II: Erdős-Rényi graphs and universality[END_REF].

• Efficient algorithms. we can solve (1.6) using iterative first-order methods such as projected gradient descent(PGD) and entropic mirror descent(EMD). Although the existence of efficient projection algorithms onto ∆ n 2 can make PGD within the simplex more expedited compared to PGD within the Birkhoff polytope, employing the EMD algorithm allows us to completely bypass this projection step, as it possesses an explicit update formula in the form of multiplicative weights updates [START_REF] Bubeck | Convex optimization: Algorithms and complexity[END_REF][START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF].

Contributions. The primary contributions of this paper are outlined below. To our knowledge, these findings mark the first results for both convex graph matching on the simplex and the mirror descent (MD) algorithm for graph matching.

• We propose an efficient iterative algorithm based on the mirror descent to solve the simplex relaxation of graph matching, which improves in O(log n) in running time with respect to the more common projected gradient descent method.

• We introduce a novel sufficient condition for achieving exact recovery, which is notably less stringent than the commonly used diagonal dominance property. We employ this condition to establish exact one-step recovery (i.e., in one iteration) of the ground truth permutation via EMD, in the noiseless case of the Correlated Gaussian Wigner model. This result holds almost surely, and for any positive step size in EMD. Furthermore, we use this condition to improve a pivotal lemma from [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF], yielding stronger results with a simpler proof.

• In the noiseless case of the Correlated Gaussian Wigner model, we prove that (1.6) has a unique solution with high probability, with the solution corresponding to a scaled version of the ground truth permutation. This, along with [START_REF] Beck | First Order Methods in Optimization[END_REF]Theorems 9.16 and 9.18], ensures that EMD and PGD converge to a unique point (under certain conditions on step size), which in the noiseless setting leads to exact recovery of the ground truth permutation.

• Through extensive experiments, we demonstrate that our proposed algorithm performs well on synthetic data, for several correlated random graph models (described in Section 2.3). In addition, we show the applicability of our algorithm to real-world datasets where it is seen to outperform other methods either in speed or accuracy.

Related work

Convex relaxations. Although many algorithms based on solving a problem of the form (1.5) have been proposed [START_REF] Ganassali | Spectral alignment of correlated gaussian matrices[END_REF][START_REF] Feizi | Spectral alignment of graphs[END_REF][START_REF] Lyzinski | Graph matching: relax at your own risk[END_REF][START_REF] Zaslavskiy | A path following algorithm for the graph matching problem[END_REF], some with very good experimental performance, there is still a lack of theoretical understanding regarding their statistical guarantees. In the case where K is the set of doubly stochastic matrices, some early theoretical work, which applies under restrictive conditions on the inputs graphs and mostly in the deterministic setting, can be found in [START_REF] Aflalo | On convex relaxation of graph isomorphism[END_REF]. A spectral estimator dubbed GRAMPA, which is the solution of a regularized version of (1.5), taking K to be the hyperplane 1 T X1 = n, was studied in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]. They prove that, after greedy rounding, their estimator perfectly recovers the ground truth permutation when the graphs are generated with a correlated Gaussian Wigner model, in the regime where the correlation is 1 -O 1 log 2 n . There is still a big gap between these guarantees and the sharp information theoretic threshold for recovery, as discusses below.

Information-theoretic limits of graph matching. The necessary and sufficient conditions for correctly estimating the matching between two graphs when they are generated at random, from either the correlated Gaussian Wigner or the correlated Erdős-Rényi models (see Section 2.3), have been investigated by several authors in [START_REF] Cullina | Exact alignment recovery for correlated Erdös-Rényi graphs[END_REF][START_REF] Hall | Partial recovery in the graph alignment problem[END_REF][START_REF] Wu | Settling the sharp reconstruction thresholds of random graph matching[END_REF][START_REF] Ganassali | Impossibility of partial recovery in the graph alignment problem[END_REF][START_REF] Ganassali | Sharp threshold for alignment of graph databases with gaussian weights[END_REF][START_REF] Mao | Testing network correlation efficiently via counting trees[END_REF][START_REF] Ganassali | Statistical limits of correlation detection in trees[END_REF]. For the correlated Gaussian Wigner model, it has been shown in [START_REF] Wu | Settling the sharp reconstruction thresholds of random graph matching[END_REF]Thm.1] that the ground truth permutation x * can be exactly recovered w.h.p. only when

σ 2 ≤ 1 -(4+ϵ) log n n . When σ 2 ≥ 1 -(4-ϵ) log n n
no algorithm can even partially recover x * . One natural question pertains to the possibility of closing the gap between the performance bound for GRAMPA as outlined in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF] and the information-theoretic threshold by utilizing convex relaxations.

MD algorithm and Multiplicative Weights Updates. The MD algorithm, introduced in [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF], has been extensively studied in convex optimization [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Bubeck | Convex optimization: Algorithms and complexity[END_REF], learning theory, and online learning [START_REF] Hazan | Introduction to online convex optimization[END_REF]. On the other hand, the EMD algorithm, also known as Exponentiated Gradient Descent (EGD) [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] or exponential Multiplicative Weights Updates (MWU) [START_REF] Arora | The multiplicative weights update method: a meta-algorithm and applications[END_REF], has found diverse applications spanning machine learning [START_REF] Bernstein | Learning compositional functions via multiplicative weight updates[END_REF], optimization [START_REF] Arora | A combinatorial, primal-dual approach to semidefinite programs[END_REF][START_REF] Arora | Fast algorithms for approximate semidefinite programming using the multiplicative weights update method[END_REF][START_REF] Anderson | An efficient algorithm for unweighted spectral graph sparsification[END_REF], computer science [START_REF] Arora | The multiplicative weights update method: a meta-algorithm and applications[END_REF][START_REF] Allen-Zhu | Spectral sparsification and regret minimization beyond matrix multiplicative updates[END_REF], and game theory [START_REF] Daskalakis | Near-optimal no-regret algorithms for zero-sum games[END_REF], among others. Although prior work, such as [START_REF] Ding | Nonnegative matrix factorization for combinatorial optimization: Spectral clustering, graph matching, and clique finding[END_REF] and [START_REF] Jiang | Graph matching via multiplicative update algorithm[END_REF], has explored MWU strategies for graph matching, there are notable distinctions from our research. Firstly, they consider different objective functions and constraint sets, yielding fundamentally different update strategies. Secondly, their approach lacks a principled assessment of robustness against noise due to the absence of statistical generative models for inputs A and B. It is worth noting that the same exponential MWU technique we explore here was previously considered in the context of the nonconvex formulation of Graph Isomorphism in [START_REF] Pelillo | Replicator equations, maximal cliques, and graph isomorphism[END_REF]. However, it was primarily proposed as a heuristic inspired by the replicator equations from theoretical biology, and no connection with the MD algorithm was established in that work.

Other algorithms for seedless graph matching. Although less related to our work, it is worth mentioning that there exist several approaches for seedless graph matching not directly based on convex relaxations. When K is the set of orthogonal matrices in (1.5), its solution can be expressed in closed form in terms of the spectrum of the graphs, as was proven in the celebrated work [START_REF] Umeyama | An eigendecomposition approach to weighted graph matching problems[END_REF]. In [START_REF] Mao | Exact matching of random graphs with constant correlation[END_REF], the authors prove that a two-step algorithm can attain exact recovery in the context of the sparse Erdős-Rényi graphs, even in the constant correlation regime. For sparse Erdős-Rényi graphs, other approaches based on counting combinatorial structures include [START_REF] Mao | Random graph matching at otter's threshold via counting chandeliers[END_REF][START_REF] Ganassali | From tree matching to sparse graph alignment[END_REF]. Exploring the possibility of extending those approaches to the dense case remains an open question.

Notation.

We denote P n (resp. S n ) the set of permutation matrices (resp. permutation maps) of size n. We define the vectorization operator vec : R n×n → R n 2 ×1 , which, for a matrix M , vec (M ) gives the corresponding vector formed by stacking the columns of M . For any given matrix M , we define the norm ∥M ∥ 1,1 := i,j |M ij | and the usual Frobenius norm ∥M ∥

F := ( i,j |M ij | 2 ) 1 2 = ⟨M, M ⟩ F , where for matrices M, M ′ , ⟨M, M ′ ⟩ F := ⟨vec(M ), vec(M ′ )⟩ (here ⟨•, •⟩ is the usual inner product in R n 2 ×1
). Similarly, for M, M ′ ∈ R n×n , we define M ⊙ M ′ to be their Hadamard (entrywise) product and M ⊗ M ′ their usual Kronecker product. For a matrix M we define its Hadamard (entrywise) exponential exp ⊙ (M ) to be matrix with (i, j)-th entry exp(M ij ). We use 1 (resp. J) to denote the vector (resp. matrix) of all ones in R n (resp. R n×n ), where the size n will be clear from context. We define ∆ n := {x ∈ R n : 1 T x = 1} and we identify ∆ n 2 with the set of matrices M ∈ R n×n such that vec(M ) ∈ ∆ n 2 . We denote Id the identity matrix, where the dimension will be clear from context. We will denote e 1 , . . . , e n the elements of the canonical basis of R n .

Algorithms 2.1 A greedy rounding procedure and a sufficiency lemma

We consider the following greedy rounding algorithm, named Greedy Maximum Weighted Matching (GMWM) algorithm in [START_REF] Lubars | Correcting the output of approximate graph matching algorithms[END_REF][START_REF] Yu | Graph matching with partially-correct seeds[END_REF], which will serve us to accomplish Step 2 in the meta-strategy described earlier. The specific version presented below, as Algorithm 1, has been taken from [START_REF] Araya | Seeded graph matching for the correlated Wigner model via the projected power method[END_REF]. To simplify our analysis, we will consider the notion of algorithmic equivariance, which, in the Algorithm 1 GMWM (Greedy maximum weight matching) Require: A cost matrix C ∈ R n×n . Ensure: A permutation matrix X ∈ P n .

1: Select (i 1 , j 1 ) such that C i 1 ,j 1 is the largest entry in C (break ties arbitrarily). Define C (1) ∈

R n×n : C (1) ij = C ij 1 i̸ =i 1 ,j̸ =j 1 -∞ • 1 i=i 1 or j=j 1 . 2: for k = 2 to N do 3: Select (i k , j k ) such that C (k-1)
i k ,j k is the largest entry in C (k-1) (break ties arbitrarily).

4:

Define C (k) ∈ R n×n : C (k) ij = C (k-1) ij 1 i̸ =i k ,j̸ =j k -∞ • 1 i=i k or j=j k . 5: end for 6: Define X ∈ {0, 1} n×n : X ij = N k=1 1 i=i k ,j=j k . 7: return X
context of graph matching, means that if an algorithm outputs an estimator XP for inputs A and B, then the same algorithm outputs XP Π T when the inputs are A and ΠBΠ T instead, for any Π ∈ P n . The following lemma tells us that any algorithm that, first, obtains X by solving (1.5) for a permutation-invariant set K and, second, rounds X using Algorithm 1, is equivariant. We include its proof in Appendix B for completeness. Lemma 1. Let K be a convex subset of R n×n such that KΠ T = K, for any Π ∈ P n , and for A, B ∈ R n×n define S(A, B) to be the set of solutions of (1.5). Assume that the set defined by {Y ∈ R n×n : Y = GMWM( X) for X ∈ S(A, B)}, is a singleton. Then the algorithm that outputs XP = GMWM( X) for any X ∈ S(A, B) is equivariant.

The previous lemma implies that we can consider, without loss of generality, the ground truth to be X * = Id (see [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Section 1.2] for additional discussion). This motivates the following deterministic lemma. Lemma 2. Let C ∈ R n×n be any matrix. Then the following statements hold.

(i) If C satisfies C ii ∨ C jj > C ij ∨ C ji , ∀i ̸ = j ∈ [n], (2.1)
then GMWM(C) = Id.

(ii) If C is symmetric and

C ii + C jj > 2C ij , ∀i ̸ = j ∈ [n], (2.2)
then GMWM(C) = Id .

(iii) If C is symmetric and positive definite, we have GMWM(C) = Id.

In (i) and (ii), if for each i ∈ [n] the set arg max j C ij is a singleton, then the sufficient conditions (2.1) and (2.2) are also necessary.

Proof. Notice that in the first step of GMWM we select the indices (i 1 , j 1 ) corresponding to the largest entry of C. Those indices satisfy i 1 = j 1 , otherwise this would contradict (2.1). Then it suffices to notice that in the subsequent steps the matrix C (k) still satisfy (2.1) which implies that i k = j k for all k ∈ [n]. This proves (i). When C is symmetric, the property (2.2) implies (2.1). This proves (ii). On the other hand, if C is positive definite, then for all i ̸ = j ∈ [n] we have (e i -e j ) T C(e i -e j ) > 0, where e k is the k-the element of the canonical basis of R n . From this, we deduce that C ii + C jj > 2C ij , which directly implies (2.1), proving (iii). To prove the converse, in the case (i), notice that if GMWM(C) = Id and for each i ∈ [n] the set arg max j C ij is a singleton, then there are no ties in lines 1 and 3 of Algorithm 1, and it is clear that the maximum for each matrix C (k) is attained by some element in its diagonal. This implies (2.1). The proof for the converse of (ii) is analogous.

In terms of the general algorithmic strategy sketched in Section 1, the main use of Lemma 2 will be the following. If we can prove that X, a solution of the convex relaxation (1.5), satisfies (2.1) with certain probability, then Lemma 2 implies that X P := GMWM( X) is equal to X * = Id (that is, we obtain exact recovery), with the same probability. An interesting point is to understand if the properties (2.1) and (2.2) are preserved under typical matrix operations such as additions, scaling, etc. In this direction, we have the following result, whose proof is direct. Lemma 3. Let C ∈ R n×n be a matrix satisfying (2.1), then for any t > 0, α ∈ R, the matrices tC and C + αJ also satisfy (2.1). Let C ′ , C ′′ ∈ R n×n be any two symmetric matrices satisfying (2.2), then tC ′ ,C + αJ and C ′ + C ′′ satisfy (2.2).

In words, the previous lemma states that property (2.1) is closed under scaling by a positive factor and a global additive translation, while (2.2) is also closed under addition.

Mirror descent algorithm for graph matching

Our proposed method consists of two main steps: first, solving (1.6), which addresses the graph matching problem on the unit simplex and yields X; second, rounding this solution using Algorithm 1 to generate X P an estimate of the ground truth X * . To efficiently solve (1.6), we propose the utilization of the mirror descent (MD) method. This first-order iterative optimization algorithm generalizes projected gradient descent by allowing for non-Euclidean2 distances [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF]. To avoid unnecessary generality, we will focus on two primary examples within the MD family on the unit simplex: Projected Gradient descent (PGD) and Entropic Mirror Descent (EMD). As mentioned ealier, EMD has also received the name of Exponentiated Gradient Descent (EGD) [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] and can be seen as an example of the Multiplicative Weights Update (MWU) [START_REF] Arora | The multiplicative weights update method: a meta-algorithm and applications[END_REF] strategy.

EMD and PGD updates. Both methods receive an initial point X (0) as input but differ in their approach to performing updates. More specifically, the update process for these methods, in the case of problem (1.6), is as follows,

X (k+1) = arg min X∈∆ n 2 n i,j=1 γ k ∇E(X (k) ) ij -X (k) ij X ij + X 2 ij , ( PGD update)
X (k+1) = arg min X∈∆ n 2 n i,j=1 γ k ∇E(X (k) ) ij -1 -log X (k) ij X ij + X ij log X ij , ( EMD update)
where (γ k ) k≥0 is a sequence of positive learning rates (also called the step-sizes). The derivation of these updates is well-known (refer to [START_REF] Beck | First Order Methods in Optimization[END_REF]Section 9]) and, without entering into details, stems from a broader form of MD, which permits the selection of different (Bregman) distances to project onto the set of constraints. In (PGD update), the Euclidean distance is used as the Bregman distance, and it is easy to see that (PGD update) is equivalent to

X (k+1) = proj ∆ n 2 X (k) -γ k ∇E(X (k) ) , (2.3)
where proj ∆ n 2 corresponds to the Euclidean (Frobenius norm) projection onto ∆ n 2 and ∇E(X) = A 2 X + XB 2 -2AXB. In the case of (EMD update), the Bregman distance employed is the Kullback-Liebler divergence, which gives rise to the term 'entropic mirror descent'. Perhaps more importantly, from (EMD update) one can easily derive (as shown in [START_REF] Beck | First Order Methods in Optimization[END_REF]Example 3.17]) the following simplified form (which corresponds to a MWU rule)

X (k+1) ij = X (k) ij exp -γ k ∇E(X (k) ) ij + c , (2.4)
where c is a constant which ensures that X (k+1) ∈ ∆ n 2 . Our proposed algorithm, written in matrix language, is summarized in Algorithm 2.

Initialization. We chose to always start with X (0) = J/n 2 , the barycenter of the unit simplex, which can be seen as an agnostic initialization, given that the Frobenius distance from X (0) to any permutation is the same. Additionally, this choice satisfies the known boundedness hypothesis for Kullback-Leibler divergence, with respect to the initial point, that guarantees O 1 √ N convergence to the optimal solution in (1.6) (see [START_REF] Beck | First Order Methods in Optimization[END_REF]Theorem 9.16] and [START_REF] Beck | First Order Methods in Optimization[END_REF]Example 9.17]).

Algorithm 2 Entropic mirror descent for Graph Matching (EMDGM)

Require: A, B ∈ R n×n , N ∈ N and sequence {γ k } N k=0 ∈ R + . Ensure: A permutation X P ∈ P n . 1: Define the initial point X (0) = J/n 2 . 2: Initialize X best ← X (0) 3: for k = 0 to N do 4: Compute N k = ∥X (k) ⊙ exp ⊙ -γ k (A 2 X (k) + X (k) B 2 -2AX (k) B) ∥ 1,1 5: Update X (k+1) = X (k) ⊙ exp ⊙ -γ k (A 2 X (k) + X (k) B 2 -2AX (k) B) /N k 6: if E(X (k+1) ) < E(X best ) then 7:
Reassign X best ← X (k+1)

8:

end if 9: end for 10: return X P := GMWM(X best ) Remark 1 (Complexity of Algorithm 2). Each iteration within the for loop of Algorithm 2 has two main steps. First, the gradient computation, which has the same time complexity as matrix multiplication, O(n ω ) where ω ≤ 2.373 (as established in [START_REF] François | Algebraic complexity theory and matrix multiplication[END_REF]). Second, the Hadamard product of X (k) and the entrywise exponentiated gradient, which has complexity n 2 . Similarly, it is clear (see discussion in [START_REF] Yu | Graph matching with partially-correct seeds[END_REF]) that Algorithm 1 has complexity O(n 2 ). Hence, overall, Algorithm 2 has a time complexity of O(N n ω ).

Generative models for correlated random graphs

We now recall the two most popular models for the statistical study of the graph matching problem. These models generate a pair of correlated adjacency matrices A, B, based on a ground truth permutation X * ∈ P n . Noteworthy fact for these models is that the solution of the QAP formulation (1.1) coincides with the maximum likelihood estimator of X * . This is not true for other models, as shown in [START_REF] Wang | Random graph matching in geometric models: the case of complete graphs[END_REF], in the case of random graph models with an underlying geometric structure.

Correlated Gaussian Wigner(CGW) model W(n, σ, X * ). The Correlated Gaussian Wigner (CGW) model, introduced in [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF], defines a distribution on pairs of complete weighted graphs, such that each of them is (marginally) distributed as a random symmetric Gaussian matrix and their correlation is described by a single parameter σ. More specifically, we will say that A, B ∈ R n×n follow the CGW model, and write A, B ∼ W(n, σ, X * ) if B = X * T AX * + σZ where σ > 0 and A, Z are i.i.d GOE(n) and X * ∈ P n . We recall that a n × n matrix A follows the GOE(n) law if its entries are distributed as

A ij ∼ N (0, 1 n ) if i < j, N (0, 2 n ) if i = j, and A ij = A ji for all i, j ∈ [n].
Correlated Erdős-Rényi (CER) model G(n, σ, p, X * ). For σ, p ∈ [0, 1], we describe the correlated Erdős-Rényi model with latent permutation X * ∈ P n by the following two-step sampling process.

1.

A is sampled according to the classic single graph Erdős-Rényi model G(n, p), i.e. for all i < j, A ij are independent Bernoulli's random variables with parameter p, A ji = A ij and A ii = 0.

2. Conditionally on A, the entries of B are i.i.d according to the following law

(X * T BX * ) ij ∼ Bern 1 -σ 2 (1 -p) if A ij = 1, Bern σ 2 p if A ij = 0. (2.5)
Another equivalent description of this model in the literature involve first sampling an Erdős-Rényi "mother" graph, and then obtaining A, B as independent subsamples with certain density parameter (related to σ). We refer to [START_REF] Pedarsani | On the privacy of anonymized networks[END_REF] for details on this alternative formulation.

Theoretical results

This section is dedicated to presenting our theoretical results. In Section 3.1, we prove that (1.6) has a unique solution with high probability, which, in the noiseless setting, corresponds to the scaled ground truth permutation. This, along with [11, Theorems 9.16 and 9.18], ensures that EMD and PGD converge to a unique point (under certain conditions on {γ k } N -1 k=0 and N ), which in the noiseless setting leads to exact recovery of the ground truth permutation. In Section 3.2, we examine the EMD dynamic (2.4) in the noiseless setting (σ = 0). We establish, as demonstrated in Theorem 2 below, that in this scenario, a single iteration of Algorithm 2 is sufficient to almost surely recover the ground truth permutation for any positive step-size γ 0 . In Section 3.3, we prove a stronger version, as shown in Theorem 3, of a critical lemma originally found in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF] (specifically, [26, Lemma 2.3]) for the GRAMPA algorithm. This lemma forms the cornerstone of GRAMPA's theoretical guarantees for the general noisy setting.

Uniqueness of the solution of simplex Graph Matching

In this section, we will prove that under the CGW model, the solution of (1.6) is unique with high probability. As usual, we assume without loss of generality that X * = Id. Note that if σ > 0, then the objective ∥AX -XB∥ 2 F is strongly convex a.s, so the solution of (1.6) is unique a.s. Hence, we focus on the noiseless case where (1.6) reduces to min

X∈∆ n 2 ∥AX -XA∥ 2 F . (3.1)
Notice that under the CGW model, A ∼ GOE(n), and it is well-known that A has distinct eigenvalues almost surely [4, Theorem 2.5.2]. In this case, we know that X = 1 n Id ∈ ∆ n 2 is a solution of (3.1), which implies that any solution has to satisfy ∥AX -XA∥ F = 0 or, equivalently, AX -XA = 0. In vector notation, this can be written as vec (AX -XA) = H vec (X) = 0, where we recall that H = (Id ⊗A -A ⊗ Id). It is easy to see that the null space of H is given by span {v i ⊗ v i } n i=1 , where v 1 , . . . , v n are the eigenvectors of A. Therefore, all the solutions of (3.1) satisfy two conditions: first, they are of the form X = n i=1 µ i v i v T i for some reals µ 1 , . . . , µ n (i.e. X ∈ span{v 1 v T 1 , . . . , v n v T n }), and second, they belong to ∆ n 2 .

Theorem 1. Let A ∼ GOE(n), then with probability larger than 1 -n 2 n-1 , the problem (3.1) has a unique solution.

Before proving Theorem 1, we first recall some concepts that will be useful in the proof. The following definitions are taken from [START_REF] Horn | Matrix Analysis[END_REF]Section 6.2]. A matrix M ∈ R n×n is reducible if there exist

P ∈ P n , r ∈ [n -1] and matrices M 11 ∈ R r×r , M 12 ∈ R (n-r)×r , M 22 ∈ R (n-r)×(n-r) such that P M P T = M 11 M 12 0 (n-r),r M 22 ,
where 0 (n-r),r is the matrix with all zero entries in R (n-r)×r . On the other hand, a square matrix is irreducible if it is not reducible. We will denote I the set of irreducible matrices and R the set of reducible matrices (their dimension will be clear from context). It is easy to see that if M is symmetric, then

M ∈ R ⇐⇒ ∃m < n, (B kk ) m k=1 with B kk ∈ R n k ×n k and B kk ∈ I ∃P ∈ P n s.t P M P T = blockdiag(B 11 , . . . , B mm ). (3.2)
This block diagonal representation (3.2) is referred to as the Frobenius (or irreducible) normal form of M . We will use repeatedly the fact that if M ∈ I and has nonnegative entries, then the "strong" form of the Perron-Frobenius theorem (see [START_REF] Horn | Matrix Analysis[END_REF]Theorem 8.4.4]) applies. As a consequence, the largest modulus eigenvalue of M ∈ I is (algebraically) simple and its associated eigenvector can be chosen with all its entries strictly larger than zero.

Proof of Theorem 1. Let S be the (random) set of solutions3 of (3.1). By the remarks at the beginning of this section, we have

S = span {v i v T i } n i=1 ∩ ∆ n 2 ,
where v 1 , . . . , v n are the eigenvectors of A. Note that S is nonempty, because clearly 1 n Id ∈ S. Our objective is to control the probability of the event

E := S ̸ = 1 n Id .
Consider R 1 to be the class of reducible matrices with a Frobenius normal form having at least one block of size one, and similarly, R 2 will denote the class of reducible matrices with a Frobenius normal form with all of its blocks of size two or more (so clearly

{R 1 , R 2 } form a partition of R).
With this definition, it is easy to see that

E = {S ∩ I ̸ = ∅} =:E 1 ∪ S ∩ R 1 ̸ = 1 n Id =:E 2 ∪ {S ∩ R 2 ̸ = ∅} =:E 3 .
Here, the events E 1 , E 2 , and E 3 are clearly disjoint. We will bound their probabilities separately.

Bounding P(E 1 ). Notice that if there exists a matrix X in S ∩ I, then for that X, the Perron-Frobenius theorem applies, which justifies the implication

X ∈ S ∩ I =⇒ X ∈ span {v i v T i } n i=1 , and ∃i ′ ∈ [n] s.t. v i ′ (k) > 0 ∀k ∈ [n], or v i ′ (k) < 0 ∀k ∈ [n],
where v i ′ (k) denotes the k-th coordinate of v i ′ . From this we deduce that

E 1 ⊆ {∃i ′ s.t. v i ′ (k) > 0, ∀k ∈ [n]} ∪ {∃i ′ s.t. v i ′ (k) < 0, ∀k ∈ [n]}.
We will need the following lemma, whose proof can be found in Appendix C. Lemma 4. Let v 1 , v 2 , . . . , v n be a set of (not necessarily independent) vectors uniformly distributed on the unit sphere S n-1 . Then it holds

P(∃i ′ ∈ [n] s.t. v i ′ (k) > 0, ∀k ∈ [n]) ≤ n 2 n .
Given that v 1 , . . . , v n are uniformly distributed in the unit sphere S n-1 , we use Lemma 4 and the union bound, to conclude that

P(E 1 ) ≤ n 2 n-1 .
(3.3)

Bounding P(E 2 ). Fix any X ∈ S ∩ R 1 .
We have that there exists at least one block in the Frobenius normal form of X (c.f. (3.2)) with size one. For simplicity, we can assume w.l.o.g. that its first block is of size one, i.e., B 11 ∈ R. Given that entries of X are nonnegative, it is clear that B 11 ≥ 0. Then, the following chain of implications holds (denoting e 1 the first element of the canonical basis of R n )

P XP T e 1 = B 11 e 1 =⇒ XP T e 1 = B 11 P T e 1 =⇒ n i=1 µ i v i ⟨v i , P T e 1 ⟩ = B 11 P T e 1 =⇒ V Dµ = B 11 P T e 1 , (3.4)
where V has as columns v 1 , . . . , v n , D is a diagonal matrix with entries D ii = ⟨v i , P T e 1 ⟩ and µ = (µ 1 , . . . , µ n ) T . We have the following simple result.

Claim 1. D is almost surely invertible. Proof. Clearly D is invertible iff D ii ̸ = 0 for all i ∈ [n]. Note that, there exists k ′ ∈ [n] such that P T e 1 = e k ′ . Hence, D ii = ⟨v i , e k ′ ⟩ for some k ′ ∈ [n]. We claim that P ∀k, i ∈ [n], ⟨v i , e k ⟩ ̸ = 0 = 1. (3.5)
Indeed, for any given pair i, k ∈ [n], we have that P(⟨v i , e k ⟩ = 0) = 0 (because v i is uniformly distributed on S n-1 ). Then (3.5) follows from the union bound. From (3.5), we deduce that D ii ̸ = 0 a.s.

Claim 1 implies that V D is almost surely invertible, which in turn implies that there exists a unique µ that satisfies V Dµ = B 11 P T e 1 . On the other hand, it is clear that µ = B 11 1 is a solution of (3.4), so it has to be the unique solution. But, if µ = B 11 1, then X = B 11 Id. However, the fact that X ∈ S implies that B 11 = 1 n . In summary, we have demonstrated that for any X ∈ S ∩ R 1 , it holds that X = 1 n Id almost surely. In other words, with probability one, there are no elements in S ∩ R 1 other than 1 n Id. Equivalently, P(E 2 ) = 0.

Bounding P(E 3 ). Fix any X ∈ S ∩ R 2 (assuming it exists). Then for some P ∈ P n , we have

P XP T = blockdiag(B 11 , . . . , B mm ) =:B , (3.6)
with all the blocks B kk ∈ R n k ×n k having size n k ≥ 2, and B kk ∈ I. We note that AX = XA ⇐⇒ P AP T B = BP AP T .

We will use the following lemma.

Lemma 5. For A ∼ GOE(n) and B as defined in (3.6), it holds

P ∀Q ∈ P n , QAQ T B ̸ = BQAQ T = 1.
Proof. Fix Q ∈ P n . Consider the following block decomposition of à := QAQ T .

à =       Ã11 Ã12 • • • Ã1m Ã21 Ã22 . . . . . . . . . . . . Ãm1 • • • Ãmm ,      
where Ãkl belongs to R n k ×n l (we use bold notation here to distinguish Ãkl from the (k, l) entry of Ã). Observe that Ãkl = ÃT

It can be seen that the following facts hold.

(i) Ã dist. = A ∼ GOE(n).
In particular, Ã11 ∼ GOE(n 1 ), which implies that Ã11 has n 1 different eigenvalues almost surely (c.f. [4, Thm. 2.5.2]). Hence, its eigenvectors, denoted ṽ(1) 1 , . . . , ṽ

n 1 (1)
are uniquely defined up to sign change.

(ii are uniquely determined up to sign change.

) Ã12 Ã21 = Ã12 ÃT 12 ∼ Wish(n 1 , n 2 ). Here M ∼ Wish(n 1 , n 2 ) means that M ∈ R n 1 ×n 1 can be written M = M M T , where M ∈ R n 1
(iii) Ã12 Ã21 is independent from Ã11 .

Given that B 11 ∈ I, we deduce, by Perron-Frobenius theorem, that its largest eigenvalue has multiplicity one. On the other hand, by (3.7) we have that

B 11 ∈ span {ṽ (1) i ṽ(1)T i } n 1 i=1 ∩ span {ṽ (1,2) i ṽ(1,2)T i } n 1 i=1 .
In other words, {ṽ

1 , . . . , ṽ(1) n 1 } and {ṽ

(1,2) 1

, . . . ṽ(1,2) n 1 } are both basis of eigenvectors for B 11 . Assume w.l.o.g that the largest eigenvalue of B 11 (of multiplicity one) is associated with ṽ(1) 1 , in the first basis, and ṽ(1,2)

1
, in the second. This implies that either ṽ(1)

1 = ṽ(1,2) 1 or ṽ(1) 1 = -ṽ (1,2) 1
. From facts (i) and (ii), ṽ

1 and ṽ(1,2)

1 are uniformly distributed in S n 1 -1 . From (iii), ṽ (1)
1 and ṽ(1,2)

1 are independent. We deduce that P(ṽ

(1) 1 = ṽ(1,2) 1 ) = P(ṽ (1)
1 = -ṽ

(1,2) 1
) = 0. Summarizing, for our fixed permutation Q we have

P(QAQ T B = BQAQ T ) ≤ P(ṽ (1) 1 = ṽ(1,2 )
1
) + P(ṽ

(1) 1 = -ṽ (1,2) 1 ) = 0.
Then, the stated result follows from a union bound over the (finite) set P n .

Using Lemma 5, we conclude that P(E 3 ) = 0.

Concluding the proof. We have obtained

P(E) = P(E 1 ) + P(E 2 ) =0 + P(E 3 ) =0 ≤ n 2 n-1 ,
which concludes the argument.

Noiseless EMD dynamics

We now consider the "noiseless" version of the GM problem, where σ = 0 or equivalently A = B (given that we assumed X * = Id), which makes the gradient ∇E(X) that appears in (2.4) equal to A 2 X + XA 2 -2AXA. We consider the CGW model, or equivalently, when A ∼ GOE(n). Interestingly, in this case one iteration of Algorithm 2 suffices to recover the ground truth (i.e., we can fix N = 1). Theorem 2. Let A ∼ GOE(n) and let X (1) be the first iterate of the dynamic defined by (2.4) with initialization X (0) = J/n 2 . Then GMWM(X (1) ) = Id a.s., for any γ 0 > 0.

To prove the previous theorem, we will need the following auxiliary result. Lemma 6. Let A ∼ GOE(n). Then for any i, j ∈ [n] with i ̸ = j, -∇E(X (0) ) satisfies a.s. the property (2.2), i.e. the following holds almost surely

∇E(J) ii + ∇E(J) jj -2∇E(J) ij < 0.
Proof. We use the decomposition

∇E(J) ii + ∇E(J) jj -2∇E(J) ij = M ii + M jj -2M ij + M ′ ii + M ′ jj -2M ′ ij
where M := A 2 J + JA 2 and M ′ := -2AXA. Define a k as the k-th column of A, with this we have for i, j ∈ [n]

(A 2 J + JA 2 ) ij = ⟨a i , n k=1 a k J kj ⟩ + ⟨ n k=1 a k J ik , a j ⟩ = ⟨ n k=1 a k , a i + a j ⟩, which clearly implies M ii + M jj = 2M ij .
For the term M ′ we have

M ′ ii + M ′ jj -2M ′ ij = -⟨(a i -a j ), J(a i -a j )⟩ =(⟨a i -a j ,1⟩) 2 . Since ⟨a i -a j , 1⟩ = k / ∈{j,i} [A ik -A jk ], it is clear that ⟨a i -a j , 1⟩ ̸ = 0 almost surely.
Observe that in this case, it follows from its definition that X (1) is symmetric. Indeed, it is the Hadamard product of X (0) and exp ⊙ -γ 0 ∇E(X (0) ) , and both are symmetric matrices.

Proof of Theorem 2. Take i, j ∈ [n]. From (2.4) and X (0) = J/n 2 we deduce

X (1) ij = N -1 0 e -γ 0 ∇E J n 2 ij = N -1 0 e -γ 0 n 2 (A 2 J+JA 2 -2AJA) ij
where the scalar random variable N 0 = ∥e -γ 0 n 2 (A 2 J+JA 2 -2AJA) ij ∥ 1,1 is clearly almost surely strictly positive (we will only need this fact in the sequel). Define the following quantity q ij = X

ii X

(1) jj

X (1) ij ,
which is well-defined given that X

(1) ij > 0, almost surely. On the other hand,

q ij = exp (-γ 0 (∇E(J) ii + ∇E(J) jj -2∇E(J) ij )).
By Lemma 6 we deduce that for i ̸ = j, we have almost surely that q ij > 1. Which means that

X (1) ij < X (1)
ii X

jj , so we deduce that X

(1) ij < 1 2 (X (1)
ii + X

(1) jj ). Consequently, X (1) satisfies (2.2), and given the symmetry of X (1) , we deduce that GMWM(X (1) ) = Id almost surely.

Remark 2 (The Erdős-Rényi case). Clearly, if we take an Erdős-Rényi graph A ∼ G(n, p), the argument in Lemma 6 will not work, since there are, with probability one, at least two vertices (say i, j ∈ [n]) with the same degree, hence ⟨a i -a j , 1⟩ = 0. Using the degree profile distribution of Erdős-Rényi graphs, which is a well-studied problem [START_REF] Bollobás | The distribution of the maximum degree of a random graph[END_REF][START_REF] Bollobás | Degree sequences of random graphs[END_REF], one can give a bound on the number of vertices that will be correctly assigned by Algorithm 2 with N = 1. We remark that for the graph isomorphism problem, more elaborated "signatures" using the degrees have been used with success under the Erdős-Rényi model [START_REF] Babai | Random graph isomorphism[END_REF][START_REF] Czajka | Improved random graph isomorphism[END_REF][START_REF] Kiyavash Osman Emre Dai | Analysis of a canonical labeling algorithm for the alignment of correlated erdős-rényi graphs[END_REF][START_REF] Ding | Efficient random graph matching via degree profiles[END_REF].

Given a similarity matrix X, several existing theoretical results for exact recovery (see [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF][START_REF] Yu | Graph matching with partially-correct seeds[END_REF][START_REF] Mao | Exact matching of random graphs with constant correlation[END_REF][START_REF] Araya | Seeded graph matching for the correlated Wigner model via the projected power method[END_REF]) use the following sufficient condition for ensuring GMWM( X) = Id;

X ii > max j̸ =i X ij , ∀i ∈ [n], (3.8)
this is known as diagonal dominance. It is clear that (3.8) implies (2.1) (which is equivalent to (2.2) in the symmetric case). The approach of using (3.8) instead of (2.2) to prove Theorem 2 fails. Indeed, when X = X (1) as in Theorem 2, showing (3.8) is equivalent to showing ∇E(J) ii < ∇E(J) ij , for all j ̸ = i. A simple calculation shows that

E[∇E(J) ii ] = -1 + o(1), E[∇E(J) ij ] = o(1) ∀i ̸ = j, and Var[∇E(J) ij ] = 3 -o(1) ∀i, j.
Then using Chebyshev's inequality, we deduce that the event ∇E(J) ii ≥ ∇E(J) ij occurs with at least a constant probability, for any given i ̸ = j. Hence it follows that (3.8) will fail to hold with at least a constant probability.

Alternative analysis for GRAMPA in noiseless setting

To further illustrate the advantage of using (2.2) over (3.8), we will prove a stronger version of [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Lemma 2.3] which provides guarantees for exact recovery (in the noiseless case) for the GRAMPA algorithm that operates as follows. First, for a regularization parameter η, it computes the similarity matrix

Xgrampa := n i,j=1 1 η 2 + (λ i -λ j ) 2 v i v T i Jv j v T j . (3.9)
Here, A = n i=1 λ i v i v T i . Next, the algorithm rounds Xgrampa using any rounding procedure (e.g., linear assignment, greedy rounding etc.). For the result in [26, Lemma 2.3], however, conditions are derived on η which ensure that Xgrampa is diagonally-dominant 4 , with high probability.

In the following theorem, we show for the noiseless setting that exact recovery is in fact possible for any choice of η. Theorem 3. Let A be a GOE(n) distributed matrix with spectral expansion A = n i=1 λ i v i v T i . Then for any η, the matrix Xgrampa , defined in (3.9), satisfies GMWM( Xgrampa ) = Id, almost surely.

Proof. Define the matrix Λ η , with entries

(Λ η ) ij = (λ i -λ j ) 2 + η 2 .
Observe that Xgrampa can be rewritten as

Xgrampa = V Λ (-1) η ⊙ (V T JV ) V T , (3.10)
where Λ

(-1) η denotes the entrywise inverse of Λ η and V is the matrix with columns v 1 , . . . , v n . Recall that a matrix M ∈ R n×n is conditionally p.s.d (resp. conditionally p.d) if, for all v ̸ = 0 such that v T 1 = 0, we have v T M v ≥ 0 (resp. v T M v > 0); see e.g., [START_REF] Horn | Topics in matrix analysis[END_REF]Section 6.3]. We now have the following claim.

Claim 2. -Λ η is conditionally p.s.d for all η ∈ R.

Proof. We will first prove that -Λ 0 is conditionally p.s.d. For that, notice that

Λ 0 = (λ ⊙ λ)1 T + 1(λ ⊙ λ) T -2λλ T , where λ = (λ 1 , . . . , λ n ) T . Take v ∈ R n such that v ̸ = 0 and v T 1 = 0. Then v T (-Λ 0 )v = 2(v T λ) 2 -v T (λ ⊙ λ)1 T v + v T 1(λ ⊙ λ)v =0 ≥ 0.
On the other hand, it is evident that Λ η = Λ 0 + η 2 J. Since adding a scalar times the all ones matrix preserves the conditional p.s.d property, we deduce that -Λ η is conditionally p.s.d for any η.

Claim 2 is useful in light of Theorem 5 (see Appendix E), which says that the entrywise exponential of a conditionally p.s.d matrix is p.s.d. We will use this to prove the following claim. Proof. We will first prove that exp ⊙ (-tΛ η ) is p.d, for all t > 0. For that, using Claim 2 we obtain that -tΛ η is conditionally p.s.d. for all t ≥ 0. By Theorem 5 part (iii), exp ⊙ (-tΛ η ) is p.d. if and only if -tΛ η satisfies (2.2). It is clear that for all i ̸ = j and t > 0 is almost surely p.d and, on the other hand, it is evident that V T JV is p.s.d. In addition, the diagonal entries of V T JV (observe that (V T JV ) ii = (v T i 1) 2 ) are almost surely different from zero, since v i 's are uniformly distributed on the sphere (see e.g., [START_REF] Anderson | An Introduction to Random Matrices[END_REF]). Then by Lemma 11 we obtain that Λ log n for some constant c. In addition, our approach also leads to a shorter and arguably simpler proof.

-t(Λ η ) ii -t(Λ η ) jj =-2tη 2 > -2t(Λ η ) ij =-2t(λ i -λ j ) 2 -2tη 2
(-1) η ⊙ (V T JV ) is
Remark 4 (Extension to the Erdős-Rényi case). Notice that in the proof of Theorem 3, the assumption that A ∼ GOE(n) is crucial at two key points. First, to establish Claim 3, it is necessary that A almost surely possesses a simple spectrum. Second, we rely on this assumption to demonstrate that the quantities v T i 1 are almost surely non-zero for i ∈ [n]. In light of [START_REF] Tao | Random matrices have simple spectrum[END_REF]Corollary 1.4], we believe that it should be relatively straightforward to extend Theorem 3 to the case where A ∼ G(n, p), at least in a high probability sense.

Numerical experiments

We now provide numerical experiments using both synthetic and real data to assess the performance of Algorithm 2. We evaluate its performance along with other competing graph matching algorithms, considering both accuracy and running time as key metrics. In graph matching, the customary measure of accuracy is the overlap with the ground truth X * ∈ P n , also known as the recovery fraction. This measure, denoted as overlap(X * , X P ), is defined as: overlap(X * , X P ) := 1 n ⟨X * , X P ⟩ F , where X P represents the output of any graph matching algorithm, and we recall that ⟨X * , X P denotes the Frobenius inner product between X * and X P .

All the experiments were conducted using MATLAB R2021a (MathWorks Inc., Natick, MA) on an Apple M1 machine with 8 cores, 3.2 GHz clock speed, and 16GB of RAM. The reported running times are based on this setup.

Synthetic data setup

For the synthetic data experiments, we follow this general setup: we generate matrices A and B using the CER and CGW models described in Section 2.3, with X * chosen uniformly at random from P n . We then apply various graph matching methods based on convex relaxations, using A and B as inputs.

Standardization of the inputs. To ease the comparison between different models, we will consider their standardized version. Following [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF], we consider the following standardized version of the CER model described in Sec. 2.3

à = (A -E(A)) p(1 -p)n , B = (B -E(B)) p(1 -p)n ,
where A, B ∼ G(n, σ, p, X * ). This can be regarded as a simple preprocessing step in the pipeline of the three graph matching algorithms that will be experimentally tested here.

Step-size selection. An obvious question for any gradient descent-like method is the selection of a step-size updating rule. We investigate the following popular strategies for EMD and PGD.

• Fixed step-sizes. In [11, Theorem 9.16 and Example 9.17] it is proven that MD converges

in O 1 √ N
iterations, under the following constant step-size rule

γ md k = √ 2 log nL E,∞ √ N + 1 , ∀k ∈ [N ], (4.1)
where L E,∞ := max

X∈∆ n 2 ∥∇E(X)∥ ∞ and for v ∈ R d , with d ∈ N, ∥v∥ ∞ := max j∈N |v j |. Given
that L E,∞ is, in general, not observed, for each k we replace it (as a heuristic) with the following observed quantity

L(k) E,∞ = max X∈{X (0) ,...,X (k) } ∥∇E(X)∥ ∞ . (4.2)
Notice that, contrary to (4.1), the rule obtained by replacing L E,∞ by the empirical quantity L(k) E,∞ is not a fixed step strategy. Analogously, for the PGD algorithm, based on the same theoretical result we define the following update rule

γ pgd k = √ 2 L(k) E,2 √ N + 1 , ∀k ∈ [N ], (4.3)
where L(k) E,2 = max

X∈{X (0) ,...,X (k) } ∥∇E(X)∥ F . (4.4)
• Dynamic step-sizes. The O 1 √ N convergence of EMD and PGD has also been established, see [START_REF] Beck | First Order Methods in Optimization[END_REF]Theorem 9.18], under the following alternative strategy

γ md k = √ 2 ∥∇E(X (k) )∥∞ √ k+1 if ∇E(X (k) ) ̸ = 0 0 if ∇E(X (k) ) = 0, (4.5)
Notice that here γ md k can change from iteration to iteration, as it is adapting to the current gradient information. In the case of PGD the analogous rule is given by

γ pgd k = √ 2 ∥∇E(X (k) )∥ F √ k+1 if ∇E(X (k) ) ̸ = 0 0 if ∇E(X (k) ) = 0. (4.6)
In the case that ∇E(X (k) ) = 0, we have chosen to impute γ md k and γ pgd k with a zero value, but this choice is arbitrary and has no effect in the value of the iterates.

• Other heuristic rules. We also consider certain purely heuristic rules, which are not based on theoretical convergence results. In our experiments, these rules sometimes outperform the theoretically established rules described above. For PGD we consider

γ k = θ ∥∇E(X (k) )∥ 2 F ∥E(X (k) )∥ 2 F , (4.7)
where θ is a positive constant, which is either fixed to 1 (default value) or determined by line search.

Comparison with other convex optimization-based methods

We compare some of the state-of-art seedless graph matching methods based on convex relaxations with our proposed method, Algorithm 2, using a fixed number of iterations N = 125. In terms of efficient methods for convex graph matching, two of the better performing ones are Grampa [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF] and QPADMM, the latter of which solves (1.5) with K the set of doubly stochastic matrices, using the method of Alternating Direction Method of Multipliers (ADMM). The QPADMM algorithm has also been considered in the numerical experiments sections in [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF][START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF], under the name QS-DS. We recall that Grampa is a regularized spectral method, which adds the term η∥X∥ 2 F to the objective function E(•) in (1.5) with K = {1 T X1 = n}. Unless otherwise specified, we adopt the recommended value of η = 0.2 (see [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Section 4]). In our experiments, we observed that variations in this parameter have minimal impact on the overall performance of Grampa, which aligns with findings reported in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]. We also include the classic Umeyama algorithm, which was originally introduced in [START_REF] Umeyama | An eigendecomposition approach to weighted graph matching problems[END_REF], and corresponds to an unregularized spectral method. Figure 1: Comparison of the performance of seedless convex methods for CGW (Fig. 1a) and CER (Fig. 1b) models. We plot the average overlap with the ground truth (recovery fraction) over 15 Monte Carlo runs for graphs of size n = 300. We used N = 125 iterations in EMDGM and in Grampa we use the regularization parameter η = 0.2. For EMDGM we used the dynamic step-size rule (4.5).

Accuracy comparison. Figure 1 illustrates that, in terms of accuracy, the performance of EMDGM sits between Grampa and QPADMM, while Umeyama is outperformed by all three methods. While this result is not unexpected, considering the well-defined hierarchy of relaxations from the hyperplane to the simplex and, finally, to the Birkhoff polytope, these experiments aim to provide a clearer understanding of the gap between them. In this regard, taking the case of the CGW model in Fig. 1a, we see that EMDGM achieves exact recovery for noise levels as high as σ ≈ 0.5, while Grampa can recover only up to σ ≈ 0.25. On the other hand, QPADMM outperforms all the rest, with a recovery of up to σ ≈ 0.6. These experiments show that EMDGM (with N = 125) is closer, in terms of accuracy, to the tightest relaxation of QPADMM. As we will see in the sequel, the performance of EMDGM is highly correlated with the number of iterations. In the case of the CER model, the gap between methods reduce and even QPADMM can recover only up to σ ≈ 0.4. This confirms that matching CER graphs is more challenging for comparable levels of noise. One possible explanation for this fact is that the (continuous) weights in the CGW model help to discriminate the correct from incorrect matches ("breaking the symmetry") more effectively than in the binary weights case, where one relies solely on the graph structure. From the theoretical point of view, we saw a manifestation of this phenomenon in Section 3, when trying to extend Theorem 1 to the CER case (see Remark 2). Overall, these experiments show that EMDGM is closer, in terms of accuracy, to the tightest relaxation of QPADMM.

Running time comparison.

In Table 1, we summarize the average running time and standard deviation over the 15 Monte Carlo runs, for each algorithm. In these experiments, we observe that both Grampa and Umeyama are of comparable speed, while EMDGM operates about eight times slower than them. Conversely, QPADMM is the slowest, being about a thousand times slower than Grampa.

In the next subsection, we will explore how it is possible to achieve good accuracy with EMDGM using fewer iterations, thus speeding up the overall process.

Comparison with PGD

As mentioned earlier, the EMD algorithm is one of several algorithmic alternatives to solve (1.6).

In this section, we have two objectives: Firstly, to compare it with other algorithms that achieve the same task, and secondly, to evaluate its performance with different values of N . For runtime efficiency, we confine our comparison to the category of first-order methods, among which PGD is arguably the most popular one. We will call PGDGM the algorithm that solves (1.6) using a gradient descent approach, receiving the same inputs and following the same pipeline as Algorithm 2, but changing the update rule from (2.4) to (2.3). Given the existing theoretical convergence guarantees for MD and PGD, see [START_REF] Beck | First Order Methods in Optimization[END_REF]Sections 8 and 9], we know that there exists rate sequences (γ k ) k≥0 such that both methods will converge to a minimizer of (1.6). Our objective here is to provide a practical comparison of both algorithms within a unified setup, focusing on the number of iterations and the selection of step-size rules. In this section, we will demonstrate that the overall results of both approaches often diverge on synthetic data (a pattern consistent with real data, as we demonstrate later in Section 4.2 below).

Accuracy comparison. In Fig. 2 we plot the average overlap between X * and the output of EMDGM, PGDGM and Grampa. Among the step-size rules described in Section 4.1 above, the better performing ones in our experiments were the dynamic adaptive rule (4.5) for EMDGM and the heuristic rule for PGDGM with θ close to 1 (we chose to include the default value θ = 1 in the plot). We see that all three algorithms perform better under the CGW model than under the CER model. Indeed, one can see in Fig. 2a that EMDGM and PGDGM can achieve perfect recovery for σ ≈ 0.45, and Grampa for σ around 0.25. Interestingly, in the case of the CER model, we see that EMDGM performs better in terms of perfect recovery with a gap of at least 0.1 in the noise level it tolerates. When we restrict the number of the iterations to N = 25, we see that all models have an accuracy of 90% or more under a noise level as high as σ = 0.25. On the other hand, as seen in Fig. 2b, for the same noise level, only the EMDGM algorithm achieves 90%, under the setting of this experiment. Overall, EMDGM outperforms both Grampa and PGDGM in this experiments.

Running time comparison. Table 2 contains the running time information for the experiments shown in Figure 2. We see that Grampa is the fastest algorithm, running approximately five times as fast as EMDGM (with N = 25), which in turn is approximately two times faster than PGDGM. Despite EMDGM being slower than Grampa, we demonstrate here that when terminated at N = 25, EMDGM still outperforms Grampa, achieving an average speed five times faster than the case with N = 125.

Experimental comparison between (2.1) and diagonal dominance

As previously noted, both EMDGM and Grampa adhere to the meta-strategy outlined in the Introduction: their output is the rounded version of a similarity matrix (specific for each algorithm).

To prove exact recovery, a large part of the literature [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF][START_REF] Yu | Graph matching with partially-correct seeds[END_REF][START_REF] Mao | Exact matching of random graphs with constant correlation[END_REF][START_REF] Araya | Seeded graph matching for the correlated Wigner model via the projected power method[END_REF] has focused on proving that the similarity matrix X satisfies the following diagonal dominance property with high probability (under different random graph model assumptions and assuming that X * = Id), which we recall is

X ii > max j̸ =i X ij , ∀i ∈ [n].
Clearly, (3.8) implies that GMWM( X) = Id. We compare (3.8) with the property (2.1), which is weaker than (3.8) and also sufficient to ensure GMWM( X) = Id (as proven in Lemma 2 part (ii)).

The objective is to see if for moderately small values of n, we see the diagonal dominance of the similarity matrix. The setup is as follows: we run EMDGM (N = 25) and Grampa with n = 500 using the CGW model and compute the following two metrics for the similarity matrix of each method.

• Number of non-diagonally dominant rows. We consider the quantity

|{i ∈ [n] : X ii ≤ max j̸ =i X ij }| (4.8)
• Number of pairs not satisfying (2.1). In this case, we consider

|{(i, j) ∈ [n] 2 : i ̸ = j and X ii ∨ X jj ≤ X ij }| (4.9)
In Fig. 3 we show the heatmaps for the metrics (4.8) and (4.9), for Grampa and EMDGM algorithms.

For both algorithms, we do not see diagonal dominance in any of the Monte Carlo samples, however in the low noise regime, both algorithms have perfect recovery (as seen in Fig. 2c). If one compares both algorithms by the metric (4.8), the conclusion would be that Grampa performs slightly better than EMDGM, however by the results in Fig. 2c, we know the opposite is true. On the other hand, we see in Figs. 3c and 3d that the metric (4.9) is 0 for algorithms for values of σ up to around 0.2, which coincides with the results in Fig. 2c. In addition, under this metric, EMDGM outperforms Grampa, which is also the conclusion of Section 4.1.2. Overall, in this experiment, metric (2.1) seems to be much better correlated with the perfect recovery cases in both methods.

Real data

This section focuses on evaluating the performance of the EMDGM algorithm on real data in comparison to PGDGM and Grampa. Before proceeding, it is important to note a few general remarks. Firstly, since some of the datasets are relatively large (approximately ten thousand vertices per graph) and sparse, we have omitted the standardization preprocessing discussed in Section 4.1 to more effectively leverage their sparsity. Secondly, unless specified otherwise, we employ dynamic step sizes derived from (4.5) for EMDGM and (4.7) for PGDGM. This choice is based on their performance with synthetic data. The Grampa regularizer η is fixed to 0.2. and3d). The input graphs were generated using the CGW model with n = 500 and the noise level corresponds to the parameter σ. Recall that here a value of 0 means perfect recovery on that particular realization.

.

Computer vision dataset

We utilize the SHREC'16 dataset [START_REF] Lähner | Matching of Deformable Shapes with Topological Noise[END_REF], comprising 25 shapes (in both high and low resolution) representing a child in various poses. The goal is to match corresponding body parts when presented with two different shapes. The initial step in our pipeline involves converting the triangulated representation of an image (as stored in the database) into an adjacency matrix, a common procedure in computer vision (refer to [57, Section 1]). We opt for the use of low-resolution versions of each image, as they are already represented by adjacency matrices with an average size of approximately 10, 000 vertices. It is worth noting that, on average, the number of edges in these matrices is only 0.05% of the total possible edges, indicating a high degree of sparsity. To leverage this sparsity, we will omit the input standardization described in Section 4.1. Another important aspect of the SHREC'16 dataset is that all images have varying sizes, resulting in inputs of the form

A ∈ R n A ×n A and B ∈ R n B ×n B with n A ̸ = n B .
It is worth noting that both EMDGM and PGDGM naturally extend to accommodate this variation, and the same can be said for Grampa (see [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Section 4]). We will apply and compare the EMDGM, PGDGM, and Grampa algorithms using the Princeton benchmark protocol [START_REF] Kim | Blended intrinsic maps[END_REF]. This benchmark protocol essentially involves a CDF comparison, defined by the following steps.

• Normalized geodesic errors. Given inputs A and B, and the output π of a graph matching algorithm, we calculate the normalized geodesic error as follows: for each node i in A, we compute d B (π(i), π * (i)), where π * represents the ground truth matching, and d B denotes the geodesic distance on B (computed as the weighted shortest path distance using the triangulation representation of the image [START_REF] Gabriel | Numerical mesh processing[END_REF]). The normalized error is defined as

ε(i) := d B (π(i), π * (i))/ Area(B),
where Area(B) is the surface area of B (computed using the triangulation representation).

• CDF computation. The error cumulative distribution function (CDF) is defined for t ≥ 0 as

CDF(t) = n A i=1 1 ε(i)≤t , (4.10)
where n A is the number of nodes in A. Notice that CDF(0) corresponds to the recovery fraction, or overlap as defined in Section 4.1.

In Fig. 4 we show the average performance of the three algorithms on the SHREC'16 dataset. We see that the performance of EMDGM is the best of the three, while for a small value of ϵ (around 0.05), Grampa has a similar performance. It is interesting that for the other values of ϵ mirror descent has a big gap in performance compared to Grampa, because in the graph matching problem we are not trying to explicitly optimize the errors beyond ϵ = 0 (which by the definition (4.10) is related to exact recovery). In this case, EMDGM appears to introduce an additional level of regularity beyond what is specified in the optimization problem. We attribute the suboptimal performance of PGDGM to our difficulty in properly tuning the step-sizes.

Autonomous systems

We use data from the Autonomous Systems dataset, which belongs to the University of Oregon Route Views Project [START_REF]Autonomous systems peering networks[END_REF], and is also available at the Stanford Network Analysis Project repository [START_REF] Leskovec | Snap datasets: Stanford large network dataset collection[END_REF][START_REF] Leskovec | Graphs over time: Densification laws, shrinking diameters and possible explanations[END_REF]. This data has also been used to evaluate the performance of Grampa and other graph matching algorithms in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Section 4.5]. The dataset consists of a dynamically evolving network of autonomous systems (a collection of IP networks under a common routing policy) that exchange traffic through the Internet. There are a total of nine graphs, with each one representing weekly observations of the autonomous system from March 31, 2001, to May 26, 2001. The node count varies between 10, 670 and 11, 174, while the edge count fluctuates between 22, 002 and 23, 409 over time. Since the number of nodes and edges exhibits these fluctuations over time (see [START_REF] Leskovec | Snap datasets: Stanford large network dataset collection[END_REF] for more summarizing statistics), we have implemented a subsampling strategy adopted in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]Section 4.5]. We focus on a subset of 1, 000 vertices with the highest degrees, found in all nine networks. The nine resulting graphs can be seen as perturbations of each other. Fixing the first graph (from March 31), as input A, and comparing it to the subsequent ones, we observe that the correlation approximately decreases over time (although is not fully monotone). This provides a comparison framework similar to the synthetic data experiments, where the correlation was controlled by the parameter σ.

In Figure 5 we compare the output of the three methods EMDGM, Grampa, PGDGM under two different metrics. The first is the recovery fraction, plotted in Figure 5a and the other is the rescaled objective ⟨A, X T P B X P ⟩ F (here X P represent the output of a graph matching algorithm), plotted in Figure 5b. As noted in [START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF], the latter metric is more instructive for this dataset due to the presence of numerous degree-one vertices connected to high-degree vertices, resulting in symmetries that render the ground truth unidentifiable. In this experiment, we apply all three graph matching algorithms to nine instances created by fixing input A as the first graph. Input B is obtained by random node relabeling for each of the nine graphs. In Figures 5c and5d we plot the same metrics, but with standardized inputs. For that, we subtract to each matrix their average value and divide them by their standard deviation. Notably, the performance of EMDGM and PGDGM improves, while that of Grampa declines. This standardization alters the objective function and appears to introduce a 'regularization' effect for the methods on the simplex. Further investigation into this phenomenon is warranted, but we leave it for future research.

While we do not have a full explanation for why PGDGM outperforms EMDGM in this case, we can provide some partial insights. In certain scenarios, the "regularity" of the objective function may be more influential for achieving better results, whereas in other cases, the performance disparity may be primarily attributed to the selection of an appropriate step-size strategy. We believe that the latter plays a greater role in this case, based on the following experiment. We fix the step-size rule, using (4.1) for MD and (4.3), and define the efficiency ratio of the function E (as seen in [START_REF] Beck | First Order Methods in Optimization[END_REF]Example 9.17]) as:

ρ E = √ log n, L E,∞ L E,2 ,
where it is clear that

√ log n √ n ≤ ρ E ≤ √ log n.
When ρ E is closer to the lower bound, MD tends to outperform PGD, while the reverse is true when ρ E is closer to the upper bound. We compute an empirical version of the unobserved quantities L E,∞ and L E,2 (which can be regarded as regularity measures for E), by taking 10, 000 uniformly distributed samples on ∆ n 2 and computing the maximum of the ℓ ∞ and ℓ 2 norms of the gradient evaluated on those points. The value of ρ E obtained varies from 0.46 to 0.53, while

Facebook networks

We use Facebook's friendship network data from [START_REF] Traud | Social structure of facebook networks[END_REF], consisting of 11, 621 individuals affiliated with Stanford University. In contrast to the previous examples, where multiple networks were available, here, we have a single large network. As a result, we need to employ a subsampling strategy to generate correlated graphs for evaluating the graph matching methods. In the context of seeded graph matching, [START_REF] Yu | Graph matching with partially-correct seeds[END_REF] employed a specific sampling strategy to evaluate methods based on seed quality. However, since our objective is different, our sampling strategy also differs from theirs. The complete dataset contains 136, 660 (unweighted) edges, with a degree distribution closely resembling a power-law pattern (density proportional to d -1 ), and an average degree of approximately 100. We obtain two correlated graphs A, B as follows. First, we create graph H by selecting 1, 000 vertices uniformly at random, forming the induced subgraph. From the parent graph H, we retain each edge with a probability of s to create graph A. We independently repeat this process to generate graph B (as usual, we shuffle the labels of B uniformly before applying the graph matching algorithms). Intuitively, a higher value of s corresponds to a stronger correlation between the graphs, with s = 1 implying isomorphism. Notice that our sampling scheme closely resembles the "parent graph" definition in [START_REF] Pedarsani | On the privacy of anonymized networks[END_REF] for the CER model. In Figure 6 we plot the recovery fraction for EMDGM, PGDGM and Grampa. As expected, all three methods achieve exact recovery in the isomorphic case s = 1. However, for smaller values of s, there is a significant gap in performance between EMDGM and PGDGM (which are comparable) and Grampa. The more pronounced drop in performance, compared to the synthetic data experiments presented in Section 4.1, is also expected given the power-law behavior exhibited by this network. An interesting observation can be made from the confidence intervals in Figure 6: EMDGM shows less variance across the Monte Carlo subsamples. To quantify this, we computed the average standard deviation over the 6 instances of s. For EMDGM, this value is 0.0126, while for PGDGM, it is significantly higher at 0.0604.

Concluding remarks

In this work, we present a novel graph matching algorithm based on a convex relaxation of the permutation set to the unit simplex. We introduce an iterative mirror descent scheme that offers closed-form expressions for the iterates, resulting in enhanced efficiency compared to competing projecting gradient methods. Theoretical analysis demonstrates that, when coupled with a standard greedy rounding procedure, our method can exactly recover the ground truth matching in the noiseless case of the correlated Gaussian Wigner model. Experimental results support this claim, showing that our proposed algorithm can achieve exact recovery even in the presence of significant noise. In terms of statistical performance, it consistently outperforms other convex relaxation-based methods with similar time complexity.

The following are potential avenues for future research.

• An intriguing theoretical question is to extend the guarantees to the case where σ ̸ = 0 in the CGW model. This presents two main challenges. First, understanding the eigenvector structure of the iterates in (2.4) is nontrivial, given its Hadamard product definition and the complexity of the gradient expression. Second, using Lemma 2 becomes less straightforward when not all pairs satisfy the condition (2.1). To address this, one possible approach is to demonstrate that the number of pairs not satisfying (2.1), as given in (4.9), converges to zero. However, this seems to present significant technical challenges.

• Extending recovery guarantees to other models is also an important theoretical objective. A common approach in the literature of iterative methods involves using uniform concentration bounds, which is essentially a worst-case scenario approach. To follow this strategy, it will be crucial to gain a deeper understanding of the structure of the iterates in (2.4) beyond X (1) .

In particular, we need to improve our ability to "localize" X (k) . Otherwise, there is a risk that, in the worst-case scenarios, the gradient may point in the wrong direction. It is not hard to provide examples of such scenarios.

• From an experimental perspective, it is intriguing to explore acceleration in the mirror descent dynamic. Although accelerated MD schemes have been considered in the literature, as discussed in [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF], they involve Bregman projections that are fast in theory but lack a closed form, making their implementation more involved. It would be interesting to assess the practical performance of these schemes.

• Another experimental question involves finding a rate selection method for EMD that enhances its statistical performance, similar to the heuristic (4.7) in the case of PGD. In general, any other techniques used to improve the performance of first-order methods can also be tested here, including adaptive restart, early stopping, and more.

for all j ≤ k -1. This in turn implies that r k ≥ e -2γ 0 /n 2 k-1 j=1 e -4γ j /n(n-1) ≥ e Given that the sequence of rates satisfies

n(n -1) 4 log 2 > N -1 j=0 γ j ≥ k-1 j=0 γ j ,
we have that

r k > 1 2 ≥ a σ -2 a σ ,
where the last inequality follows from the fact that σ ≤ 1. Noticing that k was arbitrary in [N -1], we conclude that r j > aσ-2 aσ , for all j ∈ [N -1]. From this, and from the fact that γ j > 0, we deduce that for all 2 ≤ j ≤ N -1 we have that e -γ j aσ r j -aσ +2 n(n-1)r j +n < 1, which combined with (A.7) gives that r k < 1, for 2 ≤ k ≤ N (the fact that r 1 < 1 is evident).

B Proof of Lemma 1

Notice that for any Π ∈ P n and X ∈ R n×n we have

∥AX -XΠBΠ T ∥ 2 F = ∥AXΠ -XΠB∥ 2 F ,
by the unitary invariance of the Frobenius norm. Then, given the permutation invariance of the set K, if X ∈ S(A, B), then XΠ T ∈ S(A, ΠBΠ T ). It remains to show that if XP is the output of the algorithm described in the lemma with input A, B, then the XP Π T will its output for the inputs A, ΠBΠ T . For this is sufficient to notice that (given the mechanism of Algorithm 1) that if the largest element of X is the entry (i, j) then the largest of XΠ T is the entry (i, π(j)), where π is the permutation map associated to Π.

C Proof of Lemma 4

Note that for all i ∈ [n], we have v i is marginally distributed uniformly on S n-1 . This is equivalent

to v i dist = g ∥g∥ 2 ,
where g is a standard normal vector in R n . It is easy to see that sign(v i (k)) dist = sign(g(k)), and the variables sign(g(1)), . . . , sign(g(n)) are i.i.d Rademacher distributed. Then, for any given i ∈ [n], we have

P(∀k, v i (k) > 0) = 1 2 n .
The result follows by applying the union bound.

  lk . With this, it is easy to verify that ÃB = B Ã ⇐⇒ Ãkl B ll = B kk Ãkl , ∀k, l ∈ [m], which implies that Ãlk Ãkl B ll = B ll Ãlk Ãkl , ∀k, l ∈ [m] with k ̸ = l, and Ãll B ll = B ll Ãll , ∀l ∈ [m]. In the sequel, we will assume w.l.o.g. that n 1 ≤ n 2 and only use the following implication ÃB = B Ã =⇒ Ã12 Ã21 B 11 = B 11 Ã12 Ã21 , and Ã11 B 11 = B 11 Ã11 .

  p.d for all η ∈ R.

  almost surely, since λ i ̸ = λ j almost surely [4, Theorem 2.5.2]. Hence, exp ⊙ (-tΛ η ) is almost surely p.d. for all t > 0. Given that 1 s = ∞ 0 e -st dt, we have the representation Λ -tΛ η )dt,where the integral of a matrix is understood entrywise. From Lemma 12, we have the integral of positive definite matrices is positive definite. Given that exp ⊙ (-tΛ η ) is almost surely p.d. for all t > 0, we conclude that Λ

  (a) CGW model. (b) CER model with p = 1 2 .

  (a) N = 125 in the CGW model. (b) N = 125 in the CER model with p = 1 2 . (c) N = 25 in the CGW model. (d) N = 25 in the CER model with p = 1 2 .

Figure 2 :

 2 Figure 2:We compare EMDGM, Grampa and PGD under the CGW and CER models for n = 500. For EMDGM we used the dynamic step-size rule (4.5) and PGD we used the heuristic rule (4.7) with θ = 1. In both cases, we fix the number of iterations to N = 25 (Figs.2c and 2d) and N = 125 (Figs.2a and 2b). The Grampa regularization parameter γ is set to 0.2 as recommended in[START_REF] Zhou Fan | Spectral graph matching and the quadratic relaxation I: the Gaussian model[END_REF]. We present the average overlap (recovery fraction) for 25 Monte Carlo runs. The shaded area indicates where 90% of the data falls (excluding the top and bottom 5%).

  (a) Metric (4.8) for Grampa. (b) Metric (4.8) for EMDGM.(c) Metric (4.9) for Grampa.(d) Metric (4.9) for EMDGM.

Figure 3 :

 3 Figure3: Heatmap for the metrics (4.8) and (4.9) for the similarity matrices of Grampa (Figs.3a and 3c) and EMDGM with N = 125 (Figs.3b and 3d). The input graphs were generated using the CGW model with n = 500 and the noise level corresponds to the parameter σ. Recall that here a value of 0 means perfect recovery on that particular realization.

Figure 4 :

 4 Figure 4: Performance comparison of EMDGM, PGDGM, and Grampa on the SHREC'16 dataset, based on the error cumulative distribution function (CDF) defined in (4.10). We conducted 100 iterations each for EMDGM and PGDGM.

  (a) Recovery fraction in a high-degree subsampled graph. (b) Common edges in a high-degree subsampled graph. (c) Recovery fraction in a high-degree standardized graph. (d) Rescaled objective in a high-degree standardized graph.

Figure 5 :

 5 Figure5: Comparison of EMDGM, Grampa, and PGDGM using a high-degree subsampled graph with 1, 000 vertices for each time instant. The input graph A remains fixed as the graph from March 31, while the input graph B varies (with random permutations of its rows and columns), as indicated on the x-axis. We used the dynamic step rule for EMDGM, derived from (4.5) and the heuristic rule (4.7) for PGDGM, with 125 iterations in both cases.

Figure 6 :

 6 Figure6: Comparison of graph matching algorithms on subsamples of 1, 000 edges of a Facebook friendship network. We plot the average over 15 Monte Carlo (over the uniform subsamples). We used the dynamic step rule for EMDGM, derived from (4.5) and the heuristic rule (4.7) for PGDGM, with 125 iterations in both cases.

  ×n 2 is a matrix with independent Gaussian entries with

	variance 1 n . From [4, Prop. 4.1.3], we see that Ã12 Ã21 has n 1 different eigenvalues almost
	surely. Hence, its eigenvectors, denoted	ṽ(1,2) 1	, . . . ,	ṽ(1,2)

n 1

  almost surely p.d. From (3.10) it is easy to see that Xgrampa is almost surely p.d. Combining this with Lemma 2 part (iii), the conclusion follows.

	Remark 3 (Comparison with Lemma 2.3 in [26]). Theorem 3 is stronger than [26, Lemma 2.3]
	in two aspects. First, it provides an almost sure guarantee instead of the high probability bound
	in [26, Lemma 2.3]. Second, it holds for any η ∈ R, while in [26, Lemma 2.3] the result requires
	that 1 n 0.1 < |η| < c

Table 1 :

 1 Average (over 15 Monte Carlos) running time in seconds (standard deviation in parentheses).

		Umeyama	EMDGM	Grampa	QPADMM
	CGW 0.0356(0.0211) 0.8609(0.0282) 0.0490(0.0048) 34.2666(1.1065)
	CER	0.0382(0.015) 0.9226(0.0168) 0.0432(0.0063) 29.1306(1.1302)

Table 2 :

 2 Average (over 25 Monte Carlos) running time in seconds (standard deviation in parentheses).

Allowing E to take vectors in R n

as arguments, using the canonical matrix-vector identification.

Here a non-Euclidean metric means simply a metric different from ∥ • ∥2 = ⟨•, •⟩.

This set corresponds to S(A) introduced in Section 2. Here we drop the dependence on A to ease the notation.

In fact, they consider an even stronger version of diagonal-dominance than (3.8).

√log n √ n = 0.0831 and √ log n = 2.62. Using this metric, there is not a clear advantage between the two methods. However, when applying rule (4.3), we observe that PGDGM significantly underperforms EMDGM (approximately ten times worse). It appears that the heuristic rule (4.7) plays a crucial role in accelerating the convergence of PGDGM. Unfortunately, we currently lack a similar rule for EMDGM, but this is an avenue for future exploration.

A The population dynamics

We introduce and analyze the population (or average) version of the mirror descent dynamics. For that, notice that all randomness in the dynamics defined by (2.4) is given by the input matrices A, B. In addition, the dependence on these matrices in the iterative process is captured by the term ∇E(X) = A 2 X + XB 2 -2AXB (as it is clear from line 3 of Algorithm 2). We define the population dynamics, by replacing the (random) loss function E(X) by the (deterministic) population loss, defined as E A,B E(X). This is equivalent to replace the term ∇E(X) in (2.4) by its expectation ∇E(X) := E A,B ∇E(X), which is proven as part of Lemma 10. More formally, this dynamic, which is deterministic, is defined given a sequence (γ k ) ∞ k=1 of positive rates, by the recursion

(A.1)

For each k, we define the normalization factor by

which is a positive real number.

A.1 One iteration (N = 1)

The following result shows that Algorithm 2 succeeds in finding the ground truth X * in one step (i.e., when N = 1). We will assume without loss of generality that X * = Id, and consequently in the definition of the population gradient, we have∇E(X) = E A,B ∇E(X), where the joint law of A, B is W(n, σ, X * = Id). The following result shows that in the case of the population dynamic, one step is sufficient for recovering the ground truth if we applied Algorithm 1 directly after the first iteration.

Proposition 1. Let X (1) ∈ R n×n be the first iterate in the dynamic defined by (A.1) and (A.2), defined by γ 0 ∈ R + . Then GMWM(X (1) ) = Id, for any γ 0 > 0.

Proof. From (A.1) and (A.2), we have

where N 0 > 0 is the normalization term. Using the definition for the population gradient, Lemma 10 and by a simple calculation, we get

where we define

Given that e 2γ 0 n 2 > 1 for γ 0 > 0, we deduce that X (1) satisfies (3.8), which implies that GMWM(X (1) ) = Id.

A.2 Multiple iterations (N > 1)

In practice, we do not have a priori information on how many iterations we need to perform before the rounding step, but we expect Algorithm 2 to be relatively robust with respect to the choice of the number of iterations N (at least in some range), for some choice of the step size sequence (γ k ) N k=1 . The next result proves that this is the case for the population dynamic, in the sense that all its iterates can be rounded to the ground truth, provided that a certain condition on the rate sequence (γ k ) k≥0 holds.

Theorem 4. Consider X (0) , X (1) , • • • , X (N ) ∈ R n×n the iterates of the dynamic defined by (A.1) and (A.2), defined by the rate sequence

Before proving this theorem, we will prove two lemmas that characterize the structure of iterates in the population dynamic.

k=0 be a trajectory of the population dynamic given by some step size sequence

Proof. The proof is by induction. Clearly X (0) is symmetric. Assume that X (k) is symmetric. Notice that

The proof follows by the inductive hypothesis, by the fact that ∇E(X) is symmetric if X is symmetric (see (D.5)) and the fact that the Hadamard product of two symmetric matrices is symmetric. We also prove (A.4) by induction. The base case k = 0 is evident. We assume that X (k) satisfies the property (A.4). Using (A.5), the fact that ∇E(X) satisfies (A.4) if X does, and the fact that (A.4) is closed under Hadamard products, the conclusion follows.

The previous lemma asserts that each entry of X (k) can take only two values, depending on its location inside the matrix (diagonal or off-diagonal). In this case, it is easy to see that the greedy rounding of X (k) , for any k, will be equal to the identity (the ground truth) if and only if

ij , for all j ̸ = i. Since the iterates of the population dynamics X (k) , up to the step K (for a given K ∈ N), are completely determined by the sequence of learning rates (γ k ) K k=0 , we can characterize them as a function of them. To prove the diagonal dominance, it will be convenient to define the ratio between the off-diagonal and diagonal entries, for k ≥ 1, as follows,

.

Note that, in the previous definition, r k is a function from R k to R. The following lemma provides a recursive characterization for these ratios.

Lemma 8. For any given sequence (γ k ) ∞ k=0 , the ratio between off-diagonal and diagonal entries r k defined above, satisfy the following

where a σ := 2 + n+1 n σ 2 . Proof. To save notation, we will write r k instead of r k (γ 0 , • • • , γ k-1 ). Using Lemma 7 and (D.5) we get that

11 . From this and the definition of the population dynamics, we see that

for k ≥ 2. Similarly,

On the other hand, given the fact that X (l) ∈ ∆ n 2 we have, for any l ∈ N,

which together with the fact

, imply the following,

Plugging this on (A.6) finishes the proof.

Proof of Thm.4. Given Lemma 7, it suffices to compare X (k) [START_REF] Beck | First Order Methods in Optimization[END_REF] and

11 or, equivalently, that r k < 1, then X (k) will satisfy the diagonal dominance property on all of its rows, thus we will get GMWM(X (k) ) = Id. We now fix a k ∈ [N -1]. From Lemma 8, it is easy to see that

On the other hand, notice that a σ < 4 (because σ ≤ 1), from which we deduce that

,

D Useful calculations

We gather here some useful calculations, used at different parts of the analysis.

Lemma 9.

Let A be a GOE(n) distributed random matrix, then

where T is a n 2 × n 2 matrix that represents the transposition operator defined on R n×n .

Proof. Take i, j ∈ [n]. We have

from which we get (D.1). To prove (D.2), we take any matrix X ∈ R n×n and compute

and from this we deduce that

Given that (A ⊗ A) vec(X) = AXA, it is clear that the previous implies the expression (D.2).

Lemma 10. Let A, B ∼ W(n, σ, X * = Id). Recall the definition E(X) = ∥AX -XB∥ 2 F or, in vector notation, E(vec(X)) = vec(X) T H vec(X), where H = (Id ⊗A -B ⊗ Id) 2 . We have

and

where T ∈ R n 2 ×n 2 is the matrix representing the transpose transformation in R n×n . In particular, we have (in matrix language)

Proof. Both equalities in (D.3) follow from E A,B (E(vec(X))) = vec(X) T E A,B (H) vec(X) and linearity of expectation. On the other hand, a simple calculation shows that

From this and Lemma 9, the expression (D.4) follows. Equality (D.5) follows by translation of (D.3) and (D.4) to the matrix language.

E Useful facts about positive semidefinite matrices

The following result proves that the entrywise exponential of a conditionally positive definite matrix is positive definite. Recall that a n × n symmetric matrix M is said to be conditionally p.s.d (resp. conditionally p.d) if for all v ∈ R n×n such that v T 1 = 0 and v ̸ = 0 we have v T M v ≥ 0 (resp. v T M v > 0).