Filippo Bonchi
email: filippo.bonchi@ens-lyon.fr

Damien Pous
email: damien.pous@ens-lyon.fr

Hacking Nondeterminism with Induction and Coinduction

. We compare our approach to the recently introduced antichain algorithms and we give concrete examples where we exponentially improve over antichains. Experimental results show significant improvements.

INTRODUCTION

Checking language equivalence of finite automata is a classic problem in computer science, with many applications in areas ranging from compilers to model checking.

Equivalence of deterministic finite automata (DFA) can be checked either via minimization [START_REF] Hopcroft | An n log n algorithm for minimizing in a finite automaton[END_REF] or through Hopcroft and Karp's algorithm [START_REF] Hopcroft | A linear algorithm for testing equivalence of finite automata[END_REF], which exploits an instance of what is nowadays called a coinduction proof principle [START_REF] Milner | Communication and Concurrency[END_REF][START_REF] Sangiorgi | Introduction to Bisimulation and Coinduction[END_REF][START_REF] Rutten | Automata and coinduction (an exercise in coalgebra)[END_REF]: two states are equivalent if and only if there exists a bisimulation relating them. In order to check the equivalence of two given states, Hopcroft and Karp's algorithm creates a relation containing them and tries to build a bisimulation by adding pairs of states to this relation: if it succeeds then the two states are equivalent, otherwise they are different.

On the one hand, minimization algorithms have the advantage of checking the equivalence of all the states at once, while Hopcroft and Karp's algorithm only checks a given pair of states. On the other hand, they have the disadvantage of needing the whole automata from the beginning, while Hopcroft and Karp's algorithm can be executed "onthe-fly" [START_REF] Fernandez | On-the-fly verification of finite transition systems[END_REF], on a lazy DFA whose transitions are computed on demand.

This difference is essential for our work and for other recently introduced algorithms based on antichains [START_REF] Wulf | Antichains: A new algorithm for checking universality of finite automata[END_REF][START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF]. Indeed, when starting from non-deterministic finite automata (NFA), determinization induces an exponential factor. In contrast, the algorithm we introduce in this work Extended Abstract, a full version of this paper is available in Proc. POPL, 2013, ACM. Work partially funded by the PiCoq (ANR-10-BLAN-0305) and PACE (ANR-12IS02001) projects.

http://doi.acm.org/10.1145/2713167. for checking equivalence of NFA (as well as those using antichains) usually does not build the whole deterministic automaton, but just a small part of it. We write "usually" because in few cases, the algorithm can still explore an exponential number of states.

Our algorithm is grounded on a simple observation on DFA obtained by determinizing an NFA: for all sets X and Y of states of the original NFA, the union (written +) of the language recognized by X (written [[X]]) and the language recognized by Y ([[Y ]]) is equal to the language recognized by the union of X and Y ([[X + Y ]]). In symbols:

[[X + Y ]] = [[X]] + [[Y ]] (1)
This fact leads us to introduce a sound and complete proof technique for language equivalence, namely bisimulation up to context, that exploits both induction (on the operator +) and coinduction: if a bisimulation R relates the set of states X1 with Y1, and X2 with Y2, then

[[X1]] = [[Y1]] and [[X2]] = [[Y2]

] and, by [START_REF] Abdulla | When simulation meets antichains[END_REF], we can immediately conclude that X1 + X2 and Y1 + Y2 are language equivalent as well. Intuitively, bisimulations up to context are bisimulations which do not need to relate X1 + X2 with Y1 + Y2 when X1 is already related with Y1 and X2 with Y2.

To illustrate this idea, let us check the equivalence of states x and u in the following NFA. (Final states are overlined, labeled edges represent transitions.)

