Pierre Peterlongo
email: pierre.peterlongo@irisa.fr

Nadia Pisanti
email: pisanti@di.unipi.it

Frédéric Boyer
email: frederic.boyer@inrialpes.fr

Alair Pereira Do Lago

Marie-France Sagot
email: .sagot@inria.fr

Lossless filter for multiple repetitions with Hamming distance

Keywords: approximate repetitions, k-factors, multiple local alignment, bi-factors, bi-factor array

Similarity search in texts, notably in biological sequences, has received substantial attention in the last few years. Numerous filtration and indexing techniques have been created in order to speed up the solution of the problem. However, previous filters were made for speeding up pattern matching, or for finding repetitions between two strings or occurring twice in the same string. In this paper, we present an algorithm called Nimbus for filtering strings prior to finding repetitions occurring twice or more in a string, or in two or more strings. Nimbus uses gapped seeds that are indexed with a new data structure, called a bi-factor array, that is also presented in this paper. Experimental results show that the filter can be very efficient: preprocessing with Nimbus a data set where one wants to find functional elements using a multiple local alignment tool such as Glam, the overall execution time can be reduced from 7.5 hours to 2 minutes.

Introduction

Finding approximate repetitions (motifs) in strings is one of the most challenging tasks in text mining. Its relevance grew recently because of its application to biological sequences. Although several algorithms have been designed to address this task, and have been extensively used, the problem still deserves investigation for certain types of repetitions. Indeed, when the latter are long and the number of allowed differences among them grows proportionally to their length, there is no exact tool that can efficiently manage to detect such repetitions. Widely used efficient algorithms for multiple alignment are heuristic, and offer no guarantee that false negatives are avoided. On the other hand, exhaustive inference methods cannot handle queries where the differences allowed among the occurrences of a motif are as many as 5-10% of the length of the motif, and the latter is as "long" as, say, 100 DNA bases. Indeed, exhaustive inference is done by extending or assembling in all possible ways shorter motifs that satisfy certain conditions. When the number of differences allowed is relatively high, this can therefore result in too many false positives (of intermediate length) that saturate the memory. In this paper, we introduce a preprocessing filter, called Nimbus, where most of the data containing such false positives are discarded in order to perform a more efficient exhaustive inference. Our filter is designed for finding repetitions distant pairwise by a maximal Hamming distance, occurring in r ≥ 2 input strings, or occurring possibly more than twice in one string. To our knowledge, one finds in the literature filters for local alignment between two strings [START_REF] Rasmussen | Efficient q-Gram Filters for Finding All ǫ-Matches Over a Given Length[END_REF][START_REF] Ma | Patternhunter: Faster and More Sensitive Homology Search[END_REF][START_REF] Li | PatternHunter II: Highly Sensitive and Fast Homology Search[END_REF], or for approximate pattern matching [START_REF] Navarro | Indexing Text with Approximate q-grams[END_REF][START_REF] Burkhardt | q-Gram Based Database Searching Using a Suffix Array (QUASAR)[END_REF] only (that is a different task). Heuristic methods such as Blast [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF][START_REF] Altschul | Gapped BLAST and PSI-BLAST: a New Generation of Protein Database Search Programs[END_REF] and Fasta [START_REF] Lipman | Rapid and Sensitive Protein Similarity Searches[END_REF] filter the input data and extend only seeds that are repeated short fragments satisfying some constraints. Nimbus is based on similar ideas but uses different requirements concerning the seeds; among the requirements are frequency of occurrence of the seeds, concentration and relative position. Similarly to [START_REF] Ma | Patternhunter: Faster and More Sensitive Homology Search[END_REF][START_REF] Li | PatternHunter II: Highly Sensitive and Fast Homology Search[END_REF], we use also a concept related to gapped seeds that has been shown in [START_REF] Burkhardt | Better Filtering with Gapped q-Grams[END_REF] to be particularly efficient for pattern matching. The filter we designed is lossless: unlike Blast or Fasta, Nimbus guarantees not to discard any repetitions meeting the input parameters. It uses necessary conditions based on combinatorial properties of multiple and approximate repetitions, and an algorithm that checks such properties in an efficient way. The efficiency of the filter relies on an original data structure, the bi-factor array, that is also introduced in this paper, and on a labelling of the seeds similar to the one employed in [START_REF] Iliopoulos | A First Approach to Finding Common Motifs With Gaps[END_REF]. This new data structure can be used to speedup other tasks such as the inference of structured motifs [START_REF] Marsan | Algorithms for Extracting Structured Motifs Using a Suffix Tree With Application to Promoter and Regulatory Site Consensus Identification[END_REF] or for improving other filters [START_REF] Kucherov | Multiseed Lossless Filtration[END_REF].

It is worth mentioning that Nimbus is meant as a preprocessing to any tool that searches for long motifs with mismatches shared by several strings, or that performs multiple alignments based on the detection of such repetitions. The goal of Nimbus is to reduce the input size for these tools in order to obtain an overall significant speed up. For example, we shall see in Section 6 that preprocessing a data set where one wants to find a multiple local alignment with a tool such as Glam ( [START_REF] Frith | Finding Functionnal Sequence Elements by Multiple Local Alignement[END_REF]), the overall execution time can be reduced from 7.5 hours to 2 minutes. The latter shortened time includes both the running time of the filtering of Nimbus, and that of the tool applied to the filtered strings. Notice that we chose Glam as a representatively good software for multiple alignment, but the filter could be used as preprocessing step for any other multiple alignment algorithm.

A preliminary version of this paper appeared in [START_REF] Peterlongo | Lossless Filter for Finding Long Multiple Approximate Repetitions Using a New Data Structure, the Bifactor Array[END_REF].

Necessary Conditions for Long Repetitions

A string is a sequence of zero or more symbols from an alphabet Σ. A string s of length n on Σ is represented also by s[0]s [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF] . . . s[n -1], where s[i] ∈ Σ for 0 ≤ i < n. The length of s is denoted by |s|. We denote by s[i, j] the substring, or factor, s[i]s[i + 1] . . . s[j] of s. In this case, we say that the string s[i, j] occurs at position i in s. We call k-factor a factor of length k. If s = uv for u, v ∈ Σ * , we say that v is a suffix of s.

Definition 1 Given r input strings s 1 , . . . , s r , a length L, and a distance d, we call a (L, r, d)-repetition a set {δ 1 , . . . , δ r } such that 0 ≤ δ i ≤ |s i | -L, and for all i, j ∈ [1, r] we have that

d H (s i [δ i , δ i + L -1], s j [δ j , δ j + L -1]) ≤ d.
where by d H we mean the Hamming distance between two strings, that is, the minimum number of letter substitutions that transform one into the other.

Given m input strings, the goal is to find the substrings of length L that are repeated in at least r ≤ m strings with at most d substitutions between each pair of the r repetitions, with L, r and d given. In other words, we want to extract all the (L, r, d)-repetitions from a set of r strings among m ≥ r input strings. The goal of the filter is therefore to eliminate from the input strings as many positions as possible that cannot contain (L, r, d)-repetitions. Since in most interesting applications, the value of the parameter d can be as big as 10% of L, a direct inference of repetitions of length L with L/10 substitutions based on a traditional algorithm would be unfeasible due to the big number of approximate short patterns that could contribute to a longer repetition. The main idea of our filter is based on checking necessary conditions on the number of exact k-factors that a (L, r, d)-repetition must share. Since there are only a linear number of exact factors and they are easy to find, this results in an efficient filtering. A string w of length k is called a shared k-factor for s 1 , . . . , s r if for all i ∈ [1, r] we have that w occurs in s i . Obviously, we are interested in shared k-factors that occur within substrings of length L of the input strings. Let p r be the minimum number of non-overlapping shared k-factors that a (L, r, d)-repetition must have. It is intuitive to see that a (L, 2, d)-repetition Before stating the lemma, we give an intuition for the key observations that sustain the proof. Let us consider the (virtual) alignment of r strings of length L corresponding to an (L, r, d)-repetition. If they were identical, they would clearly share L k non-overlapping k-factors. Let us now assume that the Hamming distance between any two of the r strings is d. This is the worst case scenario for the minimum number of shared k-factors when all allowed substitutions are present, which is what we need to compute. If at position i (with 0 ≤ i ≤ L -1), there is a substitution between any pair of strings among the r, then we say that position i is stained. In this case, no k-factor involving this position can be shared by all the r strings. The worst case scenario applies when the stained positions (or stains for short) are spread in such a way that they eliminate distinct (i.e. non overlapping) k-factors, that is we loose as many potential L k non-overlapping k-factors as the number of stained positions in the alignment. We therefore need only to count the maximum number of stains. If r = 2, this number is simply d, and indeed we have that

