Luca Ganassali
email: luca.ganassali@

Laurent Massoulié
email: laurent.massoulie@

Jacob Abernethy

Shivani Agarwal

33rd

Keywords: graph alignment, tree matching, Erdős-Rényi random graphs

In this paper we consider alignment of sparse graphs, for which we introduce the Neighborhood Tree Matching Algorithm (NTMA). For correlated Erdős-Rényi random graphs, we prove that the algorithm returns -in polynomial time -a positive fraction of correctly matched vertices, and a vanishing fraction of mismatches. This result holds with average degree of the graphs in O(1) and correlation parameter s that can be bounded away from 1, conditions under which random graph alignment is particularly challenging. As a byproduct of the analysis we introduce a matching metric between trees and characterize it for several models of correlated random trees. These results may be of independent interest, yielding for instance efficient tests for determining whether two random trees are correlated or independent.

Introduction

Graph alignment consists in finding an injective mapping (matching) S ⊂ V (G 1 )×V (G 2 ) between the vertex sets of two graphs G 1 and G 2 such that, for any two matched pairs (i, u), (j, v) ∈ S, then occurrence of edge {i, j} in G 1 tends to correspond to occurrence of edge {u, v} in G 2 . When this correspondence is exact for all pairs (i, u), (j, v) of matches, then the subgraphs of G 1 , G 2 induced by their nodes appearing in S are isomorphic. In this sense, graph alignment is the search for approximate graph isomorphisms. It has many applications, among which: social network de-anonymization [START_REF] Narayanan | Robust de-anonymization of large sparse datasets[END_REF], [START_REF] Narayanan | De-anonymizing social networks[END_REF]), analysis of protein-protein interaction graphs [START_REF] Kazemi | PROPER: global protein interaction network alignment through percolation matching[END_REF], [START_REF] Feizi | Spectral alignment of graphs[END_REF]), natural language processing [START_REF] Bayati | Message-passing algorithms for sparse network alignment[END_REF]), medical image processing [START_REF] Lombaert | Diffeomorphic spectral matching of cortical surfaces[END_REF].

A recent thread of research [START_REF] Pedarsani | On the privacy of anonymized networks[END_REF], Dai et al. (2019b) [START_REF] Cullina | Improved achievability and converse bounds for Erdős-Rényi graph matching[END_REF], [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF], [START_REF] Cullina | Partial recovery of Erdős-Rényi graph alignment via k-core alignment[END_REF], Fan et al. (2019a), Fan et al. (2019b)) has investigated the fundamental limits to feasibility of graph alignment in the context of a natural generative model of correlated graphs, namely the correlated Erdős-Rényi random graph model ERC(n, p, s). Specifically, it consists of two random graphs G 1 and G 2 on node set [n] = {1, . . . , n} obtained as follows. First generate two aligned graphs G 1 , G 2 with adjacency matrices A 1 , A 2 such that for each node pair {i, j}, one has P A 1 (i, j) = A 2 (i, j) = 1 = ps, P A 1 (i, j) = 1, A 2 (i, j) = 0 = P A 1 (i, j) = 0, A 2 (i, j) = 1 = p(1s), P A 1 (i, j) = A 2 (i, j) = 0 = 1p(2s), and this independently over node pairs {i, j}. Graph G 2 's adjacency matrix is then A 2 = M σ A 2 (M σ ) , where M σ is the matrix associated to a permutation σ chosen uniformly at random from S n .

Researchers have strived to determine for which parameter values (p, s), assuming n 1, one can recover the unknown permutation σ, and therefore the alignment G 2 of G 2 with G 1 . As in other high-dimensional inference tasks such as community detection, one expects such goal to be either poly-time achievable, achievable though not in poly-time, or impossible to achieve. The corresponding regions of parameter space are usually referred to as the "easy", "hard" or "Informationtheoretically (IT) impossible" phases for the problem considered. [START_REF] Cullina | Improved achievability and converse bounds for Erdős-Rényi graph matching[END_REF] have shown that it is possible to recover σ if and only if npslog(n) -→ n→∞ +∞ , thereby characterizing the "IT-impossible" phase. [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF] have proposed a polynomial-time 'degree profile matching' algorithm, and proven it to recover σ under the conditions np ≥ log α (n), 1s ≤ log -β (n) for suitable constants α, β > 0, thereby identifying a subset of the "easy" phase. More recently, Fan et al. (2019a), Fan et al. (2019b) have proposed a spectral method, and proven it to recover σ under the same conditions.

The result of [START_REF] Cullina | Improved achievability and converse bounds for Erdős-Rényi graph matching[END_REF] shows that there is no hope of recovering σ, or in other words, of perfectly re-aligning G 1 and G 2 , in the case of sparse graphs, that is graphs with average degree np of order 1. Nevertheless, their result does not rule out the possibility of partially recovering the unknown permutation σ. For the applications mentioned earlier, it is at the same time natural to assume that the graphs involved are sparse, and potentially useful to recover only a fraction of the unknown matches (i, σ(i)).

Objectives and main result

This motivates the present work, whose goal is to show that partial alignment of sparse correlated graphs is feasible, and to introduce a polynomial-time algorithm for producing such partial alignments. Our main result is the proposal of the so-called Neighborhood Tree Matching Algorithm (NTMA hereafter) together with the following Theorem I Consider the correlated Erdős-Rényi model ERC (n, p, s), where p = λ/n. For some λ 0 > 1, for all λ ∈ (1, λ 0 ], there exists s * (λ) < 1 such that, provided s ∈ (s * (λ), 1], the NTMA returns a matching S verifying the following properties with high probability:

|S ∩ {(i, σ(i)), i ∈ [n]}| = Ω(n), |S \ {(i, σ(i)), i ∈ [n]}| = o(n).
(1)

In words, our algorithm returns a set of node alignments which contains a negligible fraction of mismatches, and Ω(n) good matches. Our result covers values of λ arbitrarily close to 1, and thus applies to very sparse graphs. For λ < 1, Erdős-Rényi graphs in our correlated model have connected components of size at most logarithmic in n, so that there is no hope of recovering a positive fraction of correct matches. This result can be interpreted as follows. For partial graph alignment of sparse Erdős-Rényi correlated random graphs, there is an "easy phase" that includes the parameter range {(λ, s)

: λ ∈ (1, λ 0 ], s ∈ (s * (λ), 1]}.

Paper organization

Description of the Neighborhood Tree Matching Algorithm and the proof strategy for establishing Theorem I are given in Section 2. Our algorithm relies essentially on a tree matching operation.

To pave the way for Section 2, we therefore introduce in Section 1 a notion of matching weight between trees that is key for our algorithm, and can be computed efficiently in a recursive manner.

We further obtain probabilistic guarantees on the matching weights between random trees drawn according to some Galton-Watson branching processes. These are instrumental in the proof of Theorem I. However these may be of independent interest. Indeed we introduce in Section 1 a natural hypothesis testing problem on pairs of random trees, for which we obtain a successful test based on computation of tree matching weights.

Related work

Most relevant to the present work are the papers [START_REF] Pedarsani | On the privacy of anonymized networks[END_REF], Dai et al.

(2019b) [START_REF] Cullina | Improved achievability and converse bounds for Erdős-Rényi graph matching[END_REF], [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF]), Fan et al. (2019a), Fan et al. (2019b) already mentioned, which also focus on graph alignment in the context of the correlated Erdős-Rényi model. The main differences between the present paper and these is our focus on sparse random graph models, with average degree λ = O(1), our treatment of correlation coefficients s bounded away below 1, and our aim of partial rather than full alignment. Article [START_REF] Cullina | Partial recovery of Erdős-Rényi graph alignment via k-core alignment[END_REF] addresses a notion of partial alignment stronger than ours, and hence requires conditions under which graphs are not sparse. [START_REF] Makarychev | Maximum quadratic assignment problem: Reduction from maximum label cover and lp-based approximation algorithm[END_REF] show that graph alignment is NP-hard to solve, even approximately. This justifies the search for custom algorithms in a variety of scenarios.

The main methods proposed are: Percolation methods based on some initial seeds, i.e. matched node pairs provided a priori [START_REF] Kazemi | PROPER: global protein interaction network alignment through percolation matching[END_REF]). We remark that the matchings returned by our algorithm could be used as seeds, and then processed e.g. using percolation matching to eventually obtain an improved matching. Spectral methods are considered in [START_REF] Feizi | Spectral alignment of graphs[END_REF]), Fan et al. (2019b); Degree profile matching is introduced in [START_REF] Ding | Efficient random graph matching via degree profiles[END_REF]; Quadratic programming approaches are proposed in [START_REF] Zaslavskiy | A path following algorithm for the graph matching problem[END_REF]. Message passing methods are introduced in [START_REF] Bayati | Message-passing algorithms for sparse network alignment[END_REF]. These are structurally similar to our neighborhood tree alignment approach, which is implemented in a recursive manner and can be seen as a message-passing method. Our algorithm is however different, and comes with novel theoretical guarantees. Graph alignment is a special case of the quadratic assignment problem, reviewed in [START_REF] Panos | The quadratic assignment problem: A survey and recent developments[END_REF]. Database alignment is an important variant of graph alignment, studied in Dai et al. (2019a).

Notations

For a graph G, denote by V (G) its set of vertices, E(G) its set non-oriented edges, and

- → E (G) := {(i, j), {i, j} ∈ E(G)} its set of oriented edges. We use the notations u ∼ v if {u, v} ∈ E(G) and u → v if (u, v) ∈ - → E (G). The usual graph distance in G will be denoted δ G . For v ∈ V (G), let N G (v) denote the neighborhood of v in G, and deg G (v) its degree. For d ≥ 1 we also define B G (v, d) the set of vertices at (graph) distance at most d from v, and S G (v, d) := B G (v, d) \ B G (v, d -1
) the set of vertices at distance d from v. For a rooted tree T , we let ρ(T ) denote its root node. For any i ∈ V (T ) \ {ρ(T )}, we let π T (i) denote the parent of node i in T . For d ≥ 1, we note B d (T ) = B T (ρ(T ), d) and L d (T ) = S T (ρ(T ), d). We omit the dependencies in G or T of these notations when there is no ambiguity.

Tree matching

In this section, we introduce the matching weight between rooted trees and the related matching rate. We then establish bounds on the matching rate for specific models of random trees. We also give an application to a hypothesis testing problem on the independence between two trees.

1.1. Matching weight of two rooted trees Definition 1.1 For any d ≥ 0, let A d denote the collection of rooted trees whose leaves are all of depth d. Given two rooted trees T and T , let M d (T , T ) denote the collection of trees t ∈ A d such that there exist injective embeddings f : V (t) → V (T ), f : V (t) → V (T ) that preserve the rooted tree structure, i.e. such that

f (ρ(t)) = ρ(T ), f (ρ(t)) = ρ(T ), ∀i ∈ V (t) \ {ρ(t)}, f (π t (i)) = π T (f (i)), f (π t (i)) = π T (f (i)).
The matching weight of T and T at distance d is then defined as:

W d (T , T ) := sup t∈M d (T ,T ) |L d (t)| , (2)
i.e. the size, measured in number of leaves, of the largest tree in M d (T , T ).

ρ(T ) ρ(T )
Figure 1: Example of two trees T , T with W 3 (T , T ) = 7, where an optimal t ∈ A 3 is drawn in red.

Recursive computation of W d

We shall need the following notations and definitions. For a tree T , for i ∈ V (T ), and d ≥ 0,

T (i) d
is the sub-tree of T re-rooted at i, containing all vertices at distance less than d of i. For i, j ∈ V (T ) such that j → i, T

(i←j) d