x a z a o o a & & y a f f u a ( ( a w a g g v a o o
The determinized automaton is depicted below. The dashed lines (numbered by 1, 2, 3) form a smaller relation which is not a bisimulation, but a bisimulation up to context: the equivalence of {x, y} and {u, v, w} is deduced from the fact that {x} is related with {u} and {y} with {v, w}, without the need to further explore the automaton.

{x} a /
Bisimulations up-to, and in particular bisimulations up to context, have been introduced in the setting of concurrency theory [START_REF] Milner | Communication and Concurrency[END_REF][START_REF] Sangiorgi | On the bisimulation proof method[END_REF] as a proof technique for bisimilarity of CCS or π-calculus processes. As far as we know, they have never been used for proving language equivalence of NFA.

Among these techniques one should also mention bisimulation up to equivalence, which, as we show in this paper, is implicitly used in Hopcroft and Karp's original algorithm. This technique can be explained by noting that not all bisimulations are equivalence relations: it might be the case that a bisimulation relates X with Y , and Y with Z, but not X with Z. However, since [

[X]] = [[Y ]] and [[Y ]] = [[Z]],
we can immediately conclude that X and Z recognize the same language. Analogously to bisimulations up to context, a bisimulation up to equivalence does not need to relate X with Z when they are both related with some Y .

The techniques of up-to equivalence and up-to context can be combined, resulting in a powerful proof technique which we call bisimulation up to congruence. Our algorithm is in fact just an extension of Hopcroft and Karp's algorithm that attempts to build a bisimulation up to congruence instead of a bisimulation up to equivalence. An important property when using up to congruence is that we do not need to build the whole deterministic automata. For instance, in the above NFA, the algorithm stops after relating z with u + w and does not build the remaining states. Despite their use of the up to equivalence, this is not the case with Hopcroft and Karp's algorithm, where all accessible subsets of the deterministic automata have to be visited at least once.

The ability of visiting only a small portion of the determinized automaton is also the key feature of the antichain algorithm [START_REF] Wulf | Antichains: A new algorithm for checking universality of finite automata[END_REF] and its optimization exploiting similarity [START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF]. The two algorithms are designed to check language inclusion rather than equivalence and, for this reason, they do not exploit equational reasoning. As a consequence, the antichain algorithm usually needs to explore more states than ours. Moreover, we show how to integrate the optimization proposed in [START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF] in our setting, resulting in an even more efficient algorithm.

Outline

Section 2 recalls Hopcroft and Karp's algorithm for DFA, showing that it implicitly exploits bisimulation up to equivalence. Section 3 describes the novel algorithm, based on bisimulations up to congruence. We compare this algorithm with the antichain one in Section 4.

DETERMINISTIC AUTOMATA

A deterministic finite automaton (DFA) over the alphabet A is a triple (S, o, t), where S is a finite set of states, o : S → 2 is the output function, which determines if a state x ∈ S is final (o(x) = 1) or not (o(x) = 0), and t : S → S A is the transition function which returns, for each state x and for each letter a ∈ A, the next state ta(x). Any DFA induces a function [ This problem generalizes the familiar problem of checking whether two automata accept the same language: just take the union of the two automata as the automaton (S, o, t), and determine whether their respective starting states are language equivalent.

Language equivalence via coinduction

We first define bisimulation. We make explicit the underlying notion of progression, which we need in the sequel.

Definition 1 (Progression, Bisimulation). Given two re- lations R, R ⊆ S 2 on states, R progresses to R , denoted R R , if whenever x R y then 1. o(x) = o(y) and
2. for all a ∈ A, ta(x) R ta(y).

A bisimulation is a relation R such that R R.

As expected, bisimulation is a sound and complete proof technique for checking language equivalence of DFA: Proposition 1 (Coinduction). Two states are language equivalent iff there exists a bisimulation that relates them.

Naive algorithm

Figure 1 shows a naive version of Hopcroft and Karp's algorithm for checking language equivalence of the states x and y of a deterministic finite automaton (S, o, t). Starting from x and y, the algorithm builds a relation R that, in case of success, is a bisimulation. Proposition 2. For all x, y ∈ S, x ∼ y iff Naive(x, y).

Proof. We first observe that if Naive(x, y) returns true then the relation R that is built before arriving to step 4 is a bisimulation. Indeed, the following proposition is an invariant for the loop corresponding to step 3:

R R ∪ todo

Since todo is empty at step 4, we have R R, i.e., R is a bisimulation. By Prop. 1, x ∼ y. On the other hand, Naive(x, y) returns false as soon as it finds a word which is accepted by one state and not the other. For example, consider the DFA with input alphabet A = {a} in the left-hand side of Figure 2, and suppose we want to check that x and u are language equivalent.

During the initialization, (x,

u) is inserted in todo. At the first iteration, since o(x) = 0 = o(u), (x, u) is inserted in R and (y, v) in todo. At the second iteration, since o(y) = 1 = o(v), (y, v) is inserted in R and (z, w) in todo. At the third iteration, since o(z) = 0 = o(w), (z, w) is inserted in R and (y, v) in todo.
At the fourth iteration, since (y, v) is already in R, the algorithm does nothing. Since there are no more pairs to check in todo, the relation R is a bisimulation and the algorithm terminates returning true.

These iterations are concisely described by the numbered dashed lines in Figure 2. The line i means that the connected pair is inserted in R at iteration i. (In the sequel, when enumerating iterations, we ignore those where a pair from todo is already in R so that there is nothing to do.)

In the previous example, todo always contains at most one pair of states but, in general, it may contain several of them. We do not specify here how to choose the pair to extract in step 3.1; we discuss this point in Section 3.2.

Hopcroft and Karp's algorithm

The naive algorithm is quadratic: a new pair is added to R at each non-trivial iteration, and there are only n 2 such pairs, where n = |S| is the number of states of the DFA. To make this algorithm (almost) linear, Hopcroft and Karp actually record a set of equivalence classes rather than a set of visited pairs. As a consequence, their algorithm may stop earlier it encounters a pair of states that is not already in R but belongs to its reflexive, symmetric, and transitive closure. For instance, in the right-hand side example from Figure 2, we can stop when we encounter the dotted pair (y, w) since these two states already belong to the same equivalence class according to the four previous pairs. With this optimization, the produced relation R contains at most n pairs. Formally, ignoring the concrete data structure used to store equivalence classes, Hopcroft and Karp's algorithm consists in replacing step 3.2 in Figure 1 with

( 3 . 2 ) i f (x , y ) ∈ e(

R) t h e n c o n t i n u e ;

where e : P(S 2 ) → P(S 2 ) is the function mapping each relation R ⊆ S 2 into its symmetric, reflexive, and transitive closure. We refer to this algorithm as HK.

Bisimulations up-to

We now show that the optimization used by Hopcroft and Karp corresponds to exploiting an "up-to technique".

Definition 2 (Bisimulation up-to). Let f : P(S 2 ) → P(S 2 ) be a function on relations. A relation R is a bisimulation up to f if R f (R), i.e., if x R y, then 1. o(x) = o(y) and 2. for all a ∈ A, ta(x) f (R) ta(y).
With this definition, Hopcroft and Karp's algorithm just consists in trying to build a bisimulation up to e. To prove the correctness of the algorithm, it suffices to show that any bisimulation up to e is contained in a bisimulation. To this end, we the notion of compatible function [START_REF] Sangiorgi | On the bisimulation proof method[END_REF][START_REF] Pous | Complete lattices and up-to techniques[END_REF]:

Definition 3 (Compatible function). A function f : P(S 2 ) → P(S 2 ) is compatible if it is monotone and it preserves pro- gressions: for all R, R ⊆ S 2 , R R entails f (R) f (R ) .
Proposition 3. Let f be a compatible function. Any bisimulation up to f is contained in a bisimulation.

We could prove directly that e is a compatible function; we however take a detour to ease our correctness proof for the algorithm we propose in Section 3.

Lemma 1. The following functions are compatible: id: the identity function; f • g: the composition of compatible functions f and g; F : the pointwise union of an arbitrary family F of compatible functions:

F (R) = f ∈F f (R); f ω : the (omega) iteration of a compatible function f , de- fined by f ω = i f i , with f 0 = id and f i+1 = f • f i ; r: the constant reflexive function: r( ) = {(x, x) | x ∈ S}; s: the converse function: s(R) = {(y, x) | x R y}; t: the squaring function: t(R) = {(x, z) | ∃y, x R y R z}.
Intuitively, given a relation R, (s∪id)(R) is the symmetric closure of R, (r ∪ s ∪ id)(R) is its reflexive and symmetric closure, and (r ∪ s ∪ t ∪ id) ω (R) is its symmetric, reflexive and transitive closure: e = (r ∪ s ∪ t ∪ id) ω . Another way to understand this decomposition of e is to recall that e(R) can be defined inductively by the following rules: Proof. Same proof as for Prop. 2, by using the invariant R e(R) ∪ todo. We deduce that R is a bisimulation up to e after the loop. We conclude with Thm. 1 and Prop. 1.

Returning to the right-hand side example from Figure 2, Hopcroft and Karp's algorithm constructs the relation

RHK = {(x, u), (y, v), (z, w), (z, v)}
which is not a bisimulation, but a bisimulation up to e: it contains the pair (x, u), whose b-transitions lead to (y, w), which is not in RHK but in its equivalence closure, e(RHK).

Naive(X, Y )

( 1 ) R i s e m p t y ; todo i s e m p t y ;

( 2 ) i n s e r t (X, Y ) i n todo; [START_REF] Bartels | Generalised coinduction[END_REF]

NON-DETERMINISTIC AUTOMATA

We now move from DFA to non-deterministic automata (NFA). A NFA over the alphabet A is a triple (S, o, t), where S is a finite set of states, o : S → 2 is the output function, and t : S → P(S) A is the transition relation: it assigns to each state x ∈ S and letter a ∈ A a set of possible successors.

The powerset construction transforms any NFA (S, o, t) into the DFA (P(S), o , t ) where o : P(S) → 2 and t : P(S) → P(S) A are defined for all X ∈ P(S) and a ∈ A as follows:

o (X) =      o(x) if X = {x} with x ∈ S 0 if X = 0 o (X1) + o (X2) if X = X1 + X2 t a (X) =      ta(x) if X = {x} with x ∈ S 0 if X = 0 t a (X1) + t a (X2) if X = X1 + X2
(Here we use the symbol + to denote both set-theoretic union and Boolean or; similarly, we use 0 to denote both the empty set and the Boolean 'false'.) Observe that in (P(S), o , t ), the states form a semi-lattice (P(S), +, 0), and o and t are, by definition, semi-lattices homomorphisms. These properties are fundamental for the up-to technique we are going to introduce. In order to stress the difference with generic DFA, which usually do not carry this structure, we use the following definition.

a B B a & & y a f f z a o o x a / / y + z a / / x + y a / / x + y + z a F F
In the determinized NFA, x makes one single a-transition into y + z. This state is final: Algorithms for NFA can be obtained by computing the determinized NFA on-the-fly [START_REF] Fernandez | On-the-fly verification of finite transition systems[END_REF]: starting from the algorithms for DFA (Figure 1), it suffices to work with sets of states, and to inline the powerset construction. The corresponding code is given in Figure 3. The naive algorithm (Naive) does not use any up to technique, Hopcroft and Karp's algorithm (HK) reasons up to equivalence in step 3.2.

o (y + z) = o (y) + o (z) = x a E E a & & y a g g z a o o u a D D x a / / 1 y + z a / / 2 x + y a / / 3 x + y + z a F F 4 u a D D

Bisimulation up to congruence

The semi-lattice structure (P(S), +, 0) carried by determinized NFA makes it possible to introduce a new up-to technique, which is not available with plain DFA: up to congruence. This technique relies on the following function.

Definition 5 (Congruence closure). Let u : P(P(S) 2 ) → P(P(S) 2 ) be the function on relations on sets of states defined for all R ⊆ P(S) 2 as:

u(R) = {(X1 + X2, Y1 + Y2) | X1 R Y1 and X2 R Y2} . The function c = (r ∪ s ∪ t ∪ u ∪ id) ω is called the congruence closure function.
Intuitively, c(R) is the smallest equivalence relation which is closed with respect to + and which includes R. It could alternatively be defined inductively using the rules r, s, t, and id from the previous section, and the following one:

X1 c(R) Y1 X2 c(R) Y2 X1 + X2 c(R) Y1 + Y2 u Definition 6 (Bisimulation up to congruence). A bisimu- lation up to congruence for an NFA (S, o, t) is a relation R ⊆ P(S) 2 , such that whenever X R Y then 1. o (X) = o (Y ) and 2. for all a ∈ A, t a (X) c(R) t a (Y ).
Lemma 2. The function u is compatible.

Theorem 2. Any bisimulation up to congruence is contained in a bisimulation.

We already gave in Introduction an example of bisimulation up to context, which is a particular case of bisimulation up to congruence (up to context means up to (r ∪ u ∪ id) ω , without closing under s and t). ∈ e(R). However, R is a bisimulation up to congruence. Indeed, we have (x + y, u) ∈ c(R):

x + y c(R) u + y ((x, u) ∈ R) c(R) y + z + y ((y + z, u) ∈ R) = y + z c(R) u ((y + z, u) ∈ R)
In contrast, we need four pairs to get a bisimulation up to equivalence containing (x, u): this is the relation depicted with both dashed and dotted lines in Figure 4.

Note that we can deduce many other equations from R; in fact, c(R) defines the following partition of sets of states: {0}, {y}, {z}, {x, u, x+y, x+z, and the 9 remaining subsets}.

Optimized algorithm for NFA

The optimized algorithm, called HKC in the sequel, relies on up to congruence: step 3.2 from Figure 3 becomes

( 3 . 2 ) i f (X , Y ) ∈ c(R ∪ todo) t h

e n c o n t i n u e ;

Observe that we use c(R∪todo) rather than c(R): this allows us to skip more pairs, and this is safe since all pairs in todo will eventually be processed.

Corollary 2. For all X, Y ∈ P(S), X ∼ Y iff HKC(X, Y ).

Proof. Same proof as for Proposition 2, by using the invariant R c(R ∪ todo) for the loop. We deduce that R is a bisimulation up to congruence after the loop. We conclude with Theorem 2 and Proposition 1.

The most important point about these three algorithms is that they compute the states of the determinized NFA lazily. This means that only accessible states need to be computed, which is of practical importance since the determinized NFA can be exponentially large. In case of a negative answer, the three algorithms stop even before all accessible states have been explored; otherwise, if a bisimulation (possibly up-to) is found, it depends on the algorithm:

• With Naive, all accessible states need to be visited, by definition of bisimulation.

• With HK, the only case where some accessible states can be avoided is when a pair (X, X) is encountered: the algorithm skips this pair so that the successors of X are not necessarily computed (this situation never happens when starting with disjoint automata). In the other cases where a pair (X, Y ) is skipped, X and Y are necessarily already related with some other states in R, so that their successors will eventually be explored.

• With HKC, accessible states are often skipped. For a simple example, let us execute HKC on the NFA from Figure 4. After two iterations, R = {(x, u), (y + z, u)}. Since x + y c(R) u, the algorithm stops without building the states x + y and x + y + z. Similarly, in the example from the Introduction, HKC does not construct the four states corresponding to pairs 4, 5, and 6.

This ability of HKC to ignore parts of the determinized NFA can bring an exponential speed-up. For an example, consider the family of NFA in Figure 5, where n is an arbitrary natural number. Taken together, the states x and y are equivalent to z: they recognize the language (a+b) (a+b) n+1 . Alone, x recognizes the language (a+b) a(a+b) n , which is known for having a minimal DFA with 2 n states. Therefore, checking x+y ∼ z via minimization (as in [START_REF] Hopcroft | An n log n algorithm for minimizing in a finite automaton[END_REF]) requires exponential time, and the same holds for Naive and HK since all accessible states must be visited. This is not the case with HKC, which requires only polynomial time on this example. Indeed, HKC(x+y, z) builds the relation

R = {(x + y, z)} ∪ {(x + Yi + yi+1, Zi+1) | i < n} ∪ {(x + Yi + xi+1, Zi+1) | i < n} ,
where Yi = y+y1+ . . . +yi, and Zi = z+z1+ . . . +zi. R only contains 2n+1 pairs and is a bisimulation up to congruence. To see this, consider the pair (x+y+x1+y2, Z2) obtained from (x+y, z) after reading the word ba. Although this pair does not belong to R , it belongs to its congruence closure:

x+y+x1+y2 c(R ) Z1+y2 (x+y+x1 R Z1) c(R ) x+y+y1+y2 (x+y+y1 R Z1) c(R ) Z2 . (x+y+y1+y2 R Z2)
Remark 1. In the above derivation, the use of transitivity is crucial: R is a bisimulation up to congruence, but not a bisimulation up to context. In fact, there exists no bisimulation up to context of linear size proving x + y ∼ z.

We now discuss the exploration strategy, i.e., how to choose the pair to extract from the set todo in step 3.1. When looking for a counter-example, such a strategy has a large influence: a good heuristic can help in reaching it directly, while a bad one might lead to explore exponentially many pairs first. In contrast, the strategy does not impact much looking for an equivalence proof (when the algorithm eventually returns true). Actually, one can prove that the number of steps performed by Naive and HK in such a case does not depend on the strategy. This is not the case with HKC: the strategy can induce some differences. However, we experimentally observed that breadth-first and depth-first strategies usually behave similarly on random automata. This behaviour is due to the fact that we check congruence w.r.t. R ∪ todo rather than just R (step 3.2): with this optimization, the example above is handled in polynomial time whatever the chosen strategy. In contrast, without this small optimization, it requires exponential time with a depth-first strategy.

Computing the congruence closure

For the optimized algorithm to be effective, we need a way to check whether some pairs belong to the congruence closure of a given relation (step 3.2). We present a simple solution based on set rewriting; the key idea is to look at each pair (X, Y ) in a relation R as a pair of rewriting rules:

X → X + Y Y → X + Y ,
which can be used to compute normal forms for sets of states. Indeed, by idempotence, X R Y entails X c(R) X + Y .

Definition 7. Let R ⊆ P(S) 2 be a relation on sets of states. We define ;R ⊆ P(S) 2 as the smallest irreflexive relation that satisfies the following rules:

X R Y X ;R X + Y X R Y Y ;R X + Y Z ;R Z U + Z ;R U + Z Lemma 3.
For all relations R, ;R is confluent and normalizing.

In the sequel, we denote by X↓R the normal form of a set X w.r.t. ;R. Intuitively, the normal form of a set is the largest set of its equivalence class. Recalling the example from Figure 4, the common normal form of x + y and u can be computed as follows (R is the relation {(x, u), (y +z, u)}):

x + y * * u w w x + y + u + + x + u s s x + y + z + u
Theorem 3. For all relations R, and for all X, Y ∈ P(S),

we have X↓R = Y ↓R iff (X, Y ) ∈ c(R).
We actually have X↓R = Y ↓R iff X ⊆ Y ↓R and Y ⊆ X↓R, so that the normal forms of X and Y do not necessarily need to be fully computed in practice. Still, the worst-case complexity of this sub-algorithm is quadratic in the size of the relation R (assuming we count the number of operations on sets: unions and inclusion tests).

Note that many algorithms were proposed in the literature to compute the congruence closure of a relation (see, e.g., [START_REF] Nelson | Fast decision procedures based on congruence closure[END_REF][START_REF] Shostak | Deciding combinations of theories[END_REF][START_REF] Bachmair | Congruence closure modulo associativity and commutativity[END_REF]). However, they usually consider uninterpreted symbols or associative and commutative symbols, but not associative, commutative, and idempotent symbols, which is what we need here.

Using HKC for checking language inclusion

For NFA, language inclusion can be reduced to language equivalence: the semantics function [[-]] is a semi-lattice homomorphism, so that for all sets of states X, Y ,

[[X+Y ]] = [[Y ]] iff [[X]] + [[Y ]] = [[Y ]] iff [[X]] ⊆ [[Y ]]. Therefore, it suffices to run HKC(X+Y, Y ) to check the inclusion [[X]] ⊆ [[Y ]].
In such a situation, all pairs that are eventually manipulated by HKC have the shape (X +Y , Y ) for some sets X , Y . Step 3.2 of HKC can thus be simplified. First, the pairs in the current relation only have to be used to rewrite from right to left. Second, the following lemma shows that we do not necessarily need to compute normal forms: Lemma 4. For all sets X, Y and for all relations R, we have X+Y c(R) Y iff X ⊆ Y ↓R.

At this point, the reader might wonder whether checking the two inclusions separately is more convenient than checking the equivalence directly. This is not the case: checking the equivalence directly actually allows one to skip some pairs that cannot be skipped when reasoning by double inclusion. As an example, consider the DFA on the right of Figure 2. The relation computed by HKC(x, u) contains only four pairs (because the fifth one follows from transitivity). Instead, the relations built by HKC(x, x+u) and HKC(u+x, u) would both contain five pairs: transitivity cannot be used since our relations are now oriented (from y ≤ v, z ≤ v and z ≤ w, we cannot deduce y ≤ w). Figure 5 shows another example, where we get an exponential factor by checking the equivalence directly rather than through the two inclusions: transitivity, which is crucial to keep the relation computed by HKC(x+y, z) small (see Remark 1), cannot be used when checking the two inclusions separately.

In a sense, the behaviour of the coinduction proof method here is similar to that of standard proofs by induction, where one often has to strengthen the induction predicate to get a (nicer) proof.

Exploiting Similarity

Looking at the example in Figure 5, a natural idea would be to first quotient the automaton by graph isomorphism. By doing so, one would merge the states xi, yi, zi, and one would obtain the following automaton, for which checking x+y ∼ z is much easier. : :

As shown in [START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF] for antichain algorithms, one can do better, by exploiting any preorder contained in language inclusion. Hereafter, we show how this idea can be embedded in HKC, resulting in an even stronger algorithm.

For the sake of clarity, we fix the preorder to be similarity [START_REF] Milner | Communication and Concurrency[END_REF], which can be computed in quadratic time [START_REF] Henzinger | Computing simulations on finite and infinite graphs[END_REF]. To exploit similarity pairs in HKC, it suffices to notice that for any similarity pair x y, we have x+y ∼ y. Let denote the relation {(x+y, y) | x y}, let r denote the constant-to-function, and let c = (r ∪s∪t∪u∪id) ω . Accordingly, we call HKC' the algorithm obtained from HKC (Figure 3) by replacing (X, Y ) ∈ c(R ∪ todo) with (X, Y ) ∈ c (R ∪ todo) in step 3.2. The latter test can be reduced to rewriting thanks to Theorem 3 and the following lemma.

ANTICHAIN ALGORITHMS

Even though the problem of deciding NFA equivalence is PSPACE-complete [START_REF] Meyer | Word problems requiring exponential time[END_REF], neither HKC nor HKC' are in PSPACE: both of them keep track of the states they explored in the determinized NFA, and there can be exponentially many such states. This also holds for HK and for the more recent antichain algorithm [START_REF] Wulf | Antichains: A new algorithm for checking universality of finite automata[END_REF] (called AC in the following) and its optimization (AC') exploiting similarity [START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF]. The latter algorithms can be explained in terms of coinductive proof techniques: we establish in [START_REF] Bonchi | Extended version of this abstract, with omitted proofs, and web appendix for this work[END_REF] that they actually construct bisimulations up to context, i.e., bisimulations up to congruence for which one does not exploit symmetry and transitivity.

Theoretical comparison. We compared the various algorithms in details in [START_REF] Bonchi | Extended version of this abstract, with omitted proofs, and web appendix for this work[END_REF]. Their relationship is summarised in Figure 6, where an arrow X→Y means that (a) Y can explore exponentially fewer states than X, and (b) Y can mimic X, i.e., the coinductive proof technique underlying Y is at least as powerful as the one of X.

In the general case, AC needs to explore much more states than HKC: the use of transitivity, which is missing in AC, allows HKC to drastically prune the exploration. For instance, to check x+y ∼ z in Figure 5, HKC only needs a linear number of states (see Remark 1), while AC needs exponentially many states. In contrast, in the special case where one checks for the inclusion of disjoint automata, HKC and AC exhibit the same behaviour. Indeed, HKC cannot make use of transitivity in such a situation, as explained in Section 3.4. Things change when comparing HKC' and AC': even for checking inclusion of disjoint automata, AC' cannot always mimic HKC': the use of similarity tends to virtually merge states, so that HKC' can use the up to transitivity technique which AC' lack.

Experimental comparison. The theoretical relationships drawn in Figure 6 are substantially confirmed by an empirical evaluation of the performance of the algorithms. Here, we only give a brief overview; see [START_REF] Bonchi | Extended version of this abstract, with omitted proofs, and web appendix for this work[END_REF] for a complete description of those experiments.

We compared our OCaml implementation [START_REF] Bonchi | Extended version of this abstract, with omitted proofs, and web appendix for this work[END_REF] for HK, HKC and HKC', and the libvata C++ library [START_REF] Lengál | Vata: A library for efficient manipulation of non-deterministic tree automata[END_REF] for AC and AC'. We use a breadth-first exploration strategy: we represent the set todo from Figure 3 as a FIFO queue. As mentioned at the end of Section 3.2, considering a depth-first strategy here does not alter the behaviour of HKC in a noticeable way.

We performed experiments using both random automata and a set of automata arising from model-checking problems.

• Random automata. We used Tabakov and Vardi's model [START_REF] Tabakov | Experimental evaluation of classical automata constructions[END_REF] to generate 1000 random NFA with two letters and a given number of states. We executed all algorithms on these NFA, and we measured the number of processed pairs, i.e., the number of required iterations (like HKC, AC is a loop inside which pairs are processed). We observe that HKC improves over AC by one order of magnitude, and AC improves over HK by two orders of magnitude. Using up-to similarity (HKC' and AC') does not improve much; in fact, similarity is almost the identity relation on such random automata. The corresponding distributions for HK, HKC, and AC are plotted on Figure 7, for automata with 100 states. Note that while HKC only improves by one order of magnitude over AC when considering the average case, it improves by several orders of magnitude when considering the worst cases.

• Model checking automata. Abdulla et al. [START_REF] Abdulla | When simulation meets antichains[END_REF][START_REF] Doyen | Antichain Algorithms for Finite Automata[END_REF] used automata sequences arising from regular model-checking experiments [START_REF] Bouajjani | Abstract regular model checking[END_REF] to compare their algorithm (AC') against AC. We reused these sequences to test HKC' against AC' in a concrete scenario. For all those sequences, we checked the inclusions of all consecutive pairs, in both directions. The timings are given in Table 1, where we report the median values (50%), the last deciles (90%), the last percentiles (99%), and the maximum values (100%). We distinguish between the experiments for which a counter-example was found, and those for which the inclusion did hold. For HKC' and AC', we display the time required to compute similarity on a separate line: this preliminary step is shared by the two algorithms. As expected, HKC and AC roughly behave the same: we test inclusions of disjoint automata. HKC' is however quite faster than AC': up to transitivity can be exploited thanks to similarity pairs. Also note that over the 546 positive answers, 368 are obtained immediately by similarity.

CONCLUSIONS

Our implementation of HKC is available online [START_REF] Bonchi | Extended version of this abstract, with omitted proofs, and web appendix for this work[END_REF], together with proofs mechanized in the Coq proof assistant and an interactive applet making it possible to test the presented algorithms online, on user-provided examples.

Several notions analogous to bisimulations up to congruence can be found in the literature. For instance, selfbisimulations [START_REF]Graphes canoniques de graphes algébriques[END_REF][START_REF] Hirshfeld | A polynomial algorithm for deciding bisimilarity of normed context-free processes[END_REF] have been used to obtain decidability and complexity results about context-free processes. The main difference with bisimulation up to congruence is that self-bisimulations are proof techniques for bisimilarity rather than language equivalence. Other approaches, that are independent from the equivalence (like bisimilarity or language) are shown in [START_REF] Lenisa | From set-theoretic coinduction to coalgebraic coinduction: some results, some problems[END_REF][START_REF] Bartels | Generalised coinduction[END_REF][START_REF] Pous | Complete lattices and up-to techniques[END_REF]. These papers propose very general frameworks into which our up to congruence technique fits as a very special case. However, to our knowledge, bisimu- 1: Timings, in seconds, for language inclusion of disjoint NFA generated from model-checking. lation up to congruence has never been proposed before as a technique for proving language equivalence of NFA.

We conclude with directions for future work.

Complexity. The presented algorithms, as well as those based on antichains, have exponential complexity in the worst case while they behave rather well in practice. For instance, in Figure 7, one can notice that over a thousand random automata, very few require to explore a large amount of pairs. This suggests that an accurate analysis of the average complexity might be promising. An inherent problem comes from the difficulty to characterise the average shape of determinised NFA [START_REF] Tabakov | Experimental evaluation of classical automata constructions[END_REF]. To avoid this problem, with HKC, we could try to focus on the properties of congruence relations. For instance, given a number of states, how long can be a sequence of (incrementally independent) pairs of sets of states whose congruence closure collapses into the full relation? (This number is an upper-bound for the size of the relations produced by HKC.) One can find ad-hoc examples where this number is exponential, but we suspect it to be rather small in average.

Model checking.

The experiments summarized in Table 1 show the efficiency of our approach for regular model-checking using automata on finite words. As in the case of antichains, our approach extends to automata on finite trees. We plan to implement such a generalisation and link it with tools performing regular tree model-checking.

In order to face other model-checking problems, it would be useful to extend up-to techniques to automata on infinite words, or trees. Unfortunately, the determinisation of these automata (the so called Safra's construction) does not seem suitable for exploiting neither antichains nor up to congruence. However, for some problems like LTL realizability [START_REF] Filiot | An antichain algorithm for ltl realizability[END_REF] that can be solved without prior determinization (the socalled Safraless approaches), antichains have been crucial in obtaining efficient procedures. We leave as future work to explore whether up-to techniques could further improve such procedures.

Figure 1 :

 1 Figure 1: Naive algorithm for checking the equivalence of states x and y of a DFA (S, o, t). The code of HK(x, y) is obtained by replacing the test in step 3.2 with (x , y ) ∈ e(R).

Figure 2 :

 2 Figure 2: Checking for DFA equivalence.

Corollary 1 .

 1 x e(R) x r x e(R) y y e(R) x s x e(R) y y e(R) zx e(R) z tx R y x e(R) y id Theorem 1. Any bisimulation up to e is contained in a bisimulation. For all x, y ∈ S, x ∼ y iff HK(x, y).

Figure 3 :

 3 Figure 3: On-the-fly naive algorithm, for checking the equivalence of sets of states X and Y of an NFA (S, o, t). HK(X, Y ) is obtained by replacing the test in step 3.2 with (X , Y ) ∈ e(R), and HKC(X, Y ) is obtained by replacing it with (X , Y ) ∈ c(R ∪ todo).

Definition 4 .

 4 A determinized NFA is a DFA (P(S), o , t ) obtained via the powerset construction of some NFA (S, o, t). Hereafter, we use a new notation for representing states of determinized NFA: in place of the singleton {x}, we just write x and, in place of {x1, . . . , xn}, we write x1 + • • • + xn. Consider for instance the NFA (S, o, t) depicted below (left) and part of the determinized NFA (P(S), o , t ) (right).

  x

Figure 4 :

 4 Figure 4: A bisimulation up to congruence.

Figure 4

 4 shows a more involved example illustrating the use of all ingredients of the congruence closure function (c). The relation R expressed by the dashed numbered lines (formally R = {(x, u), (y + z, u)}) is neither a bisimulation nor a bisimulation up to e since y + z a → x + y and u a → u, but (x + y, u) /

Figure 5 :

 5 Figure 5: Family of examples where HKC exponentially improves over AC and HK; we have x + y ∼ z.

  / / • • • a,b / / mn z

Definition 8 (

 8 Similarity). Similarity is the largest relation on states ⊆ S 2 such that x y entails: 1. o(x) ≤ o(y) and 2. for all a ∈ A, x ∈ S such that x a → x , there exists some y such that y a → y and x y .

Lemma 5 .

 5 For all relations R, c (R) = c(R ∪ ).

Theorem 4 .Corollary 3 .

 43 Any bisimulation up to c is contained in a bisimulation. For all sets X, Y , X ∼ Y iff HKC'(X, Y ).

Figure 6 :

 6 Figure 6: Relationships among the algorithms.

Figure 7 :

 7 Figure 7: Distributions of the number of processed pairs, for 1000 experiments with random NFA

  .860 4.981 5.084 0.009 0.094 1.412 2.887 HKC 0.049 0.798 6.494 6.762 0.000 0.014 0.916 2.685 sim time 0.039 0.185 0.574 0.618 0.038 0.193 0.577 0.593 AC' -sim time 0.013 0.167 1.326 1.480 0.012 0.107 1.047 1.134 HKC' -sim time 0.000 0.034 0.224 0.345 0.001 0.005 0.025 0.383

	algorithm	inclusions (546 pairs) 50% 90% 99% 100%	counter-examples (518 pairs) 50% 90% 99% 100%
	AC 0.036 0Table


Ngoc Q K Duong

Franck Thudor
email: franck.thudor@technicolor.com

MOVIE SYNCHRONIZATION BY AUDIO LANDMARK MATCHING

Keywords: Movie synchronization, landmark matching, scene addition, removal, and re-ordering

This paper addresses movie synchronization, i.e. synchronizing multiple versions of the same movie, with an objective of automatically transferring metadata available on a reference version to other ones. We first exploit audio tracks associated with two different versions and adapt an existing audio fingerprinting technique to find all temporal matching positions between them. We then propose additional steps to refine the match and eliminate outliers. The proposed approach can efficiently handle situations where temporal scene edits occur like scene addition, removal, and even the challenging scene re-ordering case. Experimental results over synthetic editorial data show the effectiveness of the proposed approach with respect to the state-of-the-art dynamic time warping (DTW) based solution.

INTRODUCTION

Nowadays, different versions of a given video content may co-exist and be accessed by end-users. A typical example is successive DVD versions of a blockbuster that can be found a couple of years after the theatrical one in extended version or in director's cut version. Another example is old movies that are brought up to date with new additional visual effects or in a colorized version. Due to local censure, some cleaned up versions may also appear in which violent, religious, sexual or political scenes are removed. Temporal edits that can occur between those versions include scene addition or removal and scene re-ordering. Thus, there is a real need for movie synchronization with an objective of automatically transferring some metadata available on the reference version to the second version where those metadata are absent. Such metadata may come from an artistic work, e.g. subtitles or chapters, but they may also be generated through a computational analysis of the audio-visual content itself, e.g. which characters are present, indoor/outdoor scene, etc. In both cases, transferring directly the metadata from one version to another avoids the time consuming task of metadata re-generation.

In this paper, we propose a synchronization algorithm grounded on a successful audio fingerprint technique [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF] with significant extensions in order to address the considered scene addition, removal and re-ordering situations. The developed method is robust to distortions of the audio stream, e.g. under compression or format conversion, and computationally very efficient thank to the advantage of the fingerprinting technique.

The rest of the paper is organized as follows. We first discuss the related work in Section 2. We then present the proposed synchronization approach in Section 3 where, for clarity, each step in the global workflow is organized into a subsection. We illustrate the experimental results to validate the performance of the proposed algorithm in Section 4. Finally, we conclude in Section 5.

RELATED WORK

For such a content synchronization purpose, either image descriptors, represented in the video frames, or audio features, represented in the audio samples, are usually exploited for matching [START_REF] Shrestha | Synchronization of multiple camera videos using audio-visual features[END_REF]. As examples, Thudor et al. addressed the video frame alignment for automatic chaptering of VoD content based on video feature matching in [START_REF] Thudor | Automatic chaptering of VoD content based on DVD content[END_REF] while Bryan et al. [START_REF] Bryan | Clustering and synchronizing multi-camera video via landmark cross-correlation[END_REF] exploited the audio landmark feature for clustering and synchronizing multi-camera video. Concerning audio based approaches, audio fingerprinting is naturally considered due to its compactness, its discriminative power, and its robustness to various kinds of distortion [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF][START_REF] Cano | A review of audio fingerprinting[END_REF]. Originally, fingerprint techniques have been developed for content identification purpose, e.g. Shazam's music identification system. Recently, they have also been exploited for other applications such as identifying repeated sound events in personal recordings [START_REF] Ogle | Fingerprinting to identify repeated sound events in long-duration personal audio recordings[END_REF], clustering and synchronizing multi-camera video [START_REF] Bryan | Clustering and synchronizing multi-camera video via landmark cross-correlation[END_REF], or retrieving large scale video for copy detection [START_REF] Jégou | Babaz: a large scale audio search system for video copy detection[END_REF]. Another state-of-the-art approach for synchronizing two related contents is Dynamic Time Warping (DTW). As examples, Macrae et al. used it for music video alignment [START_REF] Macrae | Muvisync: Realtime music video alignment[END_REF] and Anguera et al. used it for spoken word matching [START_REF] Anguera | Partial sequence matching using an unbounded dynamic time warping algorithm[END_REF]. This technique efficiently tackles the problem of scene addition and removal by finding an optimal path representing the highest similarity between two feature sequences extracted from those contents. However, DTW can not deal with the scene re-ordering situation due to the monotonicity condition. Moreover, its boundary condition requires prior knowledge of the start and end coordinates, which is not a trivial information, in order to compute an optimal path.

PROPOSED APPROACH

This section describes the details of the proposed algorithm starting from audio landmark extraction and matching steps inherited from Shazam's algorithm [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF]. Following steps are proposed for the considered movie synchronization application. The global workflow of the algorithm is shown in Fig. 1.

Landmark extraction

The audio landmark feature was originally presented in the well-known Shazam music recognition system [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF]. This feature is derived from the spectrogram, a time-frequency (T-F) representation of an audio signal by means of the short time Fourier transform, where only local energy peaks, which have higher energy than all their neighboring T-F points, are indexed by their T-F coordinates. A fixed number of peaks are chosen along each T-F region in order to assure a reasonably uniform distribution and each peak, known as anchor point, is paired with nearby ones within its target zone to form the landmarks. These landmarks are then quantized and packed into the L-bit unsigned integers and stored in the database together with the associated time offset from the beginning of the file for hashing.

Note that by choosing spectral peaks as anchor points, the derived feature is intrinsically robust to noise, and by using pairs of T-F points rather than a single one, the deriving landmark exploits additionally discriminative spectral structure of the audio signal. Although the Shazam implementation, corresponding to its original paper, has not been revealed, we use the Matlab implementation by Dan Ellis [START_REF] Ellis | Robust landmark-based audio fingerprinting[END_REF] for the development of our system.

Landmark matching and histogramming

At this step, each landmark extracted from the audio track associated with the second version of the movie is used to search in the database containing all landmarks extracted from the audio track associated with the first version to find the match. The detail of this process was described in [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF] where deriving hash representation offers a great improvement in search efficiency compared to a linear search. For each matching landmark found, the corresponding time offsets from the beginning of each audio file and their difference are stored. After all landmarks have been scanned, the scatterplot of time offsets corresponding to all matching landmarks is drawn together with the histogram representing the number of matches as a function of the time differences between the two versions. As an example, these pictures are shown by Fig. 2 and Fig. 3, respectively, for two versions of a 3 minute sequence of a movie.

What can be observed from these figures is that when a significant number of matching time pairs appears on a diagonal in Fig. 2, there is a corresponding histogram peak counting for the same difference of time offsets in Fig. 3.

Thresholding and match filtering

By scanning landmarks for all possible hits, accidental matches usually co-exist together with expected matches. The latter form a diagonal in the scatterplot of time offsets and a corresponding peak in the histogram. In order to filter out these accidental matches, local maximum peaks in the histogram are first specified by comparing the numbers of matching landmarks (values in y-axis in Fig. 3) to a predefined threshold. The threshold value is generally chosen depending on the landmark density, i.e. the approximate number of extracted landmarks per second, and the durations of the matching periods between two versions. Then the differences of time offsets between two versions corresponding to these maximum peaks (values in x-axis in Fig. 3) are used to reproduce the corresponding time offsets in two versions. Note that this reverse process is possible because a time offset was stored together with the time difference in the landmark extraction step. Finally a scatterplot accounting for only matching positions lying on considered diagonals is produced. This scaterplot, as shown in Fig. 4, is actually a filtered version of the one shown in Fig. 2 where accidental matches, which do not form a significant diagonal, are eliminated.

It is important to note that at this step Shazam [START_REF] Thudor | An industrial-strength audio search algorithm[END_REF] searches for only one maximum peak in the histogram to declare if two signals are similar or not. In our proposed approach, several maximum peaks are normally identified due to the multiple edits along the second version of the movie.

Outlier elimination

The resulting scatterplot from the previous step may contain outliers, i.e. points accidently lying on a diagonal but do not represent the actual temporal match between two versions. These outliers are shown, as an example, by points O 1 , O 2 , ..., O 8 in Fig. 4. Thus, at this stage, they need to be eliminated in order to specify the correct boundary of each resulting temporal matching segment. Additionally, this step will help to specify the valid segments, if there is more than one, in each diagonal.

We propose to use a hierarchical clustering-based approach [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] for this non-trivial work where each point in a diagonal of the filtered scatterplot is first considered as a cluster containing a single item. Euclidean distance between each pair of clusters is computed and the clusters for which their distance is smaller than a pre-defined threshold are merged. This "bottom up" process is repeated until either the distance between any pair of clusters is larger than the threshold or only one cluster remains. Note that, in our implementation we define the distance between clusters as the minimum distance between their boundary points, i.e. two points having the lowest and the highest time offsets in each cluster. Finally, clusters with a small number of points are considered to be outliers and are eliminated. The remaining clusters for every diagonal, e.g. as visualized by Fig. 5, will represent the estimated temporal matching between two versions.

In the example of Fig. 5, fives consecutive matching segments A, B, C, D, and E are identified. The gap between A and B along the y-axis indicates that some scenes do not exist in version 1, but are present in version 2. In the same manner, the gaps between B and D, and D and E indicate that some scenes are added to version 2 and removed from version 1, respectively. However these scenes are considered to be matched by segment C meaning that the scenes in version 1 appear at different time offset in version 2. This situation corresponds to the scene re-ordering edits.

EXPERIMENTAL RESULTS

In order to evaluate the synchronization performance, due to the lack of groundtruth for real edited versions of the movies, we generated a synthetic dataset from 11 movies on DVD by first extracting an audio track associated with each of them. Then for each movie sound track, we selected 6 sequences of approximately 180 second duration with different genres, i.e. action, speech, and music, resulting in a total of 66 reference audio files. For each of these files, we performed 4 different edits: removing a 5, 10, and 15 second slot in the middle of the file, and re-ordering a 10 second slot from a position at a quarter of the file length to a position at three quarters of the file length, resulting in a total of 264 edited audio versions.

Since DTW is a well-known technique for sequence alignment, we use it as a baseline for the comparison with our proposed approach. Our DTW implementation took a vector of 12 Mel-Frequency Cepstral Coefficients (MFCCs), the wellknown feature widely used in speech recognition [START_REF] Logan | Mel frequency cepstral coefficients for music modeling[END_REF], as an audio descriptor. Each MFCC sample is obtained over a 40 ms window with an overlap of 50% between two consecutive ones. In the implementation of the proposed algorithm, we extracted about 30 landmarks per second for the matching and the threshold for the number of matching landmarks described in Fig. 3 was set to 50. At the clustering step, clusters with less than 5 matching points are considered to be outliers and are eliminated.

The synchronization performance is evaluated in terms of precision, i.e. the fraction of detected synchronization parts that are correct, recall, i.e. the fraction of correct synchronization parts that are detected, and F-measure, i.e. the harmonic mean of the precision and recall. We observed that the proposed algorithm offers similar precision for all three edits of scene removal while the DTW algorithm results in lower performance with respect to larger scene removal, i.e. its precision is 0.97, 0.94, and 0.92 for 5 s, 10 s, and 15 s ed-its, respectively. This is because DTW continuity condition brings more false positive at the gap between two diagonals. At another point, both algorithms offer similar result on action, speech, and music data, meaning that audio genres do not significantly affect the synchronization performance. The overall performance, averaged for all genres and all scene removal edits, is shown in Table 1. It is not surprising that for the scene removal edits the proposed algorithm offers higher precision, but lower recall than that of the DTW algorithm. This can be explained by the fact that landmark matching offers very low false positive, but it may results in high false negative when landmarks are not extracted for a certain period. On the contrary, DTW finds the matching path frame by frame so it does not bring significant false negative but more false positive usually appears at discontinuity edges. Overall, F-measure are comparable for two methods. For the scene re-ordering edit, the proposed approach results in the average precision, recall, and F-measure of 0.95, 0.90, and 0.92, respectively, while DTW algorithm is non-applicable due to the monotonicity condition.

Edits

CONCLUSION

In this paper, we presented an efficient method for frame synchronization of different versions of a media content for a practical application of automatically transferring metadata from one version to the others. The proposed approach relies on audio tracks associated with the videos such that an existing audio landmark matching is successfully applied with substantial extension. Our experimental results indicate that a high level of synchronization accuracy can be achieved when considering practical edits such as scene addition, removal and even re-ordering. Future work will validate the performance of the proposed algorithms over real-world data.

Fig. 1 .

 1 Fig. 1. Global workflow of the proposed algorithm.

Fig. 2 .

 2 Fig. 2. Scatterplot of all matching landmarks.

Fig. 3 .

 3 Fig. 3. Histogram of the difference of time offsets.

Fig. 4 .

 4 Fig. 4. Scatterplot of filtered matching landmarks.

Fig. 5 .

 5 Fig. 5. Scatterplot of final matching landmarks after outlier elimination.

Table 1 .

 1 Average synchronization performance.

		Algorithm Precision Recall F-measure
	Scene	Proposed	0.98	0.96	0.97
	removal	DTW	0.94	1.0	0.97
	Scene	Proposed	0.95	0.90	0.92
	re-ordering	DTW		non applicable

ACKNOWLEDGMENT

The authors would like to thank Ashish Bagri for sharing the code of DTW algorithm.


Jonas B Ruh

Submarine Landslides caused by Seamounts entering Accretionary Wedge Systems

Seamounts entering active subduction zone trenches initially collide with the frontal sedimentary accretionary wedges resulting in severe deformation of the overriding plate. A typical feature of this deformation is the occurrence of submarine landslides due to gravitational instabilities. Such landslides have been reported from the Middle America and the Hikurangi trench and potentially generate tsunami waves. Yet, the dynamics of accretionary wedges during seamount indentation and landsliding as a mechanical response in particular have not been investigated quantitatively. Here, I apply 3D high-resolution numerical experiments to show that the topographical evolution of an accretionary wedge is mainly depending on the volume of the entering seamount and not on its height. Submarine landslides only occur if seamounts are not completely buried by the sedimentary sequence, where the volume of avalanches can be roughly correlated with the seamount volume overtopping the incoming sediments.

INTRODUCTION

Seamounts, or submarine volcanoes, are distributed among all oceanic plates around the globe [START_REF] Hillier | Global distribution of seamounts from ship-track bathymetry data[END_REF][START_REF] Smith | Seamount Statistics in the Pacific-Ocean[END_REF]. Although most of those are rather small, it has been suggested that around 12'000 large seamounts (> 1.5 km) cover the ocean floors [START_REF] Watts | Seamount Subduction and Earthquakes[END_REF][START_REF] Wessel | The Global Seamount Census: Oceanography[END_REF]. Seamounts differ in height and may exhibit either conical shapes (peaked), multiple branches (star shaped), and some can be flat-topped [START_REF] Schmidt | From seamount to oceanic island, Porto Santo, central East-Atlantic[END_REF]. Along active convergent margins, seamounts enter subduction zones and eventually collide with overriding plates. Observations of seamount subduction at active margins are documented from the Tonga-Kermandec trench [START_REF] Timm | Louisville seamount subduction and its implication on mantle flow beneath the central Tonga-Kermadec arc[END_REF], Nankai trench [START_REF] Bangs | Seamount subduction erosion in the Nankai Trough and its potential impact on the seismogenic zone[END_REF], Japan and Kuril trench [START_REF] Lallemand | Subduction of the Daiichi-Kashima seamount in the Japan Trench[END_REF][START_REF] Nishizawa | Seismic structure of the subducting seamounts on the trench axis: Erimo Seamount and Daiichi-Kashima Seamount, northern and southern ends of the Japan Trench[END_REF], Cascadia [START_REF] Wells | Fore-arc migration in Cascadia and its neotectonic significance[END_REF], New Hebrides (Collot and Fisher, 1994), Hikurangi trench [START_REF] Pedley | Seafloor structural geomorphic evolution of the accretionary frontal wedge in response to seamount subduction, Poverty Indentation[END_REF], Middle America trench [START_REF] Von Huene | Geophysics -When seamounts subduct[END_REF], Central Chile margin (Laursen et al., 2009), Central Mariana and Izu-Bonin trenches [START_REF] Oakley | Pacific Plate subduction beneath the central Mariana and Izu-Bonin fore arcs: New insights from an old margin[END_REF], Mediterranean ridge [START_REF] Von Huene | A subducting seamount beneath the Mediterranean Ridge[END_REF], and Aleutian trench [START_REF] Das | Effect of subducting seafloor topography on the rupture characteristics of great subduction zone earthquakes[END_REF].

Relatively large seamounts are truncated when reaching the fore arc at shallow levels and may potentially be decapitated at larger depths, related to seismogenic deformation [START_REF] Cloos | Thrust-Type Subduction-Zone Earthquakes and Seamount Asperities -a Physical Model for Seismic Rupture[END_REF]. But whether seamounts along subduction interfaces trigger large earthquakes, form natural barriers for earthquake propagation, or even decrease seismic coupling between two plates is still a matter of debate (Mochizuki et al., 1997;Scholtz and Small, 1997;[START_REF] Das | Effect of subducting seafloor topography on the rupture characteristics of great subduction zone earthquakes[END_REF][START_REF] Watts | Seamount Subduction and Earthquakes[END_REF][START_REF] Wang | Do subducting seamounts generate or stop large earthquakes?[END_REF].

Before a seamount reaches the fore arc, it traverses the frontal sedimentary accretionary wedge causing severe impact on the morphologic and tectonic evolution of the overriding plate. Topographic uplift above seamounts formes bulges [START_REF] Park | A subducting seamount beneath the Nankai accretionary prism off Shikoku, southwestern Japan[END_REF]. Slope failure of such supercritical bulges [START_REF] Lallemand | Coulomb wedge model applied to the subduction of seamounts in the Japan Trench[END_REF] occur that can trigger submarine avalanches [START_REF] Lewis | The dammed Hikurangi Trough: a channel-fed trench blocked by subducting seamounts and their wake avalanches[END_REF][START_REF] Pedley | Seafloor structural geomorphic evolution of the accretionary frontal wedge in response to seamount subduction, Poverty Indentation[END_REF]. Such submarine avalanches (or landslides) can in turn cause major human and economic impact in the form of tsunamis [START_REF] Fisher | Geology and tsunamigenic potential of submarine landslides in Santa Barbara Channel[END_REF][START_REF] Tappin | Did a submarine landslide contribute to the 2011 Tohoku tsunami?[END_REF], or demolition of offshore pipelines [START_REF] Liu | Impact forces of submarine landslides on offshore pipelines[END_REF].

Collision of seamounts into accretionary systems have been dynamically reconstructed with analogue models that could reproduce main features of observed upper plate deformation [START_REF] Lallemand | Effects of oceanic ridge subduction on accretionary wedges: Experimental modeling and marine observation[END_REF][START_REF] Dominguez | Upper plate deformation associated with seamount subduction[END_REF][START_REF] Dominguez | Deformation of accretionary wedges in response to seamount subduction: Insights from sandbox experiments[END_REF]. However, there is no published numerical modelling study so far investigating the effect of seamount subduction on deformation localization within the upper plate in general and the occurrence of gravitational submarine landslides in particular.

Due to the prevalent distribution of seamounts along subduction trenches and their geohazard potential, it is essential to understand the mechanics and dynamics of developing accretionary wedges during collision with a seamount. Therefore, I here test the influence of seamount size and shape on the tectonic and mechnical evolution of accretionary wedges with main focus on gravitational collapse and submarine landslides by applying a high-resolution three-dimensional numerical code. Results are discussed in the scope of the analytical critical wedge theory.

APPLIED GEOMETRICAL SETUP

Numerical model and setup are similar to those documented in earlier studies [START_REF] Ruh | High-resolution 3D numerical modelling of thrust wedges: Influence of decollement strength on transfer zones[END_REF]2014), and are based on a three-dimensional, high-resolution, fully-staggered grid, finite difference, marker-in-cell code with a standard visco-brittle/plastic rheology and an efficient OpenMP-parallelized multigrid solver (I3ELVIS; [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF][START_REF] Gerya | Robust characteristics method for modelling multiphase visco-elastoplastic thermo-mechanical problems[END_REF]; see supplementary material). Eulerian grid dimensions of presented simulations are 97.6 • 97.6 • 14.8 km in x-, y-, z-directions with a nodal resolution of 245 • 245 • 149, respectively; resulting cell size is 400/400/100 m, each containing 8 Lagrangian markers. Accordingly, the code calculates for ~9 million nodes and ~70 million markers. The geometrical setup resembles typical "sandbox"-type models. Initial marker distribution defines a 300 m thick lowermost rigid plate, topped by a 500 m thick décollement and 3 km of sedimentary sequence.

Above, 11 km of "sticky-air" mimic a free surface allowing wedge to thicken. At the bottom, the rigid plate moves with a velocity v x = -1 cm/yr and is pulled out below the lateral boundary at x = 0, which acts as a rigid backstop (no slip). On the opposite boundary, new stratigraphy is entering the Eulerian grid with v x . Lateral side and top boundaries are free slip. Seamounts root into the rigid bottom plate and are introduced opposite of the backstop into undeformed stratigraphy after 4 Myr runtime.

Plastic/brittle failure is based on the Drucker-Prager formulation depending on pressure P, friction an-gle φ, cohesion C, and fluid pressure ratio λ [START_REF] Ruh | 3D effects of strain vs. velocity weakening on deformation patterns in accretionary wedges[END_REF]. Décollement: φ b = 10° and C b = 0.2 MPa. Sedimentary sequence: φ = 30° and C = 10 MPa. Seamounts: φ = 30° and C = 20 MPa. A hydrostatic fluid pressure ratio λ = 0.4 is applied for décollement and sediments. Sediments and seamounts are weakened by linearly lowering φ to 20° and C to 0.2 MPa according to plastic strain ε pl between 0.5 < ε pl < 1.5 [START_REF] Ruh | High-resolution 3D numerical modelling of thrust wedges: Influence of decollement strength on transfer zones[END_REF]. A total of nine simulations are presented here: One model without seamount implementation acts as reference, eight simulations were carried out to investigate effects of size and shape of seamounts on upper plate deformation. Modeled seamounts rise 1.5, 3, 4.5, or 6 km high (above the décollement level) with either conical or flat-topped shapes. Seamounts are round with a lateral slope of 30°, which is on the upper level of naturally observed seamount slopes (Figure 1 in [START_REF] Dominguez | Upper plate deformation associated with seamount subduction[END_REF]. Flat-topped seamounts exhibit a horizontal cap with a radius of 5 km.

TOPOGRAPHIC EVOLUTION DURING SEAMOUNT SUBDUCTION

Deformation of numerical accretionary wedges localizes along the backstop and the deformation front migrates outward in-sequence with ongoing shortening (Gif 1 in supplementary material). Frontal accretion and out-of-sequence thrusting assure wedge growth in width (x-direction) and thickness, respectively, leading to a wedge body thickened at the rear (x = 0) thinning towards the toe (Fig. 1a). The topographic evolution of compressive wedges varies strongly depending on size and shape of entering seamounts. A conical seamount with a height of 1.5 km (above the décollement level) leaves a trace only in the thrust sheet that was the frontal one at the time of seamount entrance into the wedge (Fig. 1b). Simulations with higher conical seamounts exhibit topographic peaks spatially related to the seamount location and several thrust sheets are deformed due to its entrance (Fig. 1c-e). The topographic response of accretionary wedges to flat-topped seamount subduction is more prominent. A 1.5 km high semount leaves a clear trace in the accreting thrust sheets (Fig. 1f). Higher flat-topped seamounts form wide entrances into the wedge front and there are no cylindrical frontal thrust sheets developing within 9 Myr (Fig. 1g-i). Topographic peaks induced by underthrusting seamounts are always located at the rearward side of the seamount top (Fig. 1).

The elevation of these topographic peaks increases with increasing seamount height.

Nonetheless, peak elevations of flat-topped seamount experiments are substantially higher than those with equally high conical seamounts. After 9 Myr, experiments of flat-topped seamounts with 1.5 and 3 km height (Fig. 1f,g) exhibit similar surface topographic reliefs as models with 4.5 and 6 km high conical seamounts (Fig. 1d,e), respectively: Only from the topographic map is not possible to determine whether the underthrusting seamount is flat-topped with 3 km height or conical with 6 km height (Fig. 1e,g).

Yet, peak elevations strongly dependent on the volume of seamounts entering an accretionary wedge (Fig. 2). This relation between the volume of a seamount and the peak elevation of overthrusting imbricate sheets is independent of the temporal evolution of seamount collision into an accretionary wedge. The similar differences of topographic peaks between 7 and 8 Myr and between 8 and 9 Myr, respectively, indicate a linear vertical thickening of accretonary wedges over time, independent of shape and height of subducting seamounts (Fig. 2).

SUBMARINE LANDSLIDES TRIGGERED BY SEAMOUNTS

It is widely accepted that seamounts entering subduction zones are a potential trigger of submarine landslides due to gravitational collapses (e.g., [START_REF] Lallemand | Coulomb wedge model applied to the subduction of seamounts in the Japan Trench[END_REF]Hühnerbach et al., 2005). Nevertheless, the dynamics of such avalanches and their volumetric and temporal relation to subducting seamounts have yet to be investigated.

The collision of a seamount has a major impact on the mechanical evolution of accretionary wedges, regardless whether submarine landslides appear or not (Gifs 2-9 in supplementary material).

Seamounts act as barriers for the migration of the deformation front as they interrupt the flat décollement. The gravitational potential of sedimentary strata overthrusting seamounts activates the décollement level on the toeward side of the seamount and a normal fault along the seamount develops (Gifs 2-9 in supplementary material). Experiments exhibiting seamounts that exceed the initial sedimentary sequence build up large enough gravitational forces within the overthrust stratigraphy to trigger gravitational collapses (Gifs 4,5,8,9 in supplementary material). In case of a 6 km high flattopped seamount, stacked material atop of the seamount fails along a listric-like normal fault and emplaces an up to 1.5 km thick sequence in front of the seamount (Fig. 3a). Strain rates indicate that the emplacement of the landslide activates the décollement toewards from the seamount and puts the undeformed sedimentary sequence at the verge of failure (Fig. 3b). A steep main scarp, a listric fault, a thin rupture surface, and a partly undeformed landslide mass are very well constrained features of observed submarine landslides [START_REF] Hampton | Submarine landslides[END_REF]. Velocity field in x-direction indicates that the material above the active listric fault moves away from the rear (Fig. 3c; such slow velocities result from lower viscosity cutoff: 10 18 Pa•s). Low velocities of displaced material patches along the wedge toe indicate an earlier landslide.

A topographic peak above the seamount and surpress frontal accretion due to the seamount interrupting the décollement level (Fig. 3) lead to a steepening of the surface taper along the profile where the seamount subduct (y = 48.8 km), as observed in seismic profiles through the Nankai accretionary wedge [START_REF] Park | A subducting seamount beneath the Nankai accretionary prism off Shikoku, southwestern Japan[END_REF]. The temporal evolution of surface slope of a compressional wedge without seamount implementation (Gif 1 in supplementary material) indicates a taper leveling along the analytical minimum critical value within the stable regime (Fig. 4a). Dynamic surface tapers of seamount experiments are calculated from surface points located between the seamount top and the wedge toe, where a continuous décollement layer exists in between [START_REF] Davis | Mechanics of Fold-and-Thrust Belts and Accretionary Wedges[END_REF].

Submarine landslides decrease surface slopes (Fig. 4a,b). Tapers move from the supercritical into the stable field (below the maximum critical value) where they remain at 6-8°. This demonstrates that modeled wedges are behaving critical and that landslides are triggered by brittle/plastic gravitational collapses. Slope profiles along subduction zones where seamounts enter can therefore give insight in the strength relationship between basal and wedge internal material, expecting that tapers away from seamount occurrences represent the minimum analytical value, whereas profiles along seamount scars are close to maximum critical tapers [START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF].

SIZE OF SUBMARINE LANDSLIDES

The volume of the dislocated mass during submarine landslides triggered by seamount subduction differs strongly for conical and flat-topped seamounts (Fig. 4c,d). A 4.5 km high conical seamount triggers a short-living landslide with a volume of ~20 km 3 (Fig. 4c). A 6 km conical seamount launches an up to 60 km 3 large avalanche that is active over a longer period of time. Initial activation of landslides are contemporaneous with the strong decrease of surface taper between 7 and 8 Myr of experiment evolution (Fig. 4a,c). Flat-topped seamounts cause much larger landslides. Peak volumes reach ~300 km 3 for 4.5 km high and up to 500 km 3 for 6 km high seamounts (Fig. 4d). The landslide remains active from its initiation at 7.7 Myr to the leveling of surface taper within the stable area at 9 Myr (Fig. 4b,d). On the basis of these results, a first-order estimation of landslide volume can be inferred from the volume of entering seamounts overtopping the initial sedimentary sequence (Fig. 4c,d), which also agrees with the absence of landslides for seamounts shorter than the stratigraphic thickness of 3km.

Numerical submarine landslides presented here can be compared to reported natural cases and estimations for their tsunamigenic potential can furthermore be provided. In general, numerical avalaches are volumetrically comparable to ordinary submarine landslides triggered by gravitational instabilities (<800 km 3 ; [START_REF] Hampton | Submarine landslides[END_REF]and references therein). An example of landsliding caused by seamount collision is documented along the Hikurangi margin [START_REF] Lewis | Collapse in a Quarternary shelf basin off East Cape, New Zealand: evidence for passage of a subducted seamount inboard of the Ruatoria giant avalanche[END_REF]. The Giant Ruartoria avalanche off New Zealand has a volume of ~2000 km 3 and is a result of collision of the now subducted Ruatoria seamount into the Hikurangi accretionary wedge [START_REF] Collot | The giant Ruatoria debris avalanche on the northern Hikurangi margin, New Zealand: Result of oblique seamount subduction[END_REF]. The size of Ruatoria seamount is interpreted to be similar to that of Gisborne seamount offshore Hikurangi, whose overtop volume can roughly be compared to the reported avalanche volume.

Figure 5 illustrates a comparison between the Ruatoria indentation and related avalanche and the numerical collision of a flat-topped 6 km high seamount into a developing accretionary wedge presented in this study (Fig. 3). In both cases, landsliding results from oversteepening and collapse of the wedge front due to seamount indentation (Fig. 5). In the natural as well as the numerical case, volumes of avalanches and exceeding seamounts are roughly similar.

A potential of tsunami generation is apparent from a fault scar volume larger than ~10 km 3 (van Huene et al., 2004), which includes all presented experiments with a seamount height greater than the initial sedimentary sequence (Fig. 4c,d). Yet, also smaller volumes of landslides have been interpreted to be potentially tsunamigenic [START_REF] Borrero | Tsunamis within the eastern Santa Barbara Channel[END_REF][START_REF] Fisher | Geology and tsunamigenic potential of submarine landslides in Santa Barbara Channel[END_REF]. Mass movements similar to here reported can generate up to 50 m high tsunamis [START_REF] Von Huene | A large tsunamigenic landslide and debris flow along the Peru trench[END_REF]. Trench areas where rough and seamount covered oceanic crust is entering a subduction zone are in a critical stress state and the interplay between seismic events, submarine landslides, and tsunami waves was reported for the 1992 Central Nicaragua earthquake [START_REF] Von Huene | Tsunamigenic slope failure along the Middle America Trench in two tectonic settings[END_REF]. Reentrants and scars indicate seamount collision and taper angles with a high potential of gravitational failure during earthquakes. Experiments presented here show that potential landslide volumes can be directly correlated to the (overtop) volume of seamount.

FIGURE CAPTIONS

Figure 1

 1 Figure 1 Surface topography of 3D accretionary wedges above initial sediment thickness in map view after 9 Myr experiment time. Backstop is located at x = 0 on the upper side of the plates. (a) No seamount implementation. (b-e) Conical seamount implementation. Dots indicate location of seamount peaks. (f-i) Flat-topped seamount implementation. Circles indicate location of flat area of seamount tops. Height of seamounts is indicated on top for all columns.

Figure 2

 2 Figure 2 Surface topography peaks related to seamount subduction (Fig. 1) are plotted against the complete volume of according implemented seamounts for an experiment time of 7, 8 and 9 Myr. Circles and squares indicate conical and flat-topped seamounts, respectively. Seamount heights are 1.5, 3, 4.5, 6 km above décollement level, increasing towards right according to their volume.

Figure 3

 3 Figure 3 (a) Compositional map, (b) strain rate, and (c) velocity in x-direction of an accretionary wedge indented by a 6 km high flat-topped seamount after 9 Myr shortening. Wedges are cut parallel to x-direction at y = 44.8 km. Seamounts are displayed by their composition.

Figure 4

 4 Figure 4 Dynamics of numerical landslides. (a,b) Temporal surface taper evolution of reference and landslide-generating experiments. Dotted line: no seamount implementation. Grey line: 4.5 km high seamount. Black line: 6 km high seamount. Dashed lines: Minimum and maximum critical taper angles (Davis et al., 1983). (c,d) Volume of landslides triggered by conical and flat-topped seamounts, respectively. Filled areas indicate volumes between 0.5 cm/yr < v x < 1 cm/yr moving against shortening direction. Dashed lines indicate volume of seamounts overtopping the initial sedimentary stratigraphic thickness. Grey: 4.5 km high seamount. Black: 6 km high seamount.

Figure 5

 5 Figure 5 Comparison between natural and numerical cases. (a) Top view of an experiment with a 6 km high flat-topped seamount after 9 Myr. Grey scale background: surface topography. Red lines: numerical faults. Blue tones: different landslides. Rose: landslide detachment. Blue and rose stripes: detachment burried by landslide. (b) Profiles in x-direction at Ly = 0.5 of the experiment presented in (a) before and after landslide occurrence at 7.4 and 9 Myr. Lower panel: surface lines before (red) and after (green) landsliding corrected for seamount passage (measured surface lines minus height of seamounts exceeding the initial sedimentary strata). Colors as in (a). (c) Simplified map of the Ruaroria indentation at the northeastern tip of Newzealand (after[START_REF] Collot | The giant Ruatoria debris avalanche on the northern Hikurangi margin, New Zealand: Result of oblique seamount subduction[END_REF][START_REF] Lewis | Collapse in a Quarternary shelf basin off East Cape, New Zealand: evidence for passage of a subducted seamount inboard of the Ruatoria giant avalanche[END_REF].Dark grey: upper margin. Light grey: imbricated zone. Dotted area: debris flow and disturbed trench fill.Lined area: post-avalanche turbidites. Red lines: main faults. Broken red line: pre-avalanche deformation front. Arrow: direction of plate motion and location of seamount entrance. Blue and rose as in (a). (d) Profiles through the Ruatoria indentation (after[START_REF] Collot | The giant Ruatoria debris avalanche on the northern Hikurangi margin, New Zealand: Result of oblique seamount subduction[END_REF]. Profile A-A' indicates

ACKNOWLEDGMENTS

I am very thankful to Taras Gerya, who provided the numerical code. Funding was received from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant agreement n° 604713 (ZIP "Zooming In between Plates").