contains at least L k -d shared k-factors, that is, p 2 = L k -d. We now compute the value of p r for r > 2. ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✂ ✁✂ ✁✂ ✁✂ ✂ ✁✂ ✁✂ ✁✂ ✂ ✁✂ ✁✂ ✁✂ ✄ ✁✄ ✁✄ ✁✄ ✁✄ ✄ ✁✄ ✁✄ ✁✄ ✁✄ ✄ ✁✄ ✁✄ ✁✄ ✁✄ ☎ ✁☎ ✁☎ ✁☎ ☎ ✁☎ ✁☎ ✁☎ ☎ ✁☎ ✁☎ ✁☎ ✆ ✁✆ ✁✆ ✁✆ ✆ ✁✆ ✁✆ ✁✆ ✆ ✁✆ ✁✆ ✁✆ ✝ ✁✝ ✁✝ ✁✝ ✝ ✁✝ ✁✝ ✁✝ ✝ ✁✝ ✁✝ ✁✝ A A A A A G G G G G C C C C G A T T A A T T T T T C C A T G G G L = 16 d = 2 k =3
p 2 = L k -d.
This case is shown in Figure 1. The key observation is that adding a third string (then a fourth, fifth, and so on up to the r th ) creates at most d 2 new stained positions. This is surprising because one could think that a single "new" string that presents d substitutions with each one of, say, r -1 previously considered strings, can stain up to (r -1)d positions. Instead, it can stain only up to d 2 previously unstained positions independently from r (provided this is greater than 2). The reason is that, as we show below, many substitutions necessarily occur at positions that are already stained.

Let φ(r, d) denotes the maximal number of distinct positions that can be stained given that we have r strings such that any two of them have Hamming distance at most d. We have that r ≥ 2 and clearly φ(2, d) = d. In the following lemma we give tight numbers for φ(r, d).

Lemma 2 Let r ≥ 2 and d ≥ 0 be integers. We have that φ(r, d) = ⌈d/2⌉ + (r -1)⌊d/2⌋. In other words, φ(r, d)

= rd 2 if d is even and φ(r, d) = r d 2 + 1 if d is odd.
PROOF. First, consider the following set of r sequences, all of them with length ⌈d/2⌉ + (r -1)⌊d/2⌋. Let i be any integer such that 1 ≤ i ≤ r. All positions of a sequence s i are A, except those in the following positions, which are T: the first ⌈d/2⌉ letters of s 1 ; the next ⌊d/2⌋ letters of s 2 , instead of s 1 ; the next ⌊d/2⌋ letters of s 3 , instead of s 2 ; and so on until the last ⌊d/2⌋ letters of s r . More formally, s 1 [0] = T and, for every 1 ≤ i ≤ r and j > 0, we set s i [j] = T iff ⌈d/2⌉ + (i -2)⌊d/2⌋ ≤ j < ⌈d/2⌉ + (i -1)⌊d/2⌋; all other positions are set to A. One can verify that the Hamming distance from s 1 to any other sequence is d and the Hamming distance from any two sequences except s 1 is 2⌊d/2⌋. Since every position is a stain in these sequences, this proves that φ(r, d) ≥ ⌈d/2⌉ + (r -1)⌊d/2⌋ and that the provided bound is tight.

We now show that φ(r, d) ≤ ⌈d/2⌉ + (r -1)⌊d/2⌋. In order to prove it, we apply an induction on d. If, d = 0, φ(r, d) is clearly 0 and the proof is trivial. Since we are seeking an upper bound for φ(r, d), let us choose r sequences with the largest possible number of stains with the restriction that every pair of them has Hamming distance at most d. In order to evaluate the actual number φ(r, d) of stains, we have to consider two subcases. In the first case, no two of these sequences have Hamming distance exactly d, and hence we can state that φ(r, d) ≤ φ(r, d-1) ≤ ⌈(d-1)/2⌉+(r-1)⌊(d-1)/2⌋ ≤ ⌈d/2⌉+(r-1)⌊d/2⌋ due to the induction hypothesis. In the second case, at least one pair of sequences is at Hamming distance d. Without loss of generality, let us assume that -say -s 1 and s 2 are such that their Hamming distance is d. This leads to d stains at d such positions j where we have s 1 [j] = s 2 [j]. Let J be the set of such positions. If there is no other stain, the proof is finished since we have already proved that φ(r, d) = d = ⌈d/2⌉ + ⌊d/2⌋ ≤ ⌈d/2⌉ + (r -1)⌊d/2⌋. Suppose from now on that there is a stain j ′ ∈ J. Hence,

s 1 [j ′ ] = s 2 [j ′ ]. By definition of a stain, there are two integers i ′′ > i ′ such that s i ′′ [j ′ ] = s i ′ [j ′ ].
Hence, i ′′ > 2 and we can suppose that 3 ≤ i ′′ ≤ r. (This in particular implies that r ≥ 3.) If i ′ ≤ 2, we take i = i ′′ and we have that

s i ′ [j ′ ] = s 1 [j ′ ] = s 2 [j ′ ] differs from s i [j ′ ] = s i ′′ [j ′ ]. If i ′ > 2, we have that s 1 [j ′ ] = s 2 [j ′ ] differs either from s i ′′ [j ′ ] or from s i ′ [j ′ ] and we can take i ∈ {i ′′ , i ′ } such that s 1 [j ′ ] = s 2 [j ′ ] differs from s i [j ′ ]. Anyway, there is i ∈ {3, . . . , r} such that s 1 [j ′ ] = s 2 [j ′ ] differs from s i [j ′ ]. For all j ∈ J, we have that s 1 [j] = s 2 [j] and that s i [j] = s 1 [j] or s i [j] = s 2 [j].
Hence, there is i ′′′ ∈ {1, 2} such that s i differs from s i ′′′ in at least ⌈d/2⌉ positions in J and also differs in j ′ . This implies that s i cannot differ from s i ′′′ in more than d -⌈d/2⌉ = ⌊d/2⌋ positions not in J. For (r -2) possible values for i, we have (r -2)⌊d/2⌋ possible values for j ′ . Hence,

φ(r, d) ≤ d + (r -2)⌊d/2⌋ = ⌈d/2⌉ + (r -1)⌊d/2⌋.
The following theorem follows from Lemma 2.

Theorem 3 Let L > 0, r ≥ 2 and d > 0 be integers. The minimum number of non-overlapping shared k-factors that a (L, r, d)-repetition must have is

p r = L k -d -(r -2) × d 2 .
PROOF. Since there are ⌊L/k⌋ possible non-overlapping common factors in a (L, r, d)-repetition and that we have that any stain can hit at most one of them, then we have that p r = L k -φ(r, d), and hence the result follows.

Observe that the theorem applies also to the case where one is interested in finding (L, r, d)-repetitions occurring in a single string.

The Algorithm

We start by describing the algorithm for the case where we are looking for repetitions occurring in a set of strings. We call the corresponding algorithm Multi-Nimbus. Later, we describe how to modify the algorithm to find a (L, r, d)-repetition occurring in a single string thus obtaining an algorithm we call Mono-Nimbus.

Multi-Nimbus takes as input the parameters L, r and d, and m (with m ≥ r) input strings. Given such parameters, it decides automatically the best value of k to apply Theorem 3 and compute the number of k-factors that are necessarily shared by a (L, r, d)-repetition.

The goal of Multi-Nimbus is to quickly and efficiently filter the strings in order to remove regions which cannot contain a (L, r, d)-repetition applying the necessary conditions described in Section 2 and keeping only the regions which satisfy these conditions. We compute the minimum number p r of repeated kfactors each motif has to contain to possibly be part of a (L, r, d)-repetition. A set of p r k-factors contained in a region of length L is called a p r -set ≤L .

Multi-Nimbus searches for all the p r -sets ≤L that are repeated in r of the m strings. All the positions where a substring of length L contains a p r -set ≤L repeated at least once in r strings are kept by the filter, the others are rejected. The algorithm extracting the positions of all the p r -sets ≤L is presented in Figure 2.

Multi-Nimbus Initialise()

1. for g in [(p r -2)k, L -2k] 2.
for all (k, g)-bi-factors bf 3.

Multi-Nimbus Recursive(gk, positions(bf ), 2, 4.

f irstKF actor(bf )) Multi-Nimbus Recursive (gmax, positions, nbKF actors, f irstKF actor) 1. if |Strings(positions)| < r then return // not in enough strings 2. if nbKF actors = p r then save positions and return 3. for g in [(p r -(nbKF actors + 1)) × k, gmax] // possible gaps length 4.

for (k, g)-bi-factors bf starting with firstKFactor 5.

positions = intersection(positions, positions(bf )) 6.

Multi-Nimbus Recursive (gk, positions, nbKF actors + 1, 7.

f irstKF actor)

[i, i + k -1]s[i + k + g, i + 2 × k + g -1]
is a bi-factor occurring at position i in s. For the sake of simplicity, we also use the term bi-factor omitting k and g.

For instance, the (2, 1)-bi-factor occurring at position 1 in AGGAGAG is GGGA. The bi-factors occurring in at least r strings are indexed by means of a bi-factor array (presented in Section 4) that allows us to have access in constant time to the bi-factors starting with a specified k-factor. The main idea is to first find repeated bi-factors with a big gap g that may still contain

(p r -2) k-factors (g ∈ [(p r -2)k, L -2k]
). We call these border bi-factors.