denotes the sub-tree of T re-rooted at i, containing all vertices at distance less than d of i but where vertex j has been removed. By definition two vertices not connected by a path are at distance ∞. T Definition 1.2 For a given pair of trees T and T , for all pair of vertices (i, u) ∈ V (T ) × V (T ), the matching weight of (i, u) at depth d is then defined as:

W d (i, u) := sup t∈M d T (i) d ,T (u) d |L d (t)| .
(3)

Moreover, for all pairs of vertices (i, u), (j, v) ∈ V (T ) × V (T ) such that j → i and v → u, the matching weight of edges j → i and v → u at distance d is then defined as:

W d (i ← j, u ← v) := sup t∈M d T (i←j) d ,T (u←v) d |L d (t)| . (4)
Remark 1.1 This definition is compatible with the first one in the context of tree matching: one has

W d (ρ(T ), ρ(T )) = W d (T , T ). Note that W 0 (i, u) = 1 and W 1 (i, u) = max (deg(i), deg(u)). Similarly, W 0 (i ← j, u ← v) = 1 and W 1 (i ← j, u ← v) = max (deg(i), deg(u)) -1.
Now fix T and T . From these definitions, for all d ≥ 1, (i, u), (j, v) ∈ V (T )×V (T ) such that j → i and v → u, by doing a first step conditioning, we obtain the following recursion formulae:

W d (i ← j, u ← v) = sup m∈M(N T (i)\{j} , N T (u)\{v}) (k,w)∈m W d-1 (k ← i, w ← u), (5)
where M (E, F) is the set of all maximal injective (or one-to-one) mappings m : E 0 ⊆ E → F, where maximal means that they are not restrictions of another one-to-one mapping m :

E 1 → F with E 1 such that E ⊇ E 1 ⊃ E 0 .
In the same way we have

W d (i, u) = sup m∈M(N T (i) , N T (u)) (k,w)∈m W d-1 (k ← i, w ← u). (6)
Thus matching weights at depth d can be obtained by computing weights at depth d-1 and solving a linear assignment problem (LAP). Recursive formulae ( 5) and ( 6) yield simple recursive algorithms (see Algorithms 3 and 2 in A) to compute all matching weights at depth d. However for small values of n, the recursive algorithms 3 and 2 are faster, although of complexity O d 2d max , which is not polynomial for d = Ω(log n).

1.3. Matching rate of random trees Definition 1.3 Consider two random trees T , T . Their matching rate is defined as

γ(T , T ) := inf γ : ∃m, c, d 0 > 0, ∀x ≥ 0, ∀d ≥ d 0 , P W d (T , T ) ≥ mxγ d ≤ e -(x-c) + . (7)
This quantity captures the geometric rate of growth of matching weights with depth d. A simpler alternative definition could be γ(T ,

T ) := inf γ : P W d (T , T ) ≥ γ d -→ d→∞ 0 .
However definition (7) better suits our purpose.

Remark 1.3 By definition, note that for any γ > γ(T , T ), P W d (T , T ) ≥ γ d converges to 0 very fast, like O exp -c(γ) d with c(γ) > 1, so that γ(T , T ) ≤ γ(T , T ).

Models of random trees

We now describe three models of random trees that are relevant to sparse graph alignment.

GW (λ):

We consider two independent Galton-Watson trees T and T with offspring distribution Poi(λ), λ > 0. We denote (T , T ) ∼ GW (λ).

GW (λ, s, δ): For δ ≥ 1, consider a labeled tree T rooted at ρ and a tree T rooted at ρ . ρ is also a node of T , at distance δ from its root ρ. The two trees are generated as follows. First, nodes in T on the path from ρ to the parent of ρ in T have, besides their child leading to ρ , Poi(λ) children in T , themselves having offspring in T given by independent Galton-Watson trees with offspring Poi(λ). Then, the intersection between T and T is a Galton-Watson tree with offspring Poi(λs), with λ > 0 and s ∈ [0, 1]. Then, to each node in the intersection tree, we attach children in T \ T and children in T \ T , each number being independent Poi(λ(1s)) variables. These children in turn have offspring in the corresponding tree given by independent Galton-Watson trees with offspring Poi(λ). See figure 3 for an illustration. We denote (T , T ) ∼ GW (λ, s, δ).

GW (λ, s): It is the previous model with δ = 0, so that the two correlated trees T and T have same root ρ. We denote (T , T ) ∼ GW (λ, s). We now turn to the analysis of matching rates for these models.

Proposition 1.1 Let λ > 1 and s ∈ [0, 1] such that λs > 1. For (T , T ) ∼ GW (λ, s), letting γ(λ, s) := γ(T , T ), we have: γ(λ, s) ≥ λs.
Proof Let T ∩ be the intersection tree between T and T . Branching process theory implies that (λs) -d L d (T ∩ ) converges almost surely to a random variable Z as d → ∞, such that P (Z > 0) = 1p, with p the extinction probability of the branching tree T ∩ . Since p < 1 when λs > 1, and for every small enough ε > 0, lim d→∞ P W d (T , T ) ≥ (λs(1ε)) d ≥ 1p , the result follows.

1.5. Matching rate of independent Galton-Watson trees Theorem 1.1 Let T , T be two independent Galton-Watson trees from the model GW (λ). Let γ(λ) := γ(T , T ). There exists λ 0 > 1 such that for all λ ∈ (1, λ 0 ], we have

γ(λ) < λ. (8)
Evaluations of γ(λ) by simulations, confirming the Theorem, are provided in Appendix B.1.

Proof outline

The full proof of Theorem 1.1 is detailed in the appendix (C.1), but we here give the key steps. We introduce some notations. First, for a tree t, let r d (t) denote the tree obtained by suppressing nodes at depths greater than d, and then iteratively pruning leaves of depth strictly less than d. When computing W d (t, t ), the only informative sub-trees are precisely r d (t) and in r d (t ), one of these being empty if t or t doesn't survive up to depth d. In the rest of the paper, we define T d the random variable r d (T ) where T is conditioned to survive up to depth d.

Consider (T , T ) ∼ GW (λ). We let E d (respectively, E d ) denote the event that tree T (respectively, T ) becomes extinct before d generations, i.e. L d (T ) = ∅ (respectively, L d (T ) = ∅). We let p d = P(E d ) = P(E d ). It is well known that it satisfies the recursion

p 0 = 0, p d = e -λ(1-p d-1 ) .
We now establish the following lemma on the structure of T d : Lemma 1.1 For any λ > 1, T d can be constructed by first sampling the number of children D of the root ρ(T ) according to distribution

P(D = k) = 1 k>0 P(Poi(λ(1 -p d-1 )) = k) P(Poi(λ(1 -p d-1 )) > 0) =: q d,k ,
and then attaching D independent copies of T d-1 to the D children of ρ(T ).

Assume ε = λ -1 to be small enough . Fix r ∈ (0, 1), let γ = 1 + rε. We first show using exponential moments that there exist m, c > 0 and d 0 > 0 such that for all x > 0

P W d 0 T d 0 , T d 0 ≥ mx ≤ e -x+c .
Then we define the random variables

X d := γ -(d-d 0 ) m -1 W d T d , T d .
Then, considering the number D of children of the root in T d (resp. D in T d ), using the previous lemma, one can establish, for all x > 0, a recursive formula of the following form

P (X d ≥ x) ≤ k, ≥1 q d,k q d, P   ∃m ∈ M ([k], [ ]) , (i,u)∈m X d-1,i,u ≥ γx   ,
where the X d-1,i,u are i.i.d. copies of X d-1 . The union bound yields

P (X d ≥ x) ≤ k, ≥1 q d,k q d, min 1, (k ∨ ) k∧ × P k∧ i=1 X d-1,i,u ≥ γx ,
where

m p := m(m -1) . . . (m -p + 1) = m! (m-p)!
. This inequality enables, with a few more technical steps (see C.1), to propagate recursively the inequality

P (X d ≥ x) ≤ e -(x-c) + .

IMPLICATIONS FOR A HYPOTHESIS TESTING PROBLEM

Let a pair of trees (T , T ) be distributed according to GW (λ) under the null hypothesis H 0 , and according to GW (λ, s) under the alternative hypothesis H 1 . They are thus independent under H 0 , and correlated under H 1 . Consider the following test:

Decide H 0 if W d (T , T ) < γ d , H 1 otherwise.
Assume that γ(λ) < γ < λs. Then in view of Remark 1.3 and Theorem 1.1 one has for some c(γ) > 1:

P decide H 1 H 0 = O e -c(γ) d ,
thus a super-exponential decay of the probability of false positive (first type error). Conversely, in view of Proposition 1.1, noting T ∩ the intersection tree under H 1 , one has

P decide H 0 H 1 , non-extinction of T ∩ = o d (1).
The false negative probability of this test thus also goes to zero, provided the intersection tree survives. As we will see in the next section, this hypothesis testing problem on a pair of random trees is related to our original graph alignment problem much as the so-called tree reconstruction problem, reviewed in [START_REF] Mossel | Survey: Information flow on trees[END_REF], is related to community detection in sparse random graphs (see e.g. [START_REF] Bordenave | Non-backtracking spectrum of random graphs: Community detection and non-regular ramanujan graphs[END_REF]).

Matching rate of intersecting trees

Theorem 1.2 Let (T , T ) ∼ GW (λ, s, δ) with δ ≥ 1 and s ∈ [0, 1]. Let γ(λ, s, δ) := γ(T , T ).
There exists λ 0 > 1 such that for all λ ∈ (1, λ 0 ] we have

sup δ≥1 γ(λ, s, δ) < λ. (9)
Evaluations of γ(λ, s, δ) by simulations, confirming the Theorem, are provided in Appendix B.1.

Proof outline

The full proof of Theorem 1.2 is detailed in the appendix (C.2), but we here give the key steps. The proof will again be by induction on d, the initial step being established with the same argument as in the proof of Theorem 1.1. ε = λ -1 is assumed to be small enough. We fix r ∈ (0, 1), and we let γ = 1 + rε . We now work with the random variables

X d := γ -(d-d 0 ) m -1 W d T d , T d ,
conditionally on the event that the path from ρ to ρ survives down to depth d in T . Then, considering D the number of children of ρ in T d , D the number of children of ρ in T d that are in the intersection tree T d ∩ T d , and D the number of children of ρ in T d \ T d , we establish for all x > 0 a recursive formula of the following form

P X d ≥ x ≤ k, ≥1 P D + D = k, D = min 1, (k ∨ ) k∧ P X d-1 + k∧ -1 i=1 X d-1,i,u ≥ γx ,
where the X d-1,i,u are i.i.d. copies of X d-1 as defined in the proof of Theorem 1.1. Again, with a few more technical steps (see C.2), we are able to propagate recursively the inequality

P X d ≥ x ≤ e -(x-c) + .

Sparse graph alignment

We now describe our main algorithm and its theoretical guarantees. For simplicity we assume that the underlying permutation σ is the identity.

Neighborhood Tree Matching Algorithm (NTMA), main result

The main intuition for the NTMA algorithm is as follows. In order to distinguish matched pairs of nodes (i, u), we consider their neighborhoods at a certain depth d, that are close to Galton-Watson trees. In the case where the two vertices are actual matches, the largest common subtree measured in terms of children at depth (exactly) d is w.h.p. of size ≥ (λs) d . However, when the two nodes i and u are sufficiently distant, previous study of matching rates shows that the growth rate of largest common subtree will be < λs. The natural idea is thus to apply the test comparing W d (i, u) to γ d for some well-chosen γ to decide whether i is matched to u. But as the reader may have noticed, testing W d (i, u) > γ d is not enough, because two-hop neighbors would dramatically increase the number of incorrectly matched pairs, making the performance collapse. To fix this, we use the dangling trees trick: instead of just looking at their neighborhoods, we look for the downstream trees from two distinct neighbors j = j of i, and

v = v of u. The trick is now to compare both W d-1 (j ← i, v ← u) and W d-1 (j ← i, v ← u) to γ d-1
. This way, even if i = u and i and u are close by, the pairs of rooted trees that can be considered will lead to one of the four cases considered and illustrated on Figure 10, that are settled in the proof of Theorem 2.2.

Our algorithm is as follows, where matching tree weights W d-1 (j ← i, v ← u) are defined in (4):

Algorithm 1 Neighborhood Tree Matching Algorithm for sparse graph alignment Input: Two graphs G 1 and G 2 of size n, average degree λ, depth d, parameter γ.

Output: A set of pairs S ⊂ V (G 1 ) × V (G 2 ). S ← ∅ for (i, u) ∈ V (G 1 ) × V (G 2 ) do if B G 1 (i, d) and B G 2 (u, d) contain no cycle, and ∃j = j ∈ N G 1 (i), ∃v = v ∈ N G 2 (u) such that W d-1 (j ← i, v ← u) > γ d-1 and W d-1 (j ← i, v ← u) > γ d-1 then S ← S ∪ {(i, u)} end end return S Remark 2.1 For d = c log n , in view of Remark 1.2, with high probability the complexity of NTMA is O |V (G 1 )| |V (G 2 )| (log n) 2 n 2c log λ d 2 max + O |E(G 1 )| |E(G 2 )| (log n)d 3 max ,
where d max is the maximum degree in G 1 and G 2 . In the context of Theorems 2.1 and 2.2 the complexity is then O (log n) 4 n 5/2 .

The two results to follow will readily imply Theorem I.

Theorem 2.1 Let (G 1 , G 2 ) ∼ ERC(n, λ/n, s) be two s-correlated Erdős-Rényi graphs such that λs > 1. Let d = c log n with c log (λ (2s)) < 1/2. Then for γ ∈ (1, λs), with high probability,

1 n n i=1 1 {(i,i)∈S} = Ω(1). ( 10
)
In other words, a non vanishing fraction of nodes is correctly recovered by NTMA (1).

Theorem 2.2 Let (G 1 , G 2 ) ∼ ERC(n, λ/n, s) be two s-correlated Erdős-Rényi graphs. Assume that γ 0 (λ) := max γ(λ), sup δ≥1 γ(λ, s, δ) < λs, and that d = c log n with c log λ < 1/4. Then for γ ∈ (γ 0 (λ), λs), with high probability,

err(n) := 1 n n i=1 1 {∃u =i, (i,u)∈S} = o(1), (11)
i.e. only at most a vanishing fraction of nodes are incorrectly matched by NTMA (1).

Remark 2.2 The set S returned by the NTMA is not necessarily a matching. Let S be obtained by removing all pairs (i, u) of S such that i or u appears at least twice. Theorems 2.1 and 2.2 guarantee that S still contains a non-vanishing number of correct matches and a vanishing number of incorrect matches. Theorem I easily follows. Simulations of NTMA-2, a simple variant of of NTMA, are reported in Appendix B.2. These confirm our theory, as the algorithm returns many good matches and few mismatches.

Proof strategy

We start by stating Lemmas, adapted from [START_REF] Massoulié | Community detection thresholds and the weak Ramanujan property[END_REF] and [START_REF] Bordenave | Non-backtracking spectrum of random graphs: Community detection and non-regular ramanujan graphs[END_REF] and proven in Appendix D, that are instrumental in the proofs of Theorems 2.1 and 2.2.

Lemma 2.1 (Control of the sizes of the neighborhoods)

Let G ∼ ER(n, λ/n), d = c log n with c log λ < 1. For all γ > 0, there is a constant C = C(γ) > 0 such that with probability 1 -O (n -γ ), for all i ∈ [n], t ∈ [d]: |S G (i, t)| ≤ C(log n)λ t . ( 12
)
Lemma 2.2 (Cycles in the neighborhoods in an ER graph) Let G ∼ ER(n, λ/n), d = c log n with c log λ < 1/2. There exists ε > 0 such that for any vertex i ∈ [n], one has

P (B G (i, d) contains a cycle) = O n -ε . ( 13
)
Lemma 2.3 (Two logarithmic neighborhoods are typically size-independent) Let G ∼ ER(n, λ/n) with λ > 1, d = c log n with c log λ < 1/2. Then there exists ε > 0 such that for any fixed nodes i = j, the variation distance between the joint law of the neighborhoods L (S G (i, t), S G (j, t)) t≤d and the product law L (S G (i, t)) t≤d ⊗ L (S G (j, t)) t≤d tends to 0 as O (n -ε ) for some ε > 0 when n → ∞.

Lemma 2.4 (Coupling the |S G (i, t)| with a Galton-Watson process) Let G ∼ ER(n, λ/n), d = c log n with c log λ < 1/2. For a fixed i ∈ [n], the variation distance between the law of (|S G (i, t)|) t≤d and the law of (Z t ) t≤d where (Z t ) t is a Galton-Watson process of offspring distribution Poi(λ) tends to 0 as O (n -ε ) when n → ∞.

PROOF OF THEOREMS 2.1 AND 2.2

Proof [Proof of Theorem 2.1] Define the joint graph G ∪ = G 1 ∪ G 2 . For i ∈ [n], let M i denote the event that the algorithm matches i in G 1 with i in G 2 , i.e. on which B G 1 (i, d) and B G 2 (i, d) contain no cycle, and ∃j

= j ∈ N G 1 (i), ∃v = v ∈ N G 2 (i) such that W d-1 ((j, v), (i, u)) > γ d-1
and W d-1 ((j , v ), (i, u)) > γ d-1 . Denote by C ∪,i,d the event that there is no cycle in B G∪ (i, d).

Arguing as in the proof of Lemma 2.4, the two neighborhoods B G 1 (i, d) and B G 2 (i, d) can be coupled with trees distributed as GW (λ, s) of Section 1. However, we will instead consider the intersection graph G ∩ = G 1 ∩ G 2 . Obviously, G ∩ ∼ ER(n, λs/n). By Lemma 2.4, the random variables |S G∩ (i, t)| can be coupled with a Galton-Watson process with offspring distribution Poi(λs) up to depth t = d. Let P i denote the event that this coupling succeeds. Since λs > 1, there is a probability 2α > 0 that this process survives up to depth d and that the first generation has at least two children. Note S this event. On event S, the matching given by the identity on the intersection tree implies the existence of two neighbors j

= j ∈ N G 1 (i) and v = v ∈ N G 2 (i) such that with high probability W d-1 (j ← i, v ← u) > γ d-1 and W d-1 (j ← i, v ← u) > γ d-1
, by standard martingale arguments, as in Proposition 1.1. This gives the lower bound for P(M i ):

P(M i ) ≥ P (C ∪,i,d ∩ P i ∩ S) ≥ 2α -o(1) > α > 0.
It is easy to see that G ∪ ∼ ER(n, λ(2s)/n). For i = j ∈ [n], define I i,j the event on which the two neighborhoods of i and j in G ∪ coincide with their independent couplings up to depth d. By lemma 2.3, P(I i,j ) = 1o(1). Then for 0 < ε < α Markov's inequality yields

P 1 n n i=1 1 {(i,i)∈S} < α -ε ≤ P n i=1 (P(M i ) -1 M i ) > εn (14) ≤ 1 n 2 ε 2 (nVar (1 M 1 ) + n(n -1)Cov (1 M 1 , 1 M 2 )) (15) ≤ Var (1 M 1 ) nε 2 + 1 -P (I 1,2 ) ε 2 → 0. ( 16
)
Proof strategy for Theorem 2.2. Consider two distinct nodes i and u. We place ourselves on the event of high probability that B G∪ (i, 2d) has no cycle. On this event, the two neighborhoods B G 1 (i, d) and B G 2 (u, d) can be coupled with two trees rooted at i, u respectively. We then distinguish several cases that are shown on Figure 10, that require detailed analysis, and which all show that for i fixed, one has

P (∃u = i, (i, u) ∈ S) = o(1).
The full proof is deferred to Appendix D.5.

Conclusion

We have introduced NTMA, an algorithm we proved to succeed at partial alignment of sparse correlated random graphs. While our Theorem applies to a limited range of average degrees λ, we conjecture that rates γ(λ) and γ(λ, s, δ) are strictly less than λ for all λ > 1 and s < 1 and thus NTMA in fact succeeds for a much broader parameter range. This will be the object of future work. and converges monotonically to the smallest root in [0, 1] of x = e -λ (1-x) . This root, that we denote p e , is the probability of ultimate extinction. For small enough ε = λ -1, it holds that

p e = 1 -2ε + O(ε 2 ),
as can be seen by analysis of the fixed point equation satisfied by p e . Let then d 0 be such that for all d ≥ d 0 , p d = 1 -2ε + O(ε 2 ). Clearly, on the event E d ∪ E d , the set of matchings M d (T , T ) is empty, so that W d (T , T ) = 0. Recall that we define T d the random variable r d (T ) where T is conditioned to survive up to depth d. Now fix r ∈ (0, 1). We shall prove that for sufficiently small ε > 0, letting γ = 1 + rε, there exists some constants c, m, d 0 > 0 such that for all x > 0, all d ≥ d 0 , one has

P W d (T d , T d ) ≥ γ d-d 0 mx ≤ e -(x-c) + . ( 17
)
We proceed by induction over dd 0 . To initialize the induction, notice that one obviously has

W d 0 (T d 0 , T d 0 ) ≤ |L d 0 (T d 0 )| =: X. By Corollary C.
1, for all m, x, θ > 0, one has:

P W d 0 (T d 0 , T d 0 ) > mx ≤ P(X > mx) ≤ Ee θX e -θmx .
Let now θ = 1/m. By taking m sufficiently large, from dominated convergence we can make Ee (1/m)X as close to 1 as we like. Choose for instance m such that Ee (1/m)X ≤ 2. Then

P(W d 0 (T d 0 , T d 0 ) > mx) ≤ 2e -x ≤ e -x+c .
for any c ≥ ln(2). Hence, for sufficiently large m, we can initialize the induction at d = d 0 with any c ≥ ln(2).

Recall we set γ = 1 + rε. Define the random variables

X d := γ -(d-d 0 ) m -1 W d T d , T d .
Let

W d (T d , T d ) = sup m∈M([D],[D ]) (i,u)∈m W d-1 (T d-1,i , T d-1,u ), where M([D], [D ]) denotes the set of all (D ∨ D ) D∧D maximal injective mappings between E 0 ⊆ [D] and [D ]. Let X d-1,i,u := γ -(d-1-d 0 ) m -1 W d-1 (T d-1,i , T d-1,u ).
Note that conditional on D and D , for each matching m ∈ M([D], [D ]), the variables (X d-1,i,u ) (i,u)∈m are i.i.d. with the same distribution as X d-1 . The induction hypothesis states that each X d-1,i,u is less, for the strong stochastic ordering of comparison of cumulative distribution functions, than c plus an exponential random variable with parameter 1. With an easy union bound, we can derive the following bounds:

P (X d > x) ≤ 1≤k≤ <∞ P D ∧ D = k, D ∨ D = min 1, k P (E 1 + . . . + E k > γx -kc) ,
(18) where E 1 , . . . , E k are independent exponential random variables of parameter 1. Lemma 1.1 states that

P(D = k) = e -λ(1-p d-1 ) λ k (1 -p d-1 ) k k! (1 -p d ) =: q d,k .
We can increase d 0 such that for some constant κ > 0, for all d ≥ d 0 :