A border bi-factor is a 2-set ≤L that we then try to extend to a p r -set ≤L . To extend a i-set ≤L to a (i + 1)-set ≤L , we find a repeated bi-factor (called an extending bi-factor ) starting with the same k-factor as the border bi-factor of the i-set ≤L and having a gap length shorter than all the other gaps of the bi-factors already composing the i-set ≤L . The occurring positions of the (i + 1)-set ≤L are the intersection of the extending bi-factor positions and of the positions of the i-set ≤L . An example of this construction is presented in Figure 3.

Step 1

Step 2

Step 3 During the second step, we add a bi-factor starting with the same k-factor (here called A), included inside the first one, and we intersect the positions. We repeat this once again (as p r = 4) and obtain a 4-set ≤L occurring in strings 3 and 5. Actually, it is not only the string numbers that are checked during the intersection but also the positions in the strings, not represented in this figure for clarity.

In order to extract all the possible p r -set ≤L , we iterate the idea described above on all the possible border bi-factors: all bi-factors with gap length in [(p r -2)k, L -2k] are considered as a possible border of a p r -set ≤L . Furthermore, while extending a i-set ≤L to a (i+1)-set ≤L , all the possible extending bi-factors have to be tested.

If the user runs Multi-Nimbus searching for very frequent repetitions (that is, r is large) with a high rate of substitutions (d is also large), the program has to choose a small value of k to avoid p r = 0. However, a small value for k (e.g., k ≤ 4) leads to a slow and inefficient filter. In order to avoid this problem, we proceed in two steps as follows. In the first step, we start by running Multi-Nimbus assuming r = 2 (that is, with a weaker condition), which in general allows to increase the value of k while improving the sensitivity and the execution time. In the second step, the remaining strings are filtered using the initial parameters asked by the user. This results in an efficient strategy that we refer to as the double pass strategy. More in general, the best choice of the value of k and the consequent value of p r changes from one application to another. Depending upon the parameters L, r and d as well as on the input size, there exist values of k that are not useful in practice, either because they are so small that they lead to a (too strong) condition which becomes hard to check, or conversely because they lead to a too weak condition (the reader can think of p r = 0 as an extreme case) which is fast to check but useless. In between these cases, there is in general a range of "reasonable" values for k. This range in practice is not so big, as we show in Section 5.3 together with the performances of the filter the different values lead to. As we shall see, the best k to be used is hard to establish a priori, as one should find the best tradeoff between specificity and speed, which depends on the specific application. Thus, the choice for k may be set by the program or left to the user.

Finding (L, r, d)-repetitions inside a single string: Mono-Nimbus Multi-Nimbus needs only a few modifications to obtain the Mono-Nimbus algorithm finding (L, r, d)-repetitions occurring inside the same string instead of occurring distributed over r ≤ m strings. Theorem 3 remains true for a single string. The fact that a repetition is distributed over more than one string or not, does not change the value of p r . The modifications to apply on the Multi-Nimbus algorithm in order to obtain Mono-Nimbus are the following:

• index the bi-factors occurring at least r times (instead of occurring in at least r strings), • the recursive part of the algorithm stops if the number of positions is lower than r (instead of the number of strings lower than r).

The double pass strategy is the same for Mono-Nimbus as for Multi-Nimbus.

Complexity analysis

Since the complexity analysis applies either to Mono or to Multi Nimbus, we use the general Nimbus name. Let us assume Nimbus has to filter m strings3 , each of length ℓ. The total input size is then n = ℓ × m.

Memory analysis, average case. For each possible gap length of the bifactors considered by the algorithm, a bi-factor array for one starting k-factor is stored in memory (taking in average O( n |Σ| k ) space, as showed in Section (4)). The possible bi-factor gap lengths are in [0, L -2k]. The total memory used by Nimbus is therefore in

O( n |Σ| k × L).
Time analysis, worst case. We here give an upper bound of the worst case time complexity. This bound is actually far from being tight in practice and is therefore a bit rough. Some future work should focus on improving this time complexity analysis, focusing on the average case and computing a more accurate complexity, possibly by means of an amortized analysis.

The overall worst case time complexity is in O(L × ℓ × n × Z p r -1 ) with Z = L × min(|Σ| k , ℓ). We shall see later (Fig. 6) that in practice computations require much smaller time consumption.

The Bi-Factor Array

Since we make heavy use of the inference of repeated bi-factors, we have used a new data structure, called a bi-factor array (BFA), that directly indexes the bi-factors of a set of strings. This data structure is a modified suffix array. Indeed, in its original version [START_REF] Manber | uffix arrays: A new method for on-line string searches[END_REF], a suffix array does not permit to find objects such as bi-factors. We therefore adapted the algorithmic ideas developed for this data structure to create the bi-factor array. A bi-factor array is a suffix array adapted for (k, g)-bi-factors (with k and g fixed) that stores them in lexicographic order (without considering the characters composing the gap region). This data structure allows to access the bi-factors starting with a specified k-factor in constant time. Nimbus never needs to consider at a same time bi-factors that start with different k-factors. Thus one BFA is constructed for each possible starting k-factor. Notice that the same data structure can be used to index bi-factors where the two factors have different sizes (say, (k 1 , g, k 2 )-factors with k 1 = k 2 ); we restrict ourselves here to the particular case of k 1 = k 2 because this is what we need for Nimbus. For the sake of simplicity, we present the algorithm of construction of the bi-factor array for one string. The generalisation to multiple strings is straightforward.