q d,1 ≤ 1 -ε + κε 2 , q d,k ≤ (3ε) k-1 k! , k ≥ 2.
Note that for x ≤ c, there is nothing to prove in ( 17), since a probability is always upper-bounded by 1. We thus only need to consider the case x > c. We conclude the proof of this Theorem by appealing to the following Lemma C.1 Let κ, C > 0 and r ∈ (0, 1) be given constants. Then there exists c > 0 large enough and ε 0 > 0 such that, for all ε ∈ (0, ε 0 ),

letting γ = 1 + rε, q 1 = 1 -ε + κε 2 , q k = (Cε) k-1 /k! for k ≥ 2, one has ∀x > c, k, ≥1 q k q l min 1, (k ∨ ) k∧ P (E 1 + . . . + E k∧ > γx -(k ∧ )c) ≤ e -(x-c) , (19)
where the E i are independent exponential random variables of parameter 1.

Its assumptions are indeed verified here with C = 3, so (17) can be propagated by using this Lemma in (18), and the conclusion of Theorem 1.1 follows.

C.2. Proof of Theorem 1.2

Proof [Proof of Theorem 1.2] We assume that λ = 1 + ε. We fix r ∈ (0, 1), and we let γ = 1 + rε for some fixed r ∈ (0, 1). . We work with trees such that (T , T ) ∼ GW (λ, s, δ). If we assume that the path from ρ to ρ does not survive down to depth d in T , then this path is no more present in T d , and the two trees T d and T d can be coupled with two trees T d and T d where ( T , T ) ∼ GW (λ), and we are in the case of Theorem 1.1.

In the following proof, we will thus condition to the event S ρ,d that the path from ρ to ρ survives down to depth d in T . Recall that the tree T d (resp. T d ) is obtained, conditionally on the fact that T (resp. in T ) survives down to depth d, by suppressing nodes at depth greater than d in T (resp. in T ), and then pruning alternatively leaves of depth strictly less than d. As in the proof of Theorem 1.1, we shall establish that for sufficiently small ε > 0, there exist constants c, m, d 0 > 0 such that for all x > 0, all d ≥ d 0 , one has

P W d (T d , T d ) ≥ γ d-d 0 mx S ρ,d ≤ e -(x-c) + . (20)
Define the random variables Recall that p d denotes the probability that a Galton-Watson tree with offspring Poi(λ) becomes extinct before d generations. Then, conditionally on S ρ,d , the random variables D, D and D have the following distributions:

X d := γ -(d-d 0 ) m -1 W d T d+δ , T d , conditional on S ρ
D ∼ 1 + Poi (λ (1 -p d-1 )) , D ∼ Poi (λs (1 -p d-1 )) D ∼ Poi (λ(1 -s) (1 -p d-1 )) , conditional on D + D > 0.
We show an illustration on Figure 9. We condition on the values , k , k taken by D, D , D . The number of maximal one-to-one mappings between the children of ρ in T d and those of ρ in T d is given by [(k + k ) ∨ ] (k +k )∧ , and each of them is of size ∧ (k + k ). Note here again that for a fixed matching between the children of ρ and ρ , the weights of the matched sub-trees are independent. We distinguish between several cases (to help understand these cases, the reader could keep figure 9 in mind):

• For a child i of ρ that is not on the path to ρ , the corresponding sub-trees are independent so that the corresponding weight is distributed as W d-1 ( T d-1 , T d-1 ) in the independent-tree model GW (λ).

• If the child of ρ on the path to ρ is matched with a child of ρ that is not in the intersection tree, again the corresponding weight is similarly distributed.

• Finally, if the child i of ρ leading to ρ is matched to a child u of ρ in the intersection tree, setting the new root at ρ := i in T and at ρ := u in T , the corresponding weight has the same distribution as

W d-1 (T d-1 , T d-1
) in the model GW (λ, s, δ), still conditioned to S ρ,d-1 . Indeed, there is a path from ρ to ρ , and the corresponding Poisson distributions are conserved.

The induction hypothesis for case 3, together with Theorem 1.1 for cases 1 and 2, therefore give us:

P X d ≥ x ≤ k, P D + D = k, D = min 1, (k ∨ ) k∧ P (E 1 + . . . + E k∧ > γx -(k ∧ )c) ,
where the E i are independent exponential random variables of parameter 1. Assume, as in the proof of Theorem 1.1, that d 0 is chosen such that for all d ≥ d 0 ,

p λ d = 1 -2ε + O(ε 2 ).
With simple computations, we can then ensure that for some κ > 0, noting q d,• the distribution of D, one has

q d,1 ≤ 1 -ε + κε 2 , q d,k ≤ (3ε) k-1 (k -1)! ≤ (6ε) k-1 k! , k ≥ 2,
where we used k ≤ 2 k-1 in the last step. By independence of D and D , D + D follows a Poi(λ(1p d-1 )) distribution, conditional on being positive. Noting q d,• this distribution, we have, as in the previous proof,

q d,1 ≤ 1 -ε + κε 2 , q d,k ≤ (3ε) k-1 k! , k ≥ 2.
We can then invoke Lemma C.1 to conclude. Note that every control in the proof is made uniformly on δ ≥ 1 .

C.3. Proof of lemma C.1 Proof [Proof of Lemma C.1]. Let
S 1 := e x-c q 2 1 e -(γx-c) + + 4q 1 q 2 e -(γx-c) + , S 2 := 2e x-c q 1 ≥3 q min 1, e -(γx-c) + , S 3 := 2e x-c 2≤k≤ q k q min 1, k P (E 1 + . . . + E k > γxkc) .

Our goal is to show that for a suitable choice of c, for all x > c, S 1 + S 2 + S 3 ≤ 1. One has

S 1 ≤ e -rεx (1 -ε + κε 2 ) 2 + 2Cε ≤ e -rεx (1 + 2Cε), (21)
and

S 2 ≤ 2e -rεx (1 -ε + κε 2 ) ≥3 (Cε) -1 ( -1)! ≤ 2e -rεx e Cε -1 -(Cε) ≤ 2e -rεx C 2 ε 2 . ( 22
)
We let k 0 be such that γx ∈ [k 0 c, (k 0 + 1)c). We then upper-bound S 3 by A + B where

A = 2e x-c k 0 k=2 q k ≥k q ! ( -k)! P (E 1 + . . . + E k > γx -kc) , (23)
B = 2e x-c k≥(k 0 +1)∨2 ≥k q k q . ( 24
)
One readily has

B ≤ 2e -rεx e γx-c k≥(k 0 +1)∨2 (Cε) k-1 k! ≥k (Cε) -1 ! (25) ≤ 2e -rεx e γx-c k≥(k 0 +1)∨2 (Cε) 2(k-1) k! (26) ≤ 2e -rεx e k 0 c (Cε) 2((k 0 +1)∨2-1) (27) ≤ 2e -rεx (Cεe c ) 2 , ( 28
)
where in the last steps we assumed that Cεe c < 1, so that

e k 0 c (Cε) 2((k 0 +1)∨2-1) ≤ (Cεe c ) 2((k 0 +1)∨2-1) ≤ (Cεe c ) 2 .
Note that for y ≥ 0, P (E 1 + . . . + E k > y) = P(Poi(y) < k) = e -y k-1 j=0 y j /j!. Write then

A ≤ 2e x-c k 0 k=2 (Cε) 2(k-1) k! k-1 j=0 e -γx+kc (γx) j j! ≤ 2e -rεx k 0 k=2 (C 2 e c ) k k! k-1 j=0 (γxε 2 ) j j! ε 2(k-1-j) ≤ 2e -rεx k 0 k=2 (C 2 e c ) k k! ε 2(k-1) + e γxε 2 -1 ≤ 2e -rεx ε -2 e ε 2 C 2 e c -1 -ε 2 C 2 e c + e C 2 e c e γxε 2 -1 (29)
Summing the upper bounds ( 21)-( 29), the desired property will then hold if for all x > c, one has:

e -rεx 1 + 2Cε + 2C 2 ε 2 + 2[Cεe c ] 2 + 2ε -2 e ε 2 C 2 e c -1 -ε 2 C 2 e c + 2e -rεx e C 2 e c e γxε 2 -1 ≤ 1. (30)
The first term is, for any fixed c, and for sufficiently small ε, upper bounded e -rεx (1 + (2C + 1)ε).

We now distinguish three cases for x.

• Case 1:

x ∈ [c, 1/ √ ε]. The second term is then O(ε √ ε). Provided rc > 2C + 1, since e -rεx ≤ e -rεc = 1 -rεc + O(ε 2
), the left-hand side of ( 30) is then upper-bounded by 1 -(rc -2C -1)ε + O(ε √ ε), and is thus less than 1.

• Case 2: x ∈ [1/ √ ε, 1/ε]. Since e -rεx ≤ e -r √ ε = 1 -Ω( √ ε)
, and e γxε 2 -1 ≤ e γε -1 = O(ε), the left-hand side of ( 30) is upper-bounded by 1 -Ω( √ ε) and is thus less than 1.

• Case 3: x ≥ 1/ε. The first term is then bounded by e -r (1 + (2C + 1)ε), which is less than 1 -Ω(1) for ε small enough. Letting y = εx, the second term reads e -ry [e εγy -1](2e C 2 e c ).

For small ε, this function is maximized for y = 1/r + O(ε), at which point it evaluates to O(ε). Thus the left-hand side of ( 30) is upper-bounded by 1 -Ω(1) in that range.

We have thus shown that for any r > 0, provided c > (2C + 1)/r, then for all sufficiently small ε, the desired property holds with γ = 1 + rε.

C.4. Proof of lemma 1.1

Proof [Proof of Lemma 1.1] For a tree t, we write t = (t 1 , . . . , t k ) to represent the fact that its root has k children, whose offsprings are given by trees t 1 , . . . , t k . Write, noting by D the number of children of ρ(T ), fixing k ≥ 1, t 1 , . . . , t k ∈ A d-1 , and letting S = i 1 < • • • < i k run over all k subsets of [ ]:

P(T d = (t 1 , . . . , t k )) = ≥0 P(T d = (t 1 , . . . , t k ), D = ) = ≥k S P D = , r d (T i j ) = t j , j ∈ [k], r d (T v ) = ∅, v / ∈ S E d = 1 1 -p d ≥k k e -λ λ ! p -k d-1 k j=1 P(T d-1 = t j )(1 -p d-1 ) = 1 1 -p d (λ(1 -p d-1 )) k k! k j=1 P(T d-1 = t j ) ≥k e -λ (λp d-1 ) -k ( -k)! = e -λ(1-p d-1 ) 1 -p d (λ(1 -p d-1 )) k k! k j=1 P(T d-1 = t j )
The conclusion follows by noting that 1p d = 1e -λ(1-p d-1 ) .

D.1. Proof of Lemma 2.1

Proof [Proof of Lemma 2.1] Fix K > 0 to be specified later and γ > 0. Fix i ∈ [n], and define

T := inf {t ≤ d, |S G (i, t)| ≥ K log n} . If T = ∞, there is nothing to prove. Given |S G (i, T -1)|, |S G (i, T )| ∼ Bin n -|S G (i, 0)| -. . . -|S G (i, T -1)| , 1 -1 - λ n |S G (i,T -1)| . Thus |S G (i, T )| sto.
≤ Bin n, λK log n n .

Using Bennett's inequality, for K > λK:

P |S G (i, T )| ≥ K log n ≤ e -λKh K -λK λK log n , with h(u) = (1 + u) log(1 + u) -u. This probability is ≤ n -2-γ if K is large enough to ver- ify λKh K -λK λK > γ + 2. With a simple use of the union bound, one gets that |S G (i, T )| ∈ [K log n, K log n] for all i ∈ [n] with probability 1 -O(n -1-γ ).
Take ε > 0 to be specified later. We then check by induction that with high probability, for all

T ≤ t ≤ d, |S G (i, t)| ∈ K λ 2 t-T (log n) t s=T 1 -ε λ 2 -s-T 2 , K λ t-T (log n) t s=T 1 + ελ -s-T 2 .
(31) The case t = T is proved here above. We will next use the inequality

λu/(2n) ≤ λu/n -λ 2 u 2 /(2n 2 ) ≤ 1 -(1 -λ/n) u ≤ λu/n. (32)
that holds as soon as λu/n < 1.

Assuming (31) holds up to t, inequality (32

) holds for u = |S G (i, t)| for n large enough, since |S G (i, t)| < n/λ for c log λ < 1. Thus for n large enough E |S G (i, t + 1)| lies in the interval        K 2 λ λ 2 t-T (log n) t s=T 1 -ε λ 2 -s-T 2 =1-O(ε) , λK λ t-T (log n) t s=T 1 + ελ -s-T 2       
With ε > 0 to be specified later, Bennett's inequality writes This coupling is well defined, and coincides with the independent setting up to step k as long as the sets s≤k S(i, s) and s≤k S(j, s) do not intersect. On the event

P |S G (i, t + 1)| -E |S G (i, t + 1)| ≥ εE |S G (i, t + 1)| ≤ 2e -1 2 λ( λ 2 ) t-T log n(1-O(ε))h(ε) , • Now if k > 1, construct S(i, k) as follows: select a subset of [n] \ s≤k-1 S(i,
A := t≤d |S(i, t)| , |S(j, t)| < C(log n)λ t , one has E   k≤d S(i, s) ∩ k≤d S(j, s)   ≤ E d k=1 Bin C(log n)λ k , k t=1 C(log n)λ t n -k t=1 C(log n)λ t ≤ C 2 (log n) 2 λ λ -1 d k=1 λ 2k n -λ λ-1 C(log n)λ k ≤ O (log n) 2 λ 2d /n if (log n)λ d = o(n), which is the case if c log λ < 1. The expectation is upper-bounded by O (log n) 2 λ 2d /n = O (log n) 2 n -2ε if c log λ ≤ 1/2 -ε. With Lemma 2.1, choosing C such that P (A) = 1 -O n -2ε , we get d TV L (S G (i, t), S G (j, t)) t≤d , L (S G (i, t)) t≤d ⊗ L (S G (j, t)) t≤d ≤ O((log n) 2 n -2ε ) + P Ā ≤ O(n -ε ). D.4. Proof of Lemma 2.4
We work here conditionally on

A := t≤d |S G (i, t)| < C(log n)λ t .
Let's define a Galton-Watson process as follows: set Z 0 = 1, and for t > 0, L (Z .

The Stein-Chen method (see e.g. [START_REF] Barbour | An Introduction to Stein's Method. co-published with Singapore University[END_REF]) enables to bound d TV (Bin(n, λ/n), Poi(λ)) by min(1, λ -1 )λ 2 /n ≤ λ/n. We also use the classical bound d TV (Poi(λ), Poi(λ )) ≤ |λλ | together with inequality (32) (which holds for n large enough since c log λ < 1) to obtain that conditionally on F t-1 : We use the same notations as in the former proof:

G ∪ = G 1 ∪ G 2 and G ∩ = G 1 ∩ G 2 . Fix i ∈ [n].
In the rest of the proof we work conditionally to the event C ∪,i,2d that B G∪ (i, 2d) has no cycle. Since c log λ < 1/4, P (C ∪,i,2d ) = 1o(1) by Lemma 2.2. Fix another vertex u = i. The d-neighborhoods B G 1 (i, d) and B G 2 (u, d) have offspring distribution stochastically dominated by Bin(n, λ/n), which is also dominated by Poi(λ ) as soon as λ = λ + O(1/n) (see e.g. [START_REF] Klenke | Stochastic ordering of classical discrete distributions[END_REF]). We can choose λ such that γ > γ(λ , 0) still holds: indeed, by a standard coupling argument, one can see that γ : λ → γ(λ) is increasing. We now build two dominating (in the usual edge presence sense) tree-like d-neighborhoods of i and u with the following construction.

• First, if the two neighborhoods don't intersect, just sample two independent trees from model GW (λ ) rooted in i and in u.

• If the two neighborhoods intersect, condition to the event that α is the contact point in the path p ∪ (unique by conditioning on C ∪,i,2d ) from i to u in the joint graph. Then there is a path of edges of G 1 (say, blue) from i to α, then a path of edges of G 2 (say, red) from α Now assume that p ∩ exists, and let v = v two neighbors of u and j = j two neighbors of i. Case (iii.a): At least one of the edges (i, j), (i, j ), (u, v), (u, v ) is not in G ∩ (e.g. edge (i, j) on figure 10, bottom left): again, the same argument applies. Case (iii.b): Edges (i, j), (i, j ), (u, v), (u, v ) are all in G ∩ (see figure 10, bottom right). Then one pair of vertices (say (j , v ) as on figure 10) can be on p ∩ and bring a high W d-1 (j ← i, v ← u) > γ d-1 matching weight, if their descendants spread over a great part of the intersection. In that case, since j and v can't be on p ∩ , the associated downstream trees are independent, and again W d-1 (j ← i, v ← u) < γ d-1 with high probability. The remaining case to be considered is that of matches (j, v ) and (j , v), with j , v on p ∩ . All trees involved are then correlated. However, the coupling construction induces a coupling of the two pairs of (d -1)-neighborhoods (from (j, v ) and from (j , v), see figure 10) with two pairs of trees from model GW (λ , s, δ) where δ = |p ∩ |. The Theorem assumes γ(λ, s, δ) < λs so that, by Theorem 1.2, the probability that W d-1 (j ← i, v ← u) > γ d-1 and W d-1 (j ← i, v ← u) > γ d-1 is upper bounded by O (exp (-n ε )).

Thus, for i fixed, one has P (∃u = i, (i, u) ∈ S) ≤ 1 -P (C ∪,i,2d ) + n × P (C ∪,i,2d ) × d 2 max × O (exp (-n ε )) = o(1).

The Theorem then follows by appealing to Markov's inequality.

Figure 2 :

 2 Figure 2: We here show an example of a tree T and its corresponding T (i←j) 2 highlighted in blue.

Remark 1. 2

 2 The complexity of computing all matching weights at depth d can be obtained as follows. We use dynamic programming and store the W k (i ← j, u ← v) in a array of size the number of pairs (e, e ) where e and e are two oriented edges in T ,T (that is, 4 × |T | × |T |). Each time we increase k, we solve one LAP for each pair (e, e ). The size of the small matrix on which the LAP is done does not exceed d max , the maximal degree in T and T . The Hungarian algorithm solves LAP with cubic time complexity. The time complexity is thus O d × |T | × |T | × d 3 max . As d max = O(log n) with high probability, the total complexity is O(dn 2 log 3 n) where n bounds the number of nodes in T and T .

ρ ρ 1 Figure 3 :

 13 Figure 3: Random trees T (blue) and T (red) from model GW (λ, s, δ) with δ = 3.

  (a) λ = 1.2, log λ ∼ 0.18. Red dashed slope ∼ 0.12. (b) λ = 2.2, log λ ∼ 0.79. Red dashed slope ∼ 0.65.

  (c) λ = 3.2, log λ ∼ 1.16. Red dashed slope ∼ 1.03.

Figure 4 :

 4 Figure 4: Comparison of d log λ (blue) and log W d (T , T ) (red) for W d (T , T ) ∼ GW (λ) conditioned to survive (100 iterations).

  (a) λ = 2.1, log λ ∼ 0.74, δ = 1. Red dashed slope ∼ 0.63. (b) λ = 2.1, log λ ∼ 0.74, δ = 2. Red dashed slope ∼ 0.63.(c) λ = 2.1, log λ ∼ 0.74, δ = 5. Red dashed slope ∼ 0.62.

Figure 5 :

 5 Figure 5: Comparison of d log λ (blue) and log W d (T , T ) (red) for W d (T , T ) ∼ GW (λ, s, δ) with s = 1, conditioned to survive (50 iterations).

  (a) s = 0.95. (b) Isomorphism case , s = 1.0.

Figure 6 :

 6 Figure 6: Mean score of NTMA-2 for λ = 2.1, d = 5 (25 iterations per value of n).

  (a) s = 0.95. (b) Isomorphism case , s = 1.0.

Figure 7 :

 7 Figure 7: Mean score of NTMA-2 for λ = 3.1, d = 4 (25 iterations per value of n).

  (a) n = 150, λ = 1.4, d = 5. (b) n = 50, λ = 2.2, d = 3.

Figure 8 :

 8 Figure 8: Mean score of NTMA-2 with different values of s (25 iterations per value of n).

  D (resp. D ) denote the number of children of the root in T d (resp. D in T d ). Given D and D , noting T d = (T d-1,1 , . . . , T d-1,D ) and T d = T d-1,1 , . . . , T d-1,D , we have that

  ,d . The proof will again be by induction on d, the initial step being established with the same argument as in the proof of Theorem 1.1. Note that this argument does not depend on δ. Denote by D the number of children of ρ in T d , D the number of children of ρ in T d that are in the intersection tree T d ∩ T d , and D the number of children of ρ in T d \ T d . By branching property, note that these three variables are independent.

Figure 9 :

 9 Figure 9: Random trees T (blue) and T (red) from Figure 3 (left), and the results T d and T d after applying r d (right). In this example, δ = 3 and d = 6, D = 2, D = 2 and D = 1.

  s) of size S(i, k) uniformly at random. Then we construct independently S(j, k) taking a uniform subset of [n] \ s≤k-1 S(j, s) of size S(j, k) .

d

  TV (|S G (i, t)| , W t ) ≤ n -1 (n -|S G (i, 0)| -. . . -|S G (i, t -1)|) λ |S G (i, t -1)| n + (n -|S G (i, 0)| -. . . -|S G (i, t -1)|) 1 -1 -λ n |S G (i,t-1)| λ |S G (i, t -1)| ≤ λ |S G (i, t -1)| n + λ |S G (i, t -1)| -(n -|S G (i, 0)| -. . . -|S G (i, t -1)|) λ |S G (i, t -1)| n + λ 2 |S G (i, t -1)| 2 2n . Now, for ε > 0 such that c log λ ≤ 1/2ε, on the event A, all variables |S G (i, s)| are bounded by C(log n)n 1/2-ε . This leads to d TV (|S G (i, t)| , W t ) ≤ O (log n)n -1/2-ε + O (log n) 3 n -2ε + O (log n) 2 n -2ε = O (log n) 3 n -2ε .This proves by induction that the total variation distance between(|S G (i, t)|) t≤d and (Z t ) t≤d is bounded by O (log n) 4 n -2ε = O (n -ε ), taking C large enough in Lemma 2.1 so that P (A) ≥ 1 -O n -2ε . D.5. Proof of Theorem 2.2 Proof [Proof of Theorem 2.2] Define d max := max max i deg G 1 (i), max u deg G 2 (u) .

  t |G t-1 ) = Poi (λZ t-1 ), where G t = σ (Z s , s ≤ t). Fix t > 0. Conditionally on F t-1 := σ (|S G (i, s)| , s ≤ t -1), define a random variable W t with distribution Poi (λ |S G (i, t -1)|). Note that