We start by recalling the properties of a suffix array. Given a string s of length n, let s[i . In general, another information is stored in the suffix array: the length of the longest common prefix (lcp) between two consecutive suffixes (s[π(i) . . . ] and s[π(i + 1) . . . ]) in the array. The construction of the permutation π of a text of length n is done in linear time and space [START_REF] Kärkkäinen | Simple Linear Work Suffix Array Construction[END_REF][11] [START_REF] Ko | Space Efficient Linear Time Construction of Suffix Arrays[END_REF]. A linear time and space lcp row construction is presented in [START_REF] Kasai | Linear-time Longest-Common-Prefix Computation in Suffix Arrays and its Applications[END_REF].

In order to compute the BFA for bi-factors starting with a given k-factor using a suffix array and its lcp, we perform the following steps:

(1) Give every k-factor a label. For instance, in a DNA sequence with k = 2, AA has the label 0, AC has label 1 and so on. An array is created containing, for every suffix, the label of its starting k-factor. In the remaining of this paper, we call a (label 1 , label 2 )-bi-factor a bi-factor of which the two k-factors are denoted by label 1 and label 2 . (2) For each suffix, the label of the k-factor occurring k + g positions before the current position is computed and stored in a virtual4 array we refer to as the predecessor label array.

(3) Construct the BFA for the bi-factor starting with the k-factor called label 1 as follows: the predecessor label array is traversed from top to bottom, each time the predecessor label value is equal to label 1 , a new position is added to the part of the BFA where bi-factors start with label 1 . Due to the suffix array properties, two consecutive bi-factors starting with label 1 are sorted with respect to the label of their second k-factor. The creation of the BFA is done such that for each (label 1 , label 2 )-bi-factor, a list of corresponding positions is stored.

We now explain in more detail how we perform the three steps above.

Labelling the k-factors. In order to give each distinct k-factor a different label, the lcp array is read from top to bottom. The label of the k-factor corresponding to the i th suffix in the suffix array, called label[i], is created as follows for i ∈ [0, n -1]:

label[i] =              0 if i = 0 label[i -1] + 1 if lcp[i] ≤ k label[i -1] ow
Giving each suffix a predecessor label. For each suffix, the label of the k-factor occurring k + g positions before has to be known. Let pred be the array containing the label of the predecessor for each position. It is filled as

follows: ∀ i ∈ [0, |s| -1], pred[i] = label π -1 π[i]-k -g (π -1 [p]
is the index in the suffix array where the suffix s[p . . . ] occurs). Actually, the pred array is not stored in memory. Instead, each cell is computed on line in constant time.

An example of the label and pred arrays is given in Figure 4.

Creating the BFA for a bi-factor starting with a k-factor called label 1 .

The BFA contains in each cell a (label 1 , label 2 )-bi-factor. We store the label 1 and label 2 values and a list of positions of the occurrences of the (label 1 , label 2 )-bi-factor. This array is constructed on the observation that for all i, the complete suffix array contains the information that a (pred The space complexity is in O(n), as all the steps use linear arrays. Furthermore, no more than four arrays are simultaneously needed, thus the effective memory used is 16×n bytes. The first two steps are done in O(n) time (simple traversals of the suffix array). The last step is an enumeration of the bi-factors found (no more than n). The last step is therefore in O(n) as well. Hence the total time construction of the suffix array is in O(n). With the following parameters: L = 100 and k = 6, Nimbus has to construct BFAs for around 90 different g values, which means 90 different BFAs. This operation takes, for strings of length 1Mb, around 1.5 minutes on a 1.2 GHz Pentium 3.

Testing the Filter

We tested a prototype of Nimbus on a 3 GHz Pentium 4 with 1Gb of memory. We performed several kinds of tests in order to evaluate how specificity and performances of the filter depend upon the input parameters. To this purpose, we used real biological sequences as well as randomly generated ones with planted repetitions. By "planted" repetitions, we mean that we manually inserted the repetitions in the strings so that we could know all the positions corresponding to a repetition, and hence evaluate how many positions, other than these, the filter (uselessly) keeps. The results of the tests made with randomly generated sequences are actually the average computed between several distinct runs (whose exact number changes from test to test since the variance was not always the same).

Speed and space with respect to input size

The first test is an evaluation of the actual time and space required by Nimbus which, as we show here, is in practice better than the worst case analysis made in Section 3. Figure 6 shows the time and memory usage in function of the input data length and of the kind of data. One can observe that the space complexity is linear with the input length, although with a high multiplicative constant factor. The execution time depends on the type of sequences that are given as input to the algorithm. In Figure 6, we report the results obtained with the same parameters but with two distinct input files. On the left, the result refers to an experiment where the input data set consists of ten copies of the same string: in this way, all the positions must be kept by the filter, representing the worst case time complexity, which is a polynomial of degree p r . On the right, the input strings are ten distinct strings, and the time complexity is clearly linear. We should bear in mind that these first results are from a prototype version of Nimbus that could be improved to further reduce the time and space used.

Specificity, time and memory with respect to input type and r

The second test aims at evaluating the performance of Nimbus in terms of the percentage of input strings that are filtered, and how much remains thatdoes not corresponds to repetitions, while still keeping an eye on the time and space costs. These values are computed using r = 2, 3, 4 and four different types of input sequences: three randomly generated with different numbers of planted motifs (respectively 2, 5, and 100), and one real biological sequence.

In general, both time and false positive ratios (that is, the ratio, computed on random strings with planted motifs, of non filtered data that are not part of a real motif) are low, showing a good performance of the filter.

In Figure 7, we present the behaviour of the filter for the four kinds of input strings of length 1 Mb. The first three strings are random and contain respectively 2, 5 and 100 motifs of length 100, pairwise distant by 10 substitutions. For each of these three strings, we ran Nimbus in order to filter searching for motifs of length L = 100 occurring at least r = 2, 3 and 4 times with less than d = 10 substitutions. The last string is the genomic sequence of the Neisseria meningitidis strain MC58. Neisseria genomes are known for the abundance and diversity of their repetitive DNA in terms of size and structure [START_REF] Tettelin | Complete Genome sequence of Neisseria Meningitidis Serogroup B Strain MC58[END_REF]. The size of the repeated elements range from 10 bases to more than 2000 bases, and their number, depending on the type of the repeated element, may reach more than 200 copies. This fact explains why the N. meningitidis MC58 genomic sequence has already been used [START_REF] Kolpakov | MREPS: Efficient and Flexible Detection of Tandem Repeats in DNA[END_REF] as a test case for programs identifying repetitive elements. We ran Nimbus on this sequence in order to filter repetitions of length L = 100 occurring at least r = 2, 3 and 4 times with less than d = 10 substitutions.

For each input sequence, we compute the actual space costs and report (line "Memory Used") the constant multiplicative factor on the linear complexity.

For r = 2, we used k = 6 which gives a good result: around a few minutes execution time for all the random strings. One can observe that for the MC58 sequence, the execution time is longer (15 to 23 minutes) due to its high rate of repetitions. For r = 3 and 4, to avoid a long and inefficient filter, we apply the double pass strategy described earlier, and start the filtration with r = 2 and k = 6. The time results are therefore divided into two parts: the time needed for the first pass and the one needed for the second pass. The time needed for the second pass is negligible with respect to the time used for the first one. This is due to the fact that the first pass filters from 89 % to 99 % of the string, thus the second pass works on a string at least 10 times shorter than the original one. This also explains why no extra memory space is needed for the second pass. For r = 3, the second pass uses k = 5 while for r = 4, the second pass uses k = 4. With r = 4 the necessary condition (p 4 = 5 shared 4-factors) is weaker than for r = 3 (p 3 = 5 shared 5-factors). That is why for MC58, more positions are kept while searching for motifs repeated 4 times, than for motifs repeated 3 times. Without using the double pass, for instance on MC58, with r = 3 the execution time is around 6 hours (instead of 15.9 minutes). The false positive ratio observed in practice is very low (less than 1.3 %). In general, many of the false positives occur around a (L, r, d)-repetition motif and not elsewhere in the strings.

Time and false positive rate with respect to k

We then evaluated the range of possible values to be chosen for k, and the relative performances these lead to. We generated three random sequences, each of length 50000, with planted (1000, 3, 10)-repetitions. We ran Nimbus with parameters L = 100, r = 3, d = 10, and k = 3, 4, 5, 6, 7. For k = 7, we have p 3 = 0 and thus the filter is useless. For k = 3 or less there are so many shared k-factors that the filter becomes too slow. Figure 8 shows the results for the other values.

As anticipated earlier in the paper (Section 3), the observed actual range of reasonable values of k is not big. In this test, we have three such values: from the smallest k = 4 which takes more than seven hours and has a very low false positive ratio, to the fastest ( which the filter is preprocessing. On the other hand, higher values of k lead to a very fast filter which, in a negligible time, can already sensibly shrink the sequences. Finally, and not surprisingly, the memory taken by Nimbus does not depend on k. We can observe that the false positive ratio is quite low (always below 0.4%), although it grows with L. The reason for the (slow) growth is the fact that when a portion of the input sequence satisfies the necessary condition, what is kept is indeed a fragment of length L. Time grows fast with respect to L. In this case, we may actually face an exponential time complexity with respect to p r , which on its turn grows with L.

Time and false positive rate with respect to d

We now show the results of other runs where the tested parameter is d. We generated three random sequences, each of size 50000, and planted there (L = 100, r = 3, d)-repetitions with growing values of d. Hence, we ran Nimbus with parameters L = 100, r = 3, using k = 6 with d from 0 to 13 (for higher values we have p 3 = 0). The results are shown in Figure 10. We can observe that the false positive ratio, although remaining below 0.25, grows linearly with respect to d, due to the fact that the necessary condition becomes more and more weak. Concerning the time performance, the pick at d = 5 is quite interesting, especially because it was obtained in all runs for distinct randomly generated input strings. A possible explanation is that for this parameter set, d = 5 corresponds to an intermediate value between lower and higher ones. When the value of d is low, we have that p r is so high (and hence the condition strong) that very little is kept and this requires a short time to be detected (very few shared k factors). When the value of d is high, p r decreases and so does the time complexity (which asymptotically grows with p r ).

Time and false positive rate with respect to r

We now focus on the r parameter. We generated a set of R sequences. Each sequence was of length 10000 and contained an occurrence of a repetition of length 100 distant from each other by at most 3 substitutions. We let R vary from 2 to 10. For each value, we applied Nimbus using the following parameters: L = 100, d = 3, k = 5 and r = R. We did not use the double pass strategy in order to focus only on the r parameter. The results are given in Figure 11. One can observe that the parameter r has not a big influence on time and specificity. Memory usage (data not shown) is also stable. One may however notice a light increase in the execution time for r = 3. This is due to the fact that, in this case, the necessary condition is high (p 3 = 16) and the number of recursive calls (see figure 2) is important. However, the Figure may not accurately represent how the behaviour of Nimbus varies with parameter r. Indeed, the previous test would not have been possible with, for instance, d = 10. Indeed in such a case, using r = 4 would have led to p 4 = 0.

Thus, the parameter r may have a strong influence, since increasing it can lead to a useless filter if the necessary condition p reaches zero. However, for valid ranges of r, its influence on the performance of Nimbus is small.

Using the filter

In this section, we show two distinct possible applications of Nimbus. The first concerns the inference of long biased repetitions, and the second multiple alignments.

Filtering for finding long repetitions

When inferring long approximate motifs, the number of differences allowed among the occurrences of a motif is usually proportional to the length of the motif. For instance, for L = 100 and allowing for as many as L/10 substitutions, one would have d = 10 which is high. This makes the task of identifying such motifs very hard and, to the best of our knowledge, no exact method for finding such motifs with r > 2 exists. Yet such high difference rates are common in molecular biology. The Nimbus filter can efficiently be used in such cases as it heavily reduces the search space. We now show some tests that prove this claim. For testing the ability of Nimbus concerning the inference of long approximate repetitions, we ran an algorithm for extracting structured motifs called Riso [START_REF] Carvalho | A Highly Scalable Algorithm for The Extraction of Cis-regulatory Regions[END_REF] on a set of 6 sequences of total length 21 kB for finding motifs of length 40 occurring in every sequence with at most 3 substitutions pairwise. Using Riso, this test took 212 seconds. By previously filtering the data with Nimbus, the same test took 3.9 seconds. The filtering time was 0.1 seconds. The use of Nimbus thus enabled to reduce the overall time for extracting motifs from 212 seconds to 4 seconds.

Filtering for finding multiple local alignments

Multiple local alignment of r sequences of length n can be done with dynamic programming using O(2 r n r ) time and O(n r ) space. In practice, this complexity limits the application to a small number of short sequences. A few heuristics, such as MULAN [START_REF] Ovcharenko | Mulan: Multiple-Sequence Local Alignment and Visualization for Studying Function and Evolution[END_REF], exist to solve this problem. One alternative exact solution could be to run Nimbus on the input data so as to exclude the non relevant information (i.e. parts that are too distant from one another) and then to run a multiple local alignment program. The execution time is hugely reduced. For instance, on a file containing 5 random strings of cumulated size 1 Mb each containing an approximate repetition5 , we ran Nimbus in approximately 1.5 minutes. On the remaining sequences, we ran a tool for finding functional elements using exact multiple local alignment called Glam [START_REF] Frith | Finding Functionnal Sequence Elements by Multiple Local Alignement[END_REF]. This operation took about 25 seconds. Running Glam without the filtering, we obtained the same results6 in more than 7.5 hours. Thus by using Nimbus, we reduced the execution time of Glam from 7.5 hours to less than 2 minutes.

Conclusions and future work

We presented a novel lossless filtration technique for finding long multiple approximate repetitions common to several strings or inside one single string. The filter localises the parts of the input data that may indeed present repetitions by applying a necessary condition based on the number of repeated k-factors the sought repetitions have to contain. This localisation is done us-ing a new type of seeds called bi-factors. The data structure that indexes them, called a bi-factor array, has also been presented in this paper. It is constructed in linear time. This data structure may be useful for various other text algorithms that search for approximate instead of exact matches. The practical results obtained show a huge improvement in the execution time of some existing multiple sequence local alignment and pattern inference tools, by a factor of up to 100. The theoretical complexity presented in the paper is a rough upper bound of the worst case. The results show that in practice the time consumption is much smaller than the theoretical complexity. Future work thus includes obtaining a better analysis of this complexity.

Other important tasks remain, such as filtering for repetitions that present an even higher rate of substitutions, or that present insertions and deletions besides substitutions. One idea for addressing the first problem would be to use bi-factors (and the corresponding index) containing one or two mismatches inside the k-factors. In the second case, working with edit instead of Hamming distance implies only a small modification on the necessary condition and on the algorithm but could sensibly increase the execution time observed in practice.

Fig. 1 .

 1 Fig. 1. Example of stained positions and shared k-factors for r = 2. Two strings of length L = 16 distant by d = 2 substitutions: the number of stained positions (represented by black squares) is 2. Let us consider, say, k = 3 for which we have that L kd = 16 3 -2 = 3. The figure shows indeed one of the worst case scenario in which the placement of the stains limits the number of shared non-overlapping k-factors (represented by hashed rectangles) to 3.

Fig. 2 .

 2 Fig. 2. Extract the positions of all the p r -sets ≤L

Fig. 3 .

 3 Fig.3. Example of the construction of a 4-set ≤L . In the first step, we find a bi-factor occurring at least once in at least r = 2 strings among m = 8 strings. During the second step, we add a bi-factor starting with the same k-factor (here called A), included inside the first one, and we intersect the positions. We repeat this once again (as p r = 4) and obtain a 4-set ≤L occurring in strings 3 and 5. Actually, it is not only the string numbers that are checked during the intersection but also the positions in the strings, not represented in this figure for clarity.

  . . ] denote the suffix starting at position i. Thus s[i . . . ] = s[i, n -1]. The suffix array of s is the permutation π of {0, 1, . . . , n-1} corresponding to the lexicographic order of the suffixes s[i . . . ]. If ≤ l denotes the lexicographic order between two strings, then s[π(0) . . . ] ≤ l s[π(1) . . . ] ≤ l . . . ≤ l s[π(n -1) . . . ].

Fig. 4 .Fig. 5 .

 45 Fig. 4. Suffix array completed with the label and the pred arrays for k = 2 and g = 1 for the text ACAACCAC.of positions of the previous bi-factor if label 2 of the latter is equal to label[i]. An example of a BFA is given in Figure5.

Fig. 6 .

 6 Fig. 6. Time and space spent by Nimbus w.r.t. the input data length. The parameters are L = 100, k = 6, d = 7, r = 3 which implies p 3 = 6. The input file contains 10 strings of equal length. On the left, the same string is replicated 10 times and given as input, whereas on the right, we used 10 distinct random strings.

Fig. 7 .

 7 Fig.7. Nimbus behaviour on four types of strings while filtering in order to find r = 2, 3 and 4 repetitions. The first three strings are random and have length 1 Mb. They contain respectively 2, 5 and 100 motifs of length 100 distant pairwise by 10 substitutions. The last string is a segment of the genomic sequence of the strain MC58 of Neisseria meningitidis also of length 1 Mb. "FP" stands for False Positive ratio and the "Kept" lines give the number and percentage of positions kept.

Fig. 8 .

 8 Fig. 8. Speed and false positive ratio for several values of k.

5. 4 Fig. 9 .

 49 Fig. 9. Time and false positive ratio with respect to L.

Fig. 10 .

 10 Fig. 10. Time and false positive ratio with respect to d.

Fig. 11 .

 11 Fig. 11. Time and false positive ratio with respect to r.

In the case of Mono-Nimbus, m = 1.

In practice these values are not stored but computed on the fly.

Repetitions of length 100 containing 10 substitutions pairwise.

Since Glam handles edit distance and Nimbus does not, in the tests we have used randomly generated data where we planted repetitions allowing for substitutions only, in order to ensure that the output would be the same and hence the time cost comparison meaningful.

Acknowledgement

The authors would like to thank the anonymous referee for his careful reading and the very helpful comments that sensibly improved the paper.

Supported by the ACI Nouvelles Interfaces des Mathématiques π-vert project of the French Ministry of Research. 2 Supported by the ACI Nouvelles Interfaces des Mathématiques π-vert project of the French Ministry of Research, the ARC BIN project from the INRIA, and the ANR project REGLIS.


Aline Neves

Maria D Miranda

Mamadou Mboup

ALGEBRAIC PARAMETER ESTIMATION OF DAMPED EXPONENTIALS

 which suffers in the presence of noise. In this paper we propose an algebraic method for the parameter estimation. The method, differently from the modified Prony method, is considerably robust to noise. The comparison of both through simulations confirm the good performance of the algebraic method.

INTRODUCTION

The parameter estimation of exponentially damped sinusoids is a problem that arises in many data applications. For example, models that describe heat and chemical components diffusion, or time series related to biomedical signals are represented by a sum of exponentials. The problem is specially interesting since frequency estimation is fundamental in signal processing. In addition, this kind of signal is transient due to the damping factor, what makes the problem challenging.

Usually, these signals are represented by:

x(t) = p ∑ j=1
a j e α j t [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF] where a j are amplitudes and α j are constants that may be complex and usually have a negative real part. Estimation of a j and α j is well known to be numerically difficult [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF]. Since, in practice, x(t) can never be exactly observed due to noise and measurement errors, a first solution would be to use a minimum square criteria where the error is e(t) = x(t)y(t) with y(t) = x(t) + η(t) being the observed signal and η(t) a random Gaussian noise. The estimation of a j and α j through this method, however, results in a nonlin- ear problem.

The study of this kind of parameter estimation has a long and rich history. The classical technique of Prony, for example, extracts sinusoid or exponential signals from time series data by solving a set of linear equations for the coefficients of the recurrence equation satisfied by these signals [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF]. It is a technique closely related to Pisarenko's method which finds the smallest eigenvalue of an estimated covariance matrix [START_REF] Ouibrahim | Prony, Pisarenko, and the matrix pen-cil: a unified presentation[END_REF]. Prony's method, however, is well known to perform poorly in the presence of noise. On the other hand, Pisarenko's method is consistent but inefficient for estimating sinusoid signals and inconsistent for estimating damped sinusoids or exponential signals. In addition, choosing initial values, ill-conditioning when two or more α j are close and the order of the estimated values are just some of the encountered difficulties in solving the problem.

Searching for a better performance, Osborne [START_REF] Osborne | Some special nonlinear least squares problems[END_REF] proposed a modified Prony algorithm equivalent to maximum likelihood estimation. The method was generalized to estimate any function which satisfies a difference equation with constant coefficients. In [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF], the authors apply the method to fitting sums of exponential functions. The method is shown to be relatively insensitive to starting values and also solves the ill-conditioning problem as far as the convergence of the algorithm is concerned, but may return a pair of damped sinusoids in place of two exponentials which are coalescing. Moreover, it suffers in the presence of noise even though it performs better then other existing methods.

In this paper, we propose an algebraic estimation method for the exponential fitting of the observed signal that results from an identification/estimation theory based on differential algebra and operational calculus [START_REF] Fliess | An algebraic framework for linear identification[END_REF][START_REF] Fliess | Questioning some paradigms of signal processing via concrete examples[END_REF]. Application of this approach in other contexts as the deconvolution of BPSK (Binary Phase Shift Keying) and QPSK (Quadrature Phase Shift Keying) signals [START_REF] Neves | An algebraic identification method for the demodulation of QPSK signal through a convolutive channel[END_REF] and the demodulation of CPM (Continuous Phase Modulation) signals [START_REF] Neves | An algebraic receiver for full response CPM demodulation[END_REF], among others [START_REF] Mboup | Parameter estimation via differential algebra and operational calculus[END_REF], all associated with the estimation of different parameters, showed that the resulting methods are specially robust to noise, enabling a good estimation of the desired parameters even at very low signal to noise ratios [START_REF] Neves | Identification Algébrique et Déterministe de Signaux et Systèmes à Temps Continu: Application à des Problèmes de Communication Numérique[END_REF]. Since noise always largely degrades the performance of existing methods for the parameter estimation of exponentially damped sinusoids, we expect the algebraic approach to enhance the results obtained under these conditions.

The paper is organized as follows. Section 2 briefly reviews the modified Prony method. Section 3 proposes an algebraic method for the estimation of the desired parameters. Considering the same framework used in [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF], we study the case of a signal given by the sum of two exponentials and briefly show how to generalize the resulting method. Simulations and discussions are presented in section 4, where both methods are compared. The paper concludes with section 5.

THE MODIFIED PRONY METHOD

The modified Prony algorithm is equivalent to maximum likelihood estimation for Gaussian noise. It was generalized in [START_REF] Osborne | A modified Prony algorithm for fitting functions defined by difference equations[END_REF] to estimate any function that satisfies a difference equation with coefficients that are linear and homogeneous in the parameters. In [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF], the algorithm was applied to exponential fitting. The authors start considering a signal as shown in [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF], which satisfies a constant coefficient differential equation. In the discrete time case, it is known that this signal can be written as an AR (Auto-Regressive) signal which satisfies the difference equation:

p+1 ∑ k=1 ϑ k y(n -1 -k) = 0 (2)
for some suitable choice of ϑ k with y(n) being equally spaced samples of y(t). From ( 2), setting y = [y(n) y(n -1) ... y(n -N)] T , we can write

X T (ϑ )y = 0 (3)
where X(ϑ ) is the convolution matrix

X(ϑ ) =          ϑ 1 . . . . . . ϑ 1 ϑ p+1 . . . . . . ϑ p+1         
The algorithm results from the minimization of the square error φ (a, α) = (yx) T (yx) subject to the constraint ϑ T ϑ = 1. A necessary condition for it to be attained is:

(B(ϑ ) -λ I) ϑ = 0 ( 4
)
where λ is the Lagrange multiplier and

B(ϑ ) i j = y T X i X T X -1 X T j y + -y T X X T X -1 X T i X j X T X -1 X T y with X j = ∂ X(ϑ )/∂ ϑ j .
The problem is then a nonlinear eigenvalue problem. The modified Prony algorithm solves (4) using a succession of linear problems converging to λ = 0. Given an estimate ϑ k the solution θ , solve:

B(ϑ k ) -λ k+1 I ϑ k+1 = 0 ϑ k+1 T ϑ k+1 = 1 (5)
with λ k+1 the nearest to zero of such solutions. Convergence is accepted when λ k+1 is small compared to B(ϑ ) .

In order to recover α j having estimated ϑ , we need to obtain the roots of the characteristic polynomial p ϑ = ∑ p+1 k=1 ϑ k z k-1 which will result in ζ j = n(1 -e -α j /n ). The final step is then

α j = -nlog(1 -ζ j /n) = n ∞ ∑ j=1 j -1 (ζ j /n) j
The modified Prony method is relatively insensitive to starting values. It may return a pair of damped sinusoids in place of two exponentials which are coalescing. More details about this method can be found in [START_REF] Osborne | A modified Prony algorithm for fitting functions defined by difference equations[END_REF][START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF].

ALGEBRAIC ESTIMATION METHOD

Let us consider a signal given by the sum of two complex exponentials in noise:

y(t) = a 1 e α 1 t + a 2 e α 2 t + γ + η(t) (6)
where α 1 and α 2 may be complex, γ is a constant bias perturbation and η(t) is an additive noise. We suppose that the amplitudes a 1 and a 2 are irrelevant for our purpose and consider them equal to one. Our objective is, thus, to estimate α 1 and α 2 .

Considering that we will observe the signal y(t) during a finite time interval T ≤ t ≤ T + λ , the noise η(t) can be decomposed as:

η(t) = γ λ T + η λ T (t) (7)
where the constant γ λ T represents its mean (average) value and η λ T (t) is a zero-mean term. Returning to (6), we can thus consider, without any loss of generality, that in the perturbation γ + η(t), η(t) is zero-mean.

The algebraic method proposed is based on differential algebra and operational calculus. Therefore, the first step is to obtain a differential equation satisfied by [START_REF] Neves | An algebraic identification method for the demodulation of QPSK signal through a convolutive channel[END_REF]. Ignoring the zero-mean noise term η(t) for the moment, the signal

y(t) = a 1 e α 1 t + a 2 e α 2 t + γ (8)
satisfies the following differential equation

y (3) = θ 1 y (2) + θ 2 y (1) (9)
where the super-script (i) denotes de order of differentiation of y(t) with respect to time, θ 1 = α 1 + α 2 and θ 2 = -(α 1 α 2 ). Observing [START_REF] Neves | Identification Algébrique et Déterministe de Signaux et Systèmes à Temps Continu: Application à des Problèmes de Communication Numérique[END_REF], it is important to note that the differential equation has, as parameters, functions of the desired α 1 and α 2 . In addition, (9) does not depend on the structured per- turbation γ. Since it is a constant, it can easily be eliminated through differentiation.

The next step is to obtain the Laplace transform of (9):

s 3 ŷ -s 2 y 0 + s ẏ0 + ÿ0 = θ 1 (s 2 ŷ -(sy 0 + ẏ0 )) + θ 2 (s ŷ -y 0 ) (10
) where y 0 , ẏ0 and ÿ0 are the initial conditions of y(t) in the time interval being considered. These unknown constants can also be viewed as structured perturbations that can be eliminated by taking the derivative of [START_REF] Osborne | A modified Prony algorithm for fitting functions defined by difference equations[END_REF] recursively with respect to the variable s three times:

(s 3 ŷ) (3) = θ 1 (s 2 ŷ) (3) + θ 2 (s ŷ) (3) (11)
Now, we can use (11) to generate a system of equations that will allow us to estimate the desired parameters. Since we have only two unknown parameters, it suffices to differentiate (11) once more to obtain:

(s 3 ŷ) (3) (s 3 ŷ) (4) = (s 2 ŷ) (3) (s ŷ) (3) (s 2 ŷ) (4) (s ŷ) (4) θ 1 θ 2 (12)

Solving the system of equations (12) will give us an estimation of α 1 and α 2 . It is clear, however, that the signal y(t) is only known in the time-domain. Applying the inverse Laplace transform to (12), we observe that the multiplication by s in the operational domain will result in a differentiation with respect to time, in the time domain. Since differentiations are not robust operations to be calculated numerically, we can avoid them by dividing both sides of equation ( 12) by s ν . Developing the (.) k for k = 3, 4 and dividing by s ν we obtain:

   y (3) s ν-3 + 9 y (2) s ν-2 + 18 y (1) s ν-1 + 6 y s ν y (4) s ν-3 + 12 y (3) s ν-2 + 36 y (2) s ν-1 + 24 y (1) s ν    =    y (3) s ν-2 + 6 y (2) s ν-1 + 6 y (1) s ν y (3) s ν-1 + 3 y (2) s ν y (4) s ν-2 + 8 y (3) s ν-1 + 12 y (2) s ν y (4) s ν-1 + 4 y (3) s ν    θ 1 θ 2 ( 13
)
where ν is larger than the largest power of s appearing in the system. In the considered case, ν has to be larger than 3.

Returning to the time domain, we will have integral operations instead of differentiations. In general, (13) is composed by terms of the form y (i) s νj , which inverse Laplace transform is given by the νj order iterated integral

L -1 y (i) s ν-j = λ 0 dλ ν-j-1 λ ν-j-1 0 • • • λ 1 0 τ i y(τ)dτ = (-1) i (ν -j -1)! λ 0 (λ -τ) ν-j-1 τ i y(τ)dτ (14)
where λ is the integration interval.

Finally, the parameters can be estimated by solving

Pθ = Q (15) where θ = [θ 1 θ 2 ] T , Q = λ 0 1 0 0 τ q 1 q 2 κ ν-4 y(τ)dτ and P = λ 0 1 0 0 τ p 1,1 p 1,2 p 2,1 p 2,2 κ ν-3 0 0 τκ ν-2 y(τ)τdτ with κ = λ -τ and p 1,1 = -τ 2 b 1,2 + 6τκb 1,1 -6κ 2 p 2,1 = τ 2 b 1,2 -8τκb 1,1 + 12κ 2 p 1,2 = -τb 1,1 + 3κ p 2,2 = τb 1,1 -4κ q 1 = -τ 3 b 1,3 + 9τ 2 κb 1,2 -18τκ 2 b 1,1 + 6κ 3 q 2 = τ 3 b 1,3 -12τ 2 κb 1,2 + 36τκ 2 b 1,1 -24κ 3 b 1, j = ∏ j i=1 (ν -i)
Note that the estimation time λ , which depends on the observation signal y(t), has to be sufficient to the integrals to converge. In general this interval can be small, which translates into fast estimation.

In addition, when the zero-mean noise term η(t) is added to the signal, the iterated integration procedure is going to largely reduce its effect at the output. For this reason the method is significantly robust to noise.

Having estimated θ, we still have to recover α 1 and α 2 . From ( 9), α 1 and α 2 are the roots of the characteristic equa- tion:

p x = z 2 -θ 1 z -θ 2 = (z -α 1 )(z -α 2 ) (16)
The generalization of the procedure shown is quite direct. Consider a signal given by the sum of p exponentials:

y(t) = p ∑ i=1 a i e α i t + γ (17)
It will satisfy the following differential equation:

d p+1 y dt p+1 = θ 1 d p y dt p + θ 2 d p-1 y dt p-1 + ... + θ p dy dt ( 18
)
where the coefficients θ k are algebraic functions of the de- sired unknown parameters:

θ k = (-1) k-1 p ∑ i 1 , ..., i k = 1 i 1 < ... < i k α i 1 α i 2 ...α i k
Having found the respective differential equation, we can resume the necessary steps for the development of the method as follows: 1. Obtain the Laplace transform of (18), i.e., a differential equation satisfied by the observed signal, with constant coefficients that depend on the desired parameters. 2. Differentiate the resulting equation with respect to the variable s recursively p + 1 times in order to eliminate the dependence on initial conditions of y(t).

Continue to differentiate the resulting expression recur-

sively in order to generate a system of equations with a number of equations equal to the number of unknown coefficients. We will then have a system with p equations. 4. Divide all equations by a term s ν , with ν > p + 1, to obtain proper operators. 5. Return to the time-domain, calculating the necessary integrals as defined in (14). 6. Solve the resulting linear system. Once we have found the estimates of θ i , i = 1, 2, the coefficients α i can be found as the roots of the characteristic polynomial of (18).

SIMULATION RESULTS AND DISCUSSION

Comparing the proposed algebraic method with the modified Prony's method, it is interesting to observe that both start considering the differential equation that is satisfied by y(t).

The modified Prony's method, however, treats the problem in discrete time, using the corresponding difference equation, while the algebraic method treats it directly in continuous time. In addition, both methods are insensitive to initial conditions.

Let us start with the case where y(t) is given by only one exponential:

y(t) = e αt
where α is complex. The development of the algebraic method for this case follows exactly the same steps presented in section 3. Here, y(t) satisfies the following differential equation:

ÿ(t) = α ẏ(t)
where the structured perturbation γ was already eliminated through differentiation. Figures 1 and2 show the convergence of both methods with respect to time, for the estimation of the real and imaginary components of α respectively. The value of α was set to -0.2 + 0.9i, ν was considered equal to 4 and white gaussian noise was added to y(t). The presented curves are the resulting mean of 100 Monte Carlo simulations. The signal y(t) was observed during 10 seconds and we considered 1000 samples for the algebraic method and 15 for the modified Prony method. This difference comes from the fact that the algebraic method treats the signals in continuous time and needs a sufficient number of samples for the integrals to converge, while the modified Prony method uses equally spaced samples of the signal. For the later, increasing the number of points degrades the performance.

As we can see, the performances of the algebraic method are similar for SNR (Signal to Noise Ratio) of 15 dB and 0 dB. However, this is not the case for the modified Prony method. In the estimation of α at 0 dB presented here, the obtained value is out of the graphic scale. In the sequel, we will consider y(t) given by ( 6), i.e., the sum of two exponentials. Figure 3 shows the result obtained using real valued α: α 1 = -0.25 and α 2 = -0.8. The fig- ure shows the mean square error (MSE) obtained comparing the estimated values of α with the correct ones for 1000 Monte Carlo simulations. The signal was observed during an interval of 20 seconds. The Algebraic Method (AM) used 5000 samples while the Modified Prony method (MP) used 25. The parameter ν was set equal to 5.

It is clear that the modified Prony method suffers considerably with the addition of noise, differently from the algebraic method. As mentioned in [START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF], the former may give a complex value as a result of the estimation of real parameters. The error shown in figure 3 was calculated considering only the real part of the estimated parameter. The number of inaccurate estimations increases as the SNR decreases. The algebraic method does not present such a behavior, i.e., if the desired parameters are real it always returns a real valued estimation. For comparison, the figure also shows the Cramer Rao Bound (CRB). Figures 4 and5 show the error obtained in the estimation of complex valued α: α 1 = -0.25 -0.9i and α 2 = -1.2 + i. Once again we can see the good performance of the algebraic method when compared to the modified Prony method. The simulations in this case used the same scenario above: ν = 5, 5000 samples for MA and 25 for MP, 1000 Monte Carlo simulations.

It should be noted that, in the case without noise, both methods coincide estimating the desired parameters with an error that tends to zero.

CONCLUSION

We have proposed an algebraic method for the parameter estimation of a sum of exponentials. This method was compared to the modified Prony method, one of the best existing methods in the literature. The starting point of both methods is the same: a linear homogeneous differential/difference equation with constant coefficients that are algebraic functions of the desired parameters. In both cases, the estima- tions are obtained by finding the roots of a polynomial. One of the main differences is that the modified Prony is based on a nonlinear optimization while the algebraic method solves a linear system of equations.

Simulations have shown that the algebraic method is more robust to noise than the modified Prony method. The problem of yielding complex results when estimating real parameters was not observed in the simulations with the algebraic method, while for the modified Prony method it occurred in simulations where SNR ≤ 20. The number of inaccurate estimations increased with the decrease of the SNR value.

Figure 1 :

 1 Figure 1: Estimation of Re{α}

Figure 2 :

 2 Figure 2: Estimation of Im{α}

1 CRB, α 2 AM, α 1 AM, α 2 MP, α 2 MP, α 1

 121221

Figure 3 :

 3 Figure 3: Comparison of the performances of the algebraic method and the modified Prony method in identifying α 1 = -0.25 and α 2 = -0.8

Figure 4 : 2 }Figure 5 :

 425 Figure 4: Comparison of the performances of the algebraic method and the modified Prony method in identifying Re{α}


Randa Khemiri
email: .randa.khemiri@gmail.com

Nejmeddine Bahri

Fatma Belghith

Fatma Ezahra Sayadi

Mohamed Atri

Nouri Masmoudi

Fast Motion Estimation for HEVC Video Coding

Keywords: HEVC, HM16.2, Fast configuration (ECU, ESD, CFM), Motion Estimation, Diamond search pattern I

In this paper, a fast configuration for Motion Estimation (ME) is described in order to reduce the computational time of the new High Efficient Video Coding (HEVC). This configuration uses the Coded Block Flag (CBF) Fast Method (CFM), the Early Coding Unit (CU) termination (ECU) and the Early Skip Detection (ESD) modes. The Diamond Pattern is used as a search algorithm for ME in the encoding process. Compared to the latest original reference software test model (HM) 16.2 of the HEVC, experimental results had showed that the complexity is reduced, in average, by 56.75% with a small bit-rate and PSNR degradation.

INTRODUCTION

Nowadays, the proposed HEVC [START_REF]ITU-T and ISO/IEC JTC 1, High Efficient Video Coding (HEVC)[END_REF] is the most recent standard used on video compression [START_REF] Sullivan | Overview of the high efficiency video coding (HEVC) standard[END_REF]. When compared to the previous standards, it can achieve as equivalent visual quality as H.264/MPEG-4 AVC encoders, by only using 50% less bitrate [START_REF] Ohm | Comparison of the coding efficiency of video coding standardsincluding high efficiency video coding (HEVC)[END_REF]. In HEVC, the motion estimation is the most important part of video compression occupying around 70% of the total encoding time [START_REF] Jongho | An SAD-Based Selective Bi-prediction Method for Fast Motion Estimation in High Efficiency Video Coding[END_REF]. It is the major factor in the compression efficiency as it occupies a large part of the whole coding arithmetic. In order to reduce the encoding time, several researchers studied ME-algorithm to optimize the test zonal search (TZ search). For reducing the encoding time, several researchers studied the motion estimation algorithm to optimize the test zonal search (TZ search). Purnachand et al. [START_REF] Purnachand | Improvements to TZSearch motion estimation algorithm For multiview video coding[END_REF] replaced the 8-point diamond and 8-point square with a 6point hexagonal in the TZS algorithm. That had reduced the encoding time by almost 50% compared to TZS algorithm with the same PSNR and bit-rate. In addition, the determination of the most optimal CU partition for each Coding Tree Unit (CTU) and the best Prediction Unit (PU) mode for each CU causes an increase in computational complexity. In the literature, many efforts had been focusing on exploring algorithms established on fast decision. Liquan et al. [START_REF] Liquan | Correspondence:an effective CU size decision method for HEVC encoders[END_REF] had presented an algorithm based on skipping some specific depths which are rarely used in the previous frame. The proposed algorithm reduces 26% in the execution time with a slight increase in bit-rate as well as a negligible coding efficiency loss.

This paper adopts some options for fast mode decision used for CU partitioning. The rest of this paper is organized as follows. In Section 2, an overview of the TZ Search algorithm, which presents the motion estimation algorithm, is provided. Section 3 presents the proposed configuration for the Random-Access-main profile (RA) to reduce the execution time. In Section 4 a comparison between this algorithm and the original one (HEVC reference software HM16.2) [START_REF]Software reference test model HM16[END_REF] is discussed. Finally, a conclusion and some future works are presented in Section 5.

II. TZ SEARCH ALGORITHM

The TZSearch is the algorithm used in the ME process in HEVC. This algorithm is divided into four distinct main steps: The PMV (Prediction Motion Vector), the first research using the diamond or the square pattern, the raster search and the refinement. These steps are summarized as follows [START_REF] Belghith | A new fast motion estimation algorithm using fast mode decision for high-efficiency video coding standard[END_REF][START_REF] Kibeya | TZSearch Pattern Search Improvement for HEVC Motion Estimation Modules[END_REF]:

A. Prediction Motion Vector(PMV)

The TZS algorithm employs four predictors: the left predictor, the up, and the upper right ones to compute the median predictor of the corresponding block. Max( A,Max( B,C ) ). This is the starting point for the next step.

B. Initial Grid Search

After selecting the start position, the second step is to perform the first search and to determine the search range and the search pattern. This step consists of finding the search window using a diamond or square pattern as shown in Figure 2(a) and Figure 2(b) with distances ranging from one to "searchrange" (For HM16.2 , the maximum window is set to 64), the distance values are multiples of two. Thus, these twosearch patterns use 8 points at each loop. The distance of the minimum distortion point is stored in the "Bestdistance" variable. The patterns used are either square or diamond depending on the configuration file.

(a).Diamond search pattern (b). Square search pattern (c). Raster search pattern with iRaster equals to 4

C. Raster Search

When the previous search is completed, the next step is the raster search, which chooses the best distance corresponding to the best-matched point from the last search. Depending on this distance, three different cases are presented. First, if the best distance is equal to zero, then the process is stopped. Second if 1<Bestdistance<iRaster, a refinement stage is proceeded directly. In fact, the "iRaster" is a parameter that represents the central point and the first resulting from the research that it should not exceed. Third, if Bestdistance>iRaster which is set appropriately, a raster scan is achieved using the iRaster's value as the stride length. The raster search is performed when the difference obtained between the motion vectors from the first phase and the starting position is too large. This raster search is carried out on the complete search window. A FullSearch algorithm is shown in Figure 2 (c) for a 16x16 search window with iRaster equal to 4.

D. Star/Raster Refinement

When the motion vectors distance obtained in the previous steps is not equal to zero, the refinement step is performed. There are two types of refinements: * Star refinement: The starting point of this search is the best point resulting from the previous searches. This step can also be done using either diamond or square pattern with different distances from one to "searchrange". This is done by multiplying the distance by two in each iteration. Hence, at a distance equal to one, the search for two adjacent points is applied in order to choose the best estimated motion vector having the best Sum of Absolute Difference (SAD). * Raster refinement: Concerning this type of refinement, the starting point of this search is the best point resulting from the previous searches. This can be done using either diamond or square pattern by dividing the distance by two in each iteration up to having a distance equal to one. Where, the search for two adjacent points is applied and the process is stopped.

III. PROPOSED FAST CONFIGURATION

Many fast mode decisions were used in this work to reduce computational complexity of the ME algorithm using diamond search pattern, which reduces the computational complexity of all the encoder. Early CU termination (ECU), Early Skip Detection (ESD) and Coded Block Flag algorithm (CBF) are the methods used. An overview of the different configuration options adopted in HEVC test model HM16.2 is presented: ECU: Enables the use of early CU determination, skipped CUs will not be split further [START_REF] Choi | Coding Tree Pruning Based CU Early Termination[END_REF].

ESD:

Enables the use of CBF-based fast encoder mode. Once a 2Nx2N CU has been evaluated, if the RootCbf is 0, further PU splits will not be evaluated. Check inter 2Nx2N first and terminate if the motion vector difference is equal to zero and coded block flag equal to zero [START_REF] Yang | Early SKIP Detection for HEVC[END_REF].

CFM:

Enables the use of early skip detection and the skip mode will be tested before any other mode. Terminate encoding decision if the PU partitioning has a coded block flag equal to zero [START_REF] Bossen | HM Software[END_REF].

All of these additional modes influence only the inter picture prediction mode decision.

IV. EXPERIMENTAL RESULTS

A. Experimental conditions

In order to evaluate the performance of the fast configuration, the recent available HEVC reference software (HM 16.2) was implemented with the configuration Random Access Main profile (RA), with all previously described ECU, ESD, and CBF methods. The proposed implementation was compared to the original algorithm in terms of PSNR, bit-rate and encoding time speedup. The maximum search range used in this work is 64 with maximum CU Size equal to 64 and maximum CU partition depth set to 4. The simulations were carried on windows 8 OS platform with Intel ®core TM i7-3770 @ 3,4 GHz CPU and 12 GB RAM. For the different classes A, B, C and D with four resolutions [START_REF] Dıaz-Honrubia | HEVC a review, trends and challenges[END_REF] (2560x1600 /1920x1080/ 832x480 and 416x240 formats), the number of frames is set to 50. The fast algorithm was evaluated for quantization parameter (QP) 22 and 37.

B. Evaluation criteria

The result of the fast configuration used in this work is shown in Table I. The data analysis of encoding time comparison is calculated using the following formula:

 × 100 (%) (1)
As consequence the fast HEVC configuration allows saving up to 56.75% of encoding time obtained exactly for the BQTerrace 1920x1080 sequence. Further, for all test sequence used, the speedup is more important for higher QPs due to the fact that skip mode is more likely to be chosen for important values of QP [14]. The improvement is more important for videos containing low motion activities contrarily for high activity sequences such as PeopleOnStreet and RaceHorses.

Table II shows the PSNR (dB), PSNR (dB), bit-rate (kbps) and the bit-rate(%) values obtained for original HEVC configuration compared to the fast one, used in ME for QP 22 and 37. PSNR and bit-rate are calculated using the both following formulas:

 (dB) (2)
 ×100(%) (3)
The results show a slice degradation in video quality in terms of PSNR, around 0.083 dB, with a little decrease in the bit-rate, around 1.25%.

Which proves that the adopted configuration is much efficient that the original, since it offered an important time saving, in addition to the decrease in bit-rate.

V. CONCLUSION

The new HEVC standard provides a significant improvement in video quality, especially for high-resolution videos. ME is the most challenging and consuming stage that causes a strong increase in computational complexity of all the encoder. Due to its new quadtree-based coding unit (CU) block partitioning structure, variable block size and multiple reference frames. In this paper, a fast configuration based on optimizing the motion estimation process using a diamond search and a fast mode decision algorithm for CU partitioning have been presented to reduce the computational complexity of the encoder. The simulation results prove that the fast HEVC encoder configuration allows saving up to 56.75% of encoding time obtained with a negligible degradation in video quality and decreasing the bit-rate.

The objective is to additionally reduce the motion estimation time with this fast configuration and with the parallelization of the search mode used here.

Figure 1 .

 1 Figure 1. MV adjacent a current PU The median is computed using the Eq. 1: Median( A, B, C ) = A + B + C -Min( A, Min( B, C ) ) -Max( A, Max( B, C) ). This is the starting point for the next step.

Figure 2 .

 2 Figure 2. Diamond/Square/Raster search pattern

TABLE I

 I

	.	ENCODING TIME COMPARISON RESULTS BETWEEN FAST AND
			ORIGINAL HEVC CONFIGURATION
					Encoding
	Class		Sequences	QP	Time (sec)
					T(%)
	Class A	PeopleOnStreet	22	25
	2560x1600		37	49
			Average		37
	Class B		BQTerrace	22	29
	1920x1080		37	84.53
			Average		56.75
	Class C		RaceHorses	22	11.46
	832x480		37	44.17
			Average		27.81
	Class D	BQSquare	22	34.64
	416x240		37	75.26
			Average		52.25

TABLE II .

 II CODING EFFICIENCY WITH THE FAST CONFIGURATION FOR 50 FRAMES

				PSNR(dB)			bit-rate(kbps)
	Sequences	QP
			Original_Conf	Fast_Conf	PSNR	Original_Conf	Fast_Conf	bit-rate(%)
	PeopleOnStreet	22	40.217	40.167	-0.04	33017.362	32859.76	-0.47
	2560x1600	37	31.21	31.173	-0.03	4738.104	4634.712	-2.182
			Average		-0.035			-1.326
	BQTerrace	22	37.99	37.938	-0.06	41970.95	41471.98	-1.18
	1920x1080	37	31.482	31.414	-0.067	1823.357	1785.09	-2
			Average		-0.063			-1.59
	RaceHorses	22	38.47	38.44	-0.029	5646.605	5621.98	-1.43
	832x480	37	29.387	29.155	-0.23	570.238	558.058	-2.138
			Average		-0.129			-1.28
	BQSquare	22	38.81	38.735	-0.075	1783.018	1767.44	-0.87
	416x240	37	29.25	29.112	-0.137	236.846	234.912	-0.81
			Average		-0.106			-0.84