L |S G (i, t)| F t-1 = Bin n -|S G (i, 0)| -. . . -|S G (i, t -1)| , 1 -1 -λ n |S G (i,t-1)|

Acknowledgments

This work was partially supported by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute).

Appendix A. Algorithms for matching weights

We here describe algorithms to compute recursively matching weights W d (i ← j, u ← v) and W d (i, u).

Simulations for tree matching

We here present some simulations of matching rates γ(λ) (figure 4) and γ(λ, s, δ) for s = 1 (figure 5) in order to illustrate Theorems 1.1 and 1.2 and the final conjecture. For these simulations, error bars correspond to one standard deviation.

B.2. Simulations for a simple variant algorithm of NTMA

We here present some simulations of simple variant algorithm of NTMA, NTMA-2, which happens to be more efficient in practice. The Algorithm NTMA-2 is as follows.

Algorithm 4 NTMA-2 Input: Two graphs G 1 and G 2 of size n, average degree λ, depth d, parameter γ.

This algorithm only selects rows and columns weight maximums and match the corresponding pairs. The last part ensures that S is a matching. For these simulations, error bars correspond to a confidence interval for the mean value of scores. In figures 6 and 7 we compare the scores of NTMA-2 for s = 0.95 with the isomorphism case s = 1.0, for different values of n. We illustrate the fact that nearly no vertex is mismatched, whereas a non-negligible fraction of nodes is indeed recovered. In figure 8, we compare the scores of NTMA-2 for fixed n but varying s, illustrating the existence of a 'critical' parameter s * (λ).

Proof This is easily seen by induction, based on the structure of T d given in Lemma 1.1.

Recall that we let E d (respectively, E d ) denote the event that tree T (respectively, T ) becomes extinct before d generations, i.e. L d (T ) = ∅ (respectively, L d (T ) = ∅). We let p d = P(E d ). It is well known that it satisfies the recursion

Appendix D. Detailed proofs for Section 2

The following proofs are adapted from the previous work of [START_REF] Massoulié | Community detection thresholds and the weak Ramanujan property[END_REF] and [START_REF] Bordenave | Non-backtracking spectrum of random graphs: Community detection and non-regular ramanujan graphs[END_REF].

with ε small enough and K large enough such that Kε > 2 + γ. Thus (31) holds for t + 1 with probability 1 -O(n -2-γ ).

All this ensures that the desired inequality ( 12) holds for all i ∈

Note that k * ≥ 2, and that if k * = ∞ then B G (i, d) does not contain any cycle. Now assume that k * < ∞. For any k ≥ 2, k * = k if and only if there are two vertices of S G (i, k -1) that are connected, or if there is a vertex of S G (i, k) connected to two vertices of S G (i, k -1). On the event

this happens with probability at most

Taking ε > 0 such that c log λ ≤ 1/2ε, choosing C such that P (A) = 1 -O n -2ε with Lemma 2.1, the probability that B G (i, d) contains a cycle is less than

t≤d and S(j, t) t≤d denote two independent realizations of the neighborhoods (i.e. with independent underlying Bernoulli variables). We then construct recursively a coupling (S(i, t), S(j, t)) t≤k :

• For k = 1, take S(i, t) to be a set of vertices uniformly chosen among sets of [n] of size S(i, 0) . Independently, take S(j, t) to be a set of vertices uniformly chosen among sets of

[n] of size S(j, 0) . to u. Next, complete this construction: along p ∪ , propagate the blue path from α towards u with probability s on each edge, stopping at the first time when one red edge is not selected. Do the symmetrical construction to propagate the red path from α towards i. Finally, to each double-colored vertex, attach independent realizations of model GW (λ , s), and to each single-colored vertex, attach independent realizations of model GW (λ ).

Note that these constructions lead to at most one path p ∪ between i and

Denote by p ∩ this hypothetical path (cf. figure 10). We then distinguish between several cases. Case

The construction gives a coupling with two independent trees from model GW (λ). By assumption γ(λ) < λs, the probability that there exist j in N G 1 (i) and v in N G 2 (u) such that W d-1 (j ← i, v ← u) > γ d-1 is upper bounded by O d 2 max exp (-n ε ) , following Remark 1.3. Hence i is matched to u with at most this probability.

Case (ii): δ G∪ (i, u) ≤ 2d but p ∩ does not exist (see figure 10, top right). Take v = v two neighbors of u and j = j two neighbors of i. Then (at least) one of these vertices is not on p ∪ (e.g. vertex j on figure 10): the downstream tree from this vertex is independent from every other neighborhood in the other graph. They can be coupled with model GW (λ), and the same bound as in case (i) holds.


Baudouin Denis De Senneville

José V Manj Ón

Pierrick Coup
email: pierrick.coupe@u-bordeaux.fr

RegQCNET: Deep Quality Control for Image-to-template Brain MRI Affine Registration

Keywords: Deep Quality Control for Image-to-template Brain MRI Affine Registration2 Quality Control, Image-to-template registration, Deep Neural Network

Affine registration of one or several brain image(s) onto a common reference space is a necessary prerequisite for many image processing tasks, such as brain segmentation or functional analysis. Manual assessment of registration quality is a tedious and time-consuming task, especially in studies comprising a large amount of data. An automated and reliable quality control (QC) becomes mandatory. Moreover, the computation time of the QC must be also compatible with the processing of massive datasets. Therefore, an automated deep neural network approaches appear as a method of choice to automatically assess registration quality.

In the current study, a compact 3D convolutional neural network (CNN), referred to as RegQCNET, is introduced to quantitatively predict the amplitude of an affine registration mismatch between a registered image and a reference template. This quantitative estimation of registration error is expressed using metric unit system. Therefore, a meaningful task-specific threshold can be manually or automatically defined in order to distinguish usable and non-usable images.

The robustness of the proposed RegQCNET is first analyzed on lifespan brain images undergoing various simulated spatial transformations and intensity variations between training and testing. Secondly, the potential of RegQCNET to classify images as usable or non-usable is evaluated using both manual and automatic thresholds. During our experiments, automatic thresholds are estimated using several computerassisted classification models (logistic regression, support vector machine, naïve bayes and random forest) through cross-validation. To this end we used expert's visual quality control estimated on a lifespan cohort of 3953 brains. Finally, the RegQCNET accuracy is compared to usual image features such as image correlation coefficient and mutual information.

Results show that the proposed deep learning QC is robust, fast and accurate to estimate affine registration error in processing pipeline.

Introduction

A wide variety of processing pipelines have been proposed in the literature to make automatic brain image analysis possible. Spatial and intensity normalizations are usually necessary prerequisites for functional [START_REF] Cook | Camino: open-source diffusion-MRI reconstruction and processing[END_REF]) [START_REF] Song | REST: a toolkit for resting-state functional magnetic resonance imaging data processing[END_REF] or structural studies [START_REF] Wei | Quantitative analysis of MRI signal abnormalities of brain white matter with high reproducibility and accuracy[END_REF] [START_REF] Chaogan | DPARSF: a MATLAB toolbox for "pipeline" data analysis of resting-state fMRI[END_REF]. These are commonly achieved using suitable algorithms designed for image-to-template registration [START_REF] Tustison | Large-scale evaluation of ants and freesurfer cortical thickness measurements[END_REF]) [START_REF] Jenkinson | [END_REF]) [START_REF] Collins | Automatic 3d intersubject registration of MR volumetric data in standardized talairach space[END_REF], inhomogeneity correction [START_REF] Tustison | N4itk: Improved n3 bias correction[END_REF]) [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF], or intensity normalization [START_REF] Nyúl | New variants of a method of MRI scale standardization[END_REF]) [START_REF] Friston | Spatial registration and normalization of images[END_REF]. A visual human inspection of the data after each step of the processing pipeline is commonly employed to detect possible problems in the outputs. This visual quality control (QC) is unfortunately not feasible when a huge amount of imaging data is involved (typically more than several thousands scans). Consequently, with the rise of large-scale datasets, recent efforts are dedicated to the development of reliable QC methods to detect pipeline failures [START_REF] Alfaro-Almagro | Image processing and quality control for the first 10,000 brain imaging datasets from uk biobank[END_REF] [START_REF] Kim | The LONI QC system: A semi-automated, web-based and freely-available environment for the comprehensive quality control of neuroimaging data[END_REF]. Although out of the scope of the current study, an increasing interest in registration QC is noticed in radiation therapy [START_REF] Brock | Use of image registration and fusion algorithms and techniques in radiotherapy: Report of the AAPM radiation therapy committee task group no. 132[END_REF] [START_REF] Paganelli | Patient-specific validation of deformable image registration in radiation therapy: Overview and caveats[END_REF].

Fig. 1 summarizes the context of the current paper. Our study focuses on the image registration step, this step is a necessary prerequisite to co-register one or several brain scans onto a common space defined by a template image. In practice, misregistered data are inherently encountered and thus a registration QC is needed (red box in Fig. 1) [START_REF] Avants | A reproducible evaluation of ants similarity metric performance in brain image registration[END_REF]. A manual assessment is generally employed for QC and thus automatic methods have been developed to achieve this task. However, such manual strategy is time consuming. Random forest [START_REF] Hessam | Quantitative error prediction of medical image registration using regression forests[END_REF]) and convolutional neural network (CNN) [START_REF] Eppenhof | Error estimation of deformable image registration of pulmonary ct scans using convolutional neural networks[END_REF] have been proposed to quantify registration accuracy for both parametric (i.e., rigid, affine) and deformable registrations in chest CT scans. In the context of neuroimgaging, several methods have been proposed for MRI brain registration to template space. In [START_REF] Fonov | Deep learning of quality control for stereotaxic registration of human brain MRI, bioRxiv[END_REF]) a cross-entropy loss function is used as an objective function to train a deep neural network on a serie of 2D control images. This method produces qualitative estimation (i.e., good or not good) of rigid registration accuracy. In [START_REF] Bannister | Deep neural networks for quality assurance of image registration[END_REF]) and [START_REF] Dubost | Automated image registration quality assessment utilizing deep-learning based ventricle extraction in clinical data[END_REF], a DICE metric between transformed and original organ contours is proposed as a surrogate of registration quality. The use of an indirect metric (i.e., DICE) estimated on an auxiliary task (i.e., segmentation) does not provide direct quantitative information on registration accuracy. Moreover, this information can be corrupted by segmentation error that is a complex task by itself. Finally, these three methods produce metrics (binary decision or auxiliary DICE) that cannot express registration error in metric unit system. Consequently, no meaningful task-specific threshold on the misalignment amplitude can be defined by a user.

Our contribution is four-fold:

(i) A compact 3D CNN is introduced to quantitatively estimate the quality of an affine alignment between a brain MRI and a template. The proposed QC network, referred to as RegQCNET, is quantitative and can be expressed using metric system units. Moreover, an efficient and robust training procedure, based on simulated affine transformations, is proposed. The inputs of the designed CNN are: the registered image and the reference template. Besides, it is demonstrated that RegQCNET meets computational requirements related to massive processing.

(ii) The robustness of the proposed RegQCNET is analyzed on lifespan brain images undergoing various simulated spatial transformations and intensity variations between training and testing.

(iii) The potential of RegQCNET to classify images as usable or non-usable is evaluated using both manual and automatic thresholds. Automatic thresholding is evaluated using several computer-assisted classification models (logistic regression, support vector machine, naïve bayes and random forest) through a cross-validation procedure. To this end we used expert's visual quality control estimated on a lifespan cohort of 3953 brains as a gold standard.

(iv) The RegQCNET accuracy is compared to usual image features such as image correlation coefficient and mutual information.

Figure 1: General principle of the proposed Quality Control (QC) on image-to-template registration. The current study aims at providing a QC module (red bock) designed to detect misaligned images. Note that this module can potentially be fed into an additional correction module (outside the scope of the current study).

Materials and Methods

Datasets

Figure 2 details the processing sequence designed to generate the datasets involved in our experiments. Throughout this study, we used 3 datasets -1 for training and 2 for testing -referred to as "Simulated training dataset", "Simulated testing dataset" (these last two were built using a lifespan dataset, for which synthetic affine transformations were applied) and "Real testing dataset". First, all the native images have been downloaded from public databases. While considered as the native images, these images may contain specific preprocessing (i.e. NDAR includes defacing). Afterwards, all these images went through our preprocessing pipeline as described in the following. Each generated dataset is displayed as a gray block. Image processing tasks are reported with green blocks.

2.1.1. Preprocessing of input images. To ensure spatial and intensity normalization between images, we used a preprocessing pipeline. Consequently, all the scans involved in this study were first preprocessed beforehand using the volBrain pipeline [START_REF] Manjón | volbrain: An online MRI brain volumetry system[END_REF]. This pipeline is based on the following steps: i) denoising [START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF], ii) inhomogeneity correction [START_REF] Tustison | N4itk: Improved n3 bias correction[END_REF], iii) affine registration into the template space (181×217×181 voxels at 1×1×1 mm 3 , the ICBM 152 Atlas template was taken as reference for registration [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF])), iv) manual human assessment of the registration as described in [START_REF] Coupé | Towards a unified analysis of brain maturation and aging across the entire lifespan: A MRI analysis[END_REF]) [START_REF] Coupé | Lifespan changes of the human brain in alzheimer's disease[END_REF]), and v) tissue-based intensity normalization [START_REF] Manjón | Robust MRI brain tissue parameter estimation by multistage outlier rejection[END_REF]. Finally, image intensities were normalized using z-scoring using the mean and standard deviation from the complete image field-of-view.

2.1.2. Simulated training and testing datasets. 360 T1-weighted MRI of cognitively normal subjects were randomly selected under constraints from the dataset used in our previous BigData study on normal aging [START_REF] Coupé | Towards a unified analysis of brain maturation and aging across the entire lifespan: A MRI analysis[END_REF]. This dataset was based on 9 datasets publicly available (C-MIND, NDAR, ABIDE, ICBM, IXI, OASIS, AIBL, ADNI1 and ADNI2). From 1 to 90 years we selected 2 females (F) and 2 males (M) for each age (i.e., 2F and 2M of 1 year old, 2F and 2M of 2 years old and so on). Therefore, we obtained a balanced group with 50% of each gender uniformly distributed from 1 to 90 years. This balanced selection is done to limit bias introduction in training and testing datasets and to make our QC method robust to different age and gender.

All the 360 MRIs underwent a human quality control. Consequently, all these images were considered as correctly aligned with negligible residual registration mismatch. The RMSE was thus considered to be equal to 0. We are aware that these images are not perfectly aligned and that negligible errors might remain. Considering these remaining errors equal to epsilon or zero do not impact the rest of the study.

At this point, we have a set of 360 scans including spatial and intensity normalizations. This set was split into two separate datasets (stratified/no redundant data): 300 scans were used to build the "Simulated training dataset" and 60 scans were used to build the "Simulated testing dataset". This split was done under contraint to ensure well-balance of age and gender between both.

The "Simulated training dataset". One can expect that the training set has to be populated densely enough in terms of both anatomical inter-subject variability and simulated spatial deformations. During our experiments, RegQCNET was trained using N simulated scans (randomly selected with replacement from the 300 abovementioned scans). These scans were simulated using affine spatial transformations. In order not to favor large or small transformations, we force the RMSE distribution resulting from the simulated transformation to be uniform. We tested N-values in the following set: {100, 1000, 10000}. Note that the obtained training set may include several transformations for each patient when N > 300.

The "Simulated testing dataset". This set was composed of 500 scans (randomly selected with replacement from the 60 subjects selected to build the testing data). In this way, several transformations were applied to each patient (≈ 8 in average).

Simulated spatial affine transformations. To train and test the proposed method, 3D spatial transformations -composed by translation, rotation and scale variationswere simulated. Ranges for X-, Y-, Z-translations, rotations around X-, Y-, Z-axis, X-, Y-, Z-scaling factors are detailed in the experimental section below. The RMSE was calculated for each simulated transformation. Let MAX_RMSE be the upper limit of the simulated RMSE. A set of spatial transformations with a uniform RMSE distribution in the interval [0, MAX_RMSE] voxels was built. To this end, 3D spatial transformations composed by translation, rotation and scale variation, were simulated as follows:

• X-, Y-, Z-translations were randomly selected separately in the interval [-100, 100] voxels,

• Rotations around X-, Y-, Z-axis were randomly selected in the interval [-45, 45] degrees,

• X-, Y-, Z-scaling factors were randomly selected in the interval [0.5, 1.5] (a factor of 1 being equivalent to no scaling), These transformations were then applied using b-spline interpolation to the images (see Fig. 4).

2.1.3. The "Real testing dataset". The performance of the proposed RegQCNET was evaluated on a massive database (N = 3953) including cognitively normal patients, patients with Alzheimer's Disease (AD) and Mild Cognitive Impairment (MCI). These 3953 MRIs were the remaining subjects from the large-scale cohort used in [START_REF] Coupé | Lifespan changes of the human brain in alzheimer's disease[END_REF]) after removal of the 360 cognitively normal subjects used to build the simulated training and testing dataset. Consequently, the real testing dataset contained pathological alterations unseen in the training dataset. A visual assessment was done by checking screen shots of one sagittal, one coronal and one axial slice in middle of the 3D volume using the volBrain reports [START_REF] Manjón | volbrain: An online MRI brain volumetry system[END_REF]. Therefore, a human-based QC was available for all the scans and was used as qualitative ground truth.

Proposed RegQCNET

2.2.1. Implemented quantitative metric. In this study, we aim at quantifying the residual misalignment (noted T ) between two given images via the Root Mean Square Error criterion (RMSE) computed as follows [START_REF] Maurer | Registration of head volume images using implantable fiducial markers[END_REF]): The proposed RegQCNET is thus designed to predict registration RMSE using two given images: a reference template and a registered one.

RMSE = 1 |Ω| r∈Ω u( r) 2 + v( r) 2 + w( r) 2 (
2.2.2. Implemented deep neural network. Figure 3 describes the architecture of the proposed quantitative CNN-based QC for image-to-template registration. Input images were first down-sampled by a factor 4 (note that a down-sampling factor 2 was also tested and discussed). We used a convolutional encoder followed by a 3 regression layers per resolution level using a basis of 24 filters of 3×3×3 (i.e., 24 filters for the first layer, 48 for the second and so on). Each block was composed of batch normalization, convolution and ReLU activation. We employed the following parameters: batch size = 1, optimizer = Adam with default parameters, epoch = 100, loss = Mean Square Error (MSE) and dropout = 0.5 after each block. We used 2 input channels (the down-sampled T1w and the template images).

Experimental setup.

2.3.1. Assessment of RegQCNET on the "Simulated testing dataset". For this dataset, RegQCNET estimations were challenged against real RMSEs by evaluating R 2 , slope and Y-intercept of a linear regression. While the R 2 provides information about the precision of the proposed QC, the slope and Y-intercept quantifies its accuracy. We compared results obtained using different sizes for the training dataset (i.e., with N=100, 1000 and 10000, respectively).

To assess the robustness of our method to inaccuracy in image normalization and inhomogeneity correction we performed two experiments. First, we used a given uniform intensity shift to simulate normalization inconsistencies between images. Second, we used a non-uniform intensity bias to simulate inaccurate inhomogeneity correction within images.

Robustness against uniform intensity shift. To evaluate the robustness of RegQCNET to incorrect image normalization, RegQCNET was challenged against uniform intensity variation applied on all scans included in the testing dataset. For this purpose, the test scans were identically disturbed as follows: Let MAX_INTENSITY be the maximum intensity of an image I, and J be the image I after application of the spatially homogeneous intensity bias. For all voxel location r, the intensity I( r) was multiplied by a factor 2, while being restricted in the interval original intensity range [0, MAX_INTENSITY]:

J( r) = I( r) × 2 if I( r) × 2 < MAX_INTENSITY MAX_INTENSITY otherwise (2)
Robustness against non-uniform intensity bias. To evaluate the robustness of RegQCNET to error during inhomogenety correction, RegQCNET was challenged against a non-uniform intensity bias applied on all scans included in the testing dataset.

For this purpose, the test scans were identically disturbed as follows: Let r 0 = (x 0 , y 0 , z 0 ) be the voxel coordinates at the central position of an image I, and K be the image I after application of the spatially heterogeneous intensity bias. For all voxel location r, the intensity I( r) was weighted by a voxel-wise exponential decay as follow:

K( r) = I( r) × exp - r -r 0 2 2 2σ 2 (3)
Practically, intensities in voxels located close to the central position r 0 were less disturbed than those located further. In the scope of this study, we used σ = 60.

2.3.2. Assessment of RegQCNET on the "Real testing dataset". First, the 3953 brain images of the "Real testing dataset" were visually inspected to build a gold standard. 13 mis-registered brain images were detected and considered as a "negative case" for the rest of the manuscript. The "Simulated training dataset" with N = 10000 was here employed for training. The accuracies, area under the ROC curve (AUROC), sensitivity, specificity, positive predictive values (PPVs) and negative predictive values (NPVs) were recorded for the following experiments:

Manually defined threshold. RegQCNET was used to differentiate scans with RMSE higher than a user-defined threshold noted δ (δ was expressed in millimeters). We tested δ-values in the following set: {5, 10, 20, 50} mm. The AUROC was computed using the method detailed in [START_REF] Cantor | Determining the area under the ROC curve for a binary diagnostic test[END_REF].

Automatically defined threshold. A 10-fold-stratified cross-validation was used to evaluated the performance of RegQCNET when using an threshold automatically tuned by machine learning algorithm. The dataset of 3953 brains was randomly partitioned into training (90%) and testing (10%) subsets (making sure that at least one positive and one negative case were included in each subset). For this purpose, the following classification algorithms were applied using the commercial software Matlab ( 1994-2020 The MathWorks, Inc.)/"Statistics and Machine Learning" toolbox: logistic regression (LR), support vector machine (SVM), naïve bayes (NB) and random forest (RF). Default hyper-parameters in Matlab implementations were employed for RF and SVM (RF: Classification method/100 bagged decision trees, SVM: supports sequential minimal optimization/box constraint/linear kernel) [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF]. The cross-validation steps were repeated 1000 times with shuffling of the folds. Finally, average metrics, standard deviations and confident intervals were calculated.

Comparison with usual image features. Our automatically defined threshold experiment was also conducted using correlation coefficient (CC) and mutual information (MI) for comparison.

Hardware and implementation details.

We evaluated the computational burden of our proposed method using an Intel Xeon E5-2683 2.4 GHz (2 Hexadeca-core) with 256 GB of RAM equipped by a GPU Nvidia Tesla V100 with 16 GB of memory. The computation time during the testing session was evaluated without and with the use of the GPU. Our implementation was using Tensorflow 1.4 and Keras 2.2.4.

Results

Fig. 4 shows typical images generated using synthetic affine transformations. Middle transversal, coronal and sagittal slices are reported for several 3D volumes. The template used as reference for affine image registration (see section 2.1.1) is displayed in the first row. The second (scan #1) row shows 3D brain images from the original lifespan dataset (RMSE considered equal to 0 mm). Lower rows (scan # [3 -5]) show examples of 3D scans obtained after the application of simulated spatial transformations of various amplitudes, as described in section 2.1.2. Note that several images underwent different masking (see (f) on the second row) due to defacing (e.g., NDAR dataset).

3.1. Assessment of RegQCNET on the "Simulated testing dataset". Fig. 5 reports the precision and the accuracy of RegQCNET obtained on the "simulated testing dataset" (as described in section 2.1.2) using a training dataset of 100 (5a), 1000 (5b) and 10000 (5c) scans. As one can expect, the precision (rated by the R 2 of the linear fit) improves when the size of the training dataset increased. The R 2 converged slowly toward 1 along with the size of the training dataset increased (R 2 equal to 0.84, 0.95 and 0.99 were obtained using 100, 1000 and 10000 images, respectively). The accuracy (rated by the slope and the Y-intercept of the linear regression) followed the same trend. As long as N increased, the slope and the Y-intercept converged toward optimal values (i.e., 1 and 0, respectively). Fig. 6 shows the impact uniform intensity shift and spatially varying intensity bias (as described in section 2.1.2) on the performance of the proposed RegQCNET. While a uniform intensity shift did not deteriorate the precision (R 2 =0.99) and the accuracy (slope=0.98/Y-intercept=-1.08) (6c) this observation did not hold when a non-uniform bias (6e): in this last case, both precision (R 2 =0.98) and accuracy (slope=0.89/Y-intercept=7.54) were slightly worse.

Scan #1 RMSE=0 (d) (e) (f) Scan #2 RMSE=10 (g) (h) (i) Scan #3 RMSE=21 (j) (k) (l) Scan #4 RMSE=41 (m) (n) (o) Scan #5 RMSE=99 (p) (q) (r)
3.2. Assessment of RegQCNET on the "Real testing dataset".

Fig. 7 shows RMSE estimated by RegQCNET on the 3953 tested brain MRIs.

Transversal cross-section of mis-registered images, as detected by visual inspection, are reported above and below the graph. These images can be visually compared to the corresponding template cross-section reported in Fig. 4a. In two images (referred to as case #9 and case #11 in Fig. 7) of patients with Alzeihmer's Disease, very large lateral ventricle slightly disturbed RegQCNET (estimated RMSE < 20 mm). In the 9 other images, huge mis-registrations are observable, which have been detected by RegQCNET (estimated RMSE > 50 mm).

3.2.1. Manually defined threshold. Table 1 reports classification scores of RegQCNET using different manually defined thresholds.

A classification threshold δ = 10 mm provided best scores: accuracy=99.6%, AUROC=1.0, sensitivity=99.6%, specificity=100.0%, PPV=100.0%, and NPV=44.8% (note that NPV=44.8% means here that 16 good-registered images were considered as mis-registered). Good-and misregistered images were thus correctly identified in 3924/3940 and 13/13 brain images, respectively.

3.2.2. Automatically defined threshold. The RegQCNET output served as a metric in all tested machine learning classifiers (see Table 2). In particular, using a logistic regression classifier, QC scores were: accuracy=96.0%, AUROC=1.0, sensitivity=95.9%, specificity=100.0%, PPV=97.5%, and NPV=93.7%.

Comparison with usual image metrics.

Very poor scores were obtained using CC (best classifier=naïve bayes) and MI (best classifier=logistic regression): a good detection of correctly registered images was achievable by accepting a large amount of false-negative cases, as shown in Fig. 8. Conversely, a perfect detection of mis-registered images was only achievable by accepting a dramatic impact on the sensitivity (i.e., 1.7% and 30.9% for CC and MI, respectively, as shown in Table 2).

Discussion

The proposed method aims at quantifying the amplitude of the spatial affine mismatch between a brain MRI and a template. To this end, we used RMSE as criterion to evaluate affine registration quality. Our experimental results demonstrate that the proposed RegQCNET outperforms traditional intensity-based criteria. Moreover, using automatic threshold, our approach delivers reproducible results and minimizes operator dependency.

It can be observed that good scores were obtained thanks to the use of i) a well-balanced training population covering the entire lifespan and ii) a well-balanced distribution of the simulated transformations (i.e., uniform distribution of the resulting RMSEs). It is interesting to highlight that comparable precision and accuracy were obtained on simulated data (as described in section 2.1.2) using a subsampling factor 2 on the images (instead of the subsampling factor 4 used in the presented results).

One can distinguish two potential contributions on the apparent image-to-template mismatch: (i) the effective registration error that we aim to quantify and (ii) the anatomical variability between the actual image and the template (in particular, the size of the brain varies a lot along the lifespan). Any unobserved spatial transformations/brain shapes/image artifacts during the training step may disturb in turn the proposed quantitative CNN-based QC. This phenomenon can be observed in Fig. 5a where an insufficiently populated training dataset was employed (100 images). In turn, a dramatic impact arise on both precision and accuracy. Note that using correlation or mutual information as registration QC, any intensity variation between a given brain MRI and the template is attributed to an image mismatch. The tested machine learning classifiers thus provided very poor results in terms of accuracy, Another limitation arise when an intensity perturbation occurs between training and testing. While a spatially homogeneous bias did not impact the performance (Fig. 6c), an accurate correction of spatially heterogeneous intensity bias is a necessary prerequisite to obtain good performance when using RegQCNET (Fig. 6e) ( [START_REF] Nyúl | New variants of a method of MRI scale standardization[END_REF]. It has to be noted that our framework was robust to various degrees of masking due to defacing (e.g., NDAR dataset).

As one can expect, the range of spatial transformations during training has to be carefully determined. That brings us to an inherent limit of the proposed technique. Indeed, a RMSE extrapolation outside training limits is intrinsically not possible using our CNN-based approach and thus a large training range is mandatory. This limitation could be limited by training RegQCNET on a larger range of RMSE.

Using the proposed CNN-based QC, a few tenth of a seconds (700 ms and 200 ms without and with the use of GPU acceleration, respectively) is needed to provide a quantitative prediction of the image-to-template alignment accuracy. This perfectly meets our computational requirements related to the inclusion of this QC step in massive processing.

While our experimental results demonstrate that a compact network is an efficient solution to estimate the quality of affine registration between a T1-MRI and a template, the direct translation of our approach to radiation therapies is not straightforward. In the context of radiation therapies, several performance indicators and registration QC solutions have been proposed (see [START_REF] Paganelli | Scale invariant feature transform in adaptive radiation therapy: a tool for deformable image registration assessment and re-planning indication[END_REF]) and [START_REF] Brock | Use of image registration and fusion algorithms and techniques in radiotherapy: Report of the AAPM radiation therapy committee task group no. 132[END_REF]). Concerning abdominal organs, it must be underlined that deformations are non-rigid and thus extension of our framework to non-rigid registration would be required. For instance, recent works conducted in the abdomen used biomechanical criteria to assess image registration accuracy [START_REF] Zachiu | Anatomically plausible models and quality assurance criteria for online mono-and multi-modal medical image registration[END_REF] [START_REF] Zachiu | Biomechanical quality assurance criteria for deformable image registration algorithms used in radiotherapy guidance[END_REF]. Such approaches involve the estimation of mechanical stress, which would occur within the observed tissues. The calculated stress can then be compared to plausible physiological limits. We believe that a combination of these two complementary approaches (i.e., the network and the biomechanical strategies) should be investigated in future studies.

Conclusion

This study demonstrates that quantitative estimation of registration mismatch between a brain image and a template can be achieved using 3D CNNs. However, to ensure the quality of the estimation, the training dataset have to be carefully designed. To this end, in this study we used: i) a gender and age well-balanced lifespan dataset covering the entire lifespan, ii) an uniformly distributed amplitudes of random spatial transformations to cover registration error from 0 to 100 millimeters, and iii) a sufficient amplitude range of simulated spatial transformations. The proposed tool can be used as quality control for automated image registration of T1-weigthed brain onto a reference template. Future studies will include the extension of the proposed RegQCNET to complex elastics image deformations, the estimation of 3D RMSE maps, the impact of incomplete, noisy and corrupted brains, as well as the extension of the method to cross-contrast and multi-modal images.

Figure 2 :

 2 Figure 2: Processing sequence designed to generate datasets involved in experiments.Each generated dataset is displayed as a gray block. Image processing tasks are reported with green blocks.

  1) r = (x, y, z) being the voxel coordinates, Ω the image coordinates domain, |Ω| the number of voxels in Ω (|Ω| = 181 × 217 × 181 in the current study) and T = (u, v, w) the voxelwise 3D residual displacement vector field.

Figure 3 :

 3 Figure 3: Architecture of the proposed RegQCNET for image-to-template registration. Each block is composed of batch normalization, convolution and ReLU activation. The number of input channels (NC) as well as the number of 3 × 3 × 3 filters are indicated on the top of each block.

Figure 4 :

 4 Figure 4: Typical images generated using synthetic affine transformations. Transversal (left column), coronal (middle column) and sagittal (right column) central slice are reported for the template scan (a-c) and for different subjects/various RMSE values. Each generated image is referred to as "Scan #[1 -5]".

Figure 5 :

 5 Figure 5: RMSE estimated by RegQCNET on the simulated dataset using training set composed of 100 (a),1000 (b) and 10000 (c) scans. Estimated RMSEs are plotted against real RMSEs and the R 2 , slope and Y-intercept of a linear regression are reported in the insert of each graph.

Figure 6 :

 6 Figure 6: Results obtained on the simulation experiment when an intensity perturbation is applied on the testing dataset. A training dataset composed by 10000 scans (without intensity bias) was used. The axial slice of a brain scan is reported before (a) and after application of an uniform (b) and a non-uniform (d) intensity bias. Estimated RMSEs are plotted against real RMSEs for the uniform (c) and the non-uniform (e) bias, and the R 2 , slope and Y-intercept of a linear regression are reported.

Figure 7 :Figure 8 :

 78 Figure 7: RMSEs estimated by RegQCNET on the 3953 brains of the data base. Misregistered images (transversal), as detected by visual inspection, are reported above and below the graph. The corresponding template cross-section is shown in Fig. 4a.

Table 1 :

 1 Classification scores of the proposed RegQCNET on the "Real testing dataset". Quantitative scores (i.e., computer-assisted) are given, assuming the qualitative inspection (i.e., visual) as a gold standard. AUROC: area under the ROC curve; PPV: positive predictive value; NPV: negative predictive value. Accuracies, sensitivities, specificities, PPVs, and NPVs are shown in percentages. Best performance are reported in bold font.

	Tustison

et al. 2010)

Table 2 :

 2 Classification scores of the various classifiers on the "Real testing dataset". Quantitative scores were derived via evaluation of RegQCNET, correlation coefficient (CC), and mutual information (MI) (after 10-fold cross-validation) by various machinelearning algorithms (LR: logistic regression, SVM: support machine vector, NB: naïve bayes, RF: random forest). Quantitative indicators are shown with standard deviations and 95% confidence intervals in parentheses. Best performance are reported in bold font for each indicator.

Acknowledgment

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/). This work benefited from the support of the project DeepvolBrain of the French National Research Agency (ANR-18-CE45-0013). This study was achieved within the context of the Laboratory of Excellence TRAIL ANR-10-LABX-57 for the BigDataBrain project. Moreover, we thank the Investments for the future Program IdEx Bordeaux (ANR-10-IDEX-03-02, HL-MRI Project), Cluster of excellence CPU and the CNRS/INSERM for the DeepMultiBrain project. This study has been also supported by the DPI2017-87743-R grant from the Spanish Ministerio de Economia, Industria Competitividad. The authors gratefully acknowledge the support of NVIDIA Corporation with their donation of a TITAN X GPU used in this research.