Romain Tavenard
email: romain.tavenard@idiap.ch

Laurent Amsaleg

Improving the Efficiency of Traditional DTW Accelerators

Keywords: Dynamic time warping, Indexing, Lower bounds, Upper bounds, Indexing Trees

Dynamic time warping (DTW) is the most popular approach for evaluating the similarity of time series, but its computation is costly. Therefore, simple functions lower bounding DTW distances have been designed, accelerating searches by quickly pruning sequences that could not possibly be best matches. The tighter the bounds, the more they prune and the better the performance. Designing new functions that are even tighter is difficult because their computation is likely to become complex, canceling the benefits of their pruning. It is possible, however, to design simple functions with a higher pruning power by relaxing the no false dismissal assumption, resulting in approximate lower bound functions. This paper describes how very popular approaches accelerating DTW such as LB Keogh and LB PAA can be made more efficient via approximations. The accuracy of approximations can be tuned, ranging from no false dismissal to potential losses when aggressively set for great response time savings. At very large scale, indexing time series is mandatory. This paper also describes how approximate lower bound functions can be used with iSAX. Furthermore, it shows that a k-means-based quantization step for iSAX gives significant performance gains.

Introduction

Searching in very large databases of time series has received a lot of attention from researchers. Dynamic time warping (DTW) has proved to be an excellent similarity measure for time series. It has been successfully applied to domains as diverse as pattern recognition [START_REF] Munich | Continuous dynamic time warping for translation-invariant curve alignment with applications to signature verification[END_REF], life sciences and bioinformatics [START_REF] Aach | Aligning gene expression time series with time warping algorithms[END_REF][START_REF] Gavrila | Towards 3-d model-based tracking and recognition of human movement: a multi-view approach[END_REF], speech recognition [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF], etc. DTW is of great help when using any data represented as a sequence for discovering motifs or rules, for clustering or classifying sequences, or for answering query by contents.

DTW finds the optimal alignment between two given time series. One of its strength is to cope with local distortions along the time dimension. DTW, however, is slow and costly to calculate as its time complexity is quadratic. Therefore, three families of approaches have been defined to accelerate DTWbased searches.

Lower-Bounding Functions

The first family of approaches includes cheap-to-compute functions lower bounding DTW [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF][START_REF] Yi | Efficient retrieval of similar time sequences under time warping[END_REF]. They accelerate the DTW process because they quickly prune sequences that could not possibly be a best match. Then, the full DTW is computed on the list of candidate sequences that remain. Lower bound functions have two key properties: (i) They incur no false dismissals as suggested in [START_REF] Faloutsos | Fast subsequence matching in time-series databases[END_REF] and (ii) the tighter they are, the more they prune. At a very fine level, even LB Keogh, known to be the tightest lower bound in the literature, can be quite far from the true value of the DTW [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF]. There is therefore room for improvement, as a tighter bound would allow more pruning and faster searches. Yet, designing a tighter lower bound function is difficult as its computation is likely to become complex, and thus less profitable.

The first contribution of this paper is the definition of approximate lower bounding functions built on top of the popular LB Keogh and LB PAA (exact) lower bounds. These approximate lower bounding functions have a user-tunable tightness: It can range from the one achieved by their exact counterpart, and in this case, no false dismissals are observed, to much tighter bounds where false dismissals are possible as the DTW may get over estimated.

Indexing Time-Series

The second family of approaches includes schemes specifically designed to cope with very large databases of time series. Camerra, Palpanas, [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF], [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] and Vlachos, Hadjieleftheriou, Gunopulos and [START_REF] Vlachos | Indexing multidimensional time-series with support for multiple distance measures[END_REF] all index time series for improved performance. Indexing shrinks the search space to some neighborhood defined around the query sequence, which, in turn, dramatically reduces the response time. Not surprisingly, lower bounding functions are at the core of such indexing schemes for even better performance. For example, [START_REF] Lin | Experiencing SAX: a novel symbolic representation of time series[END_REF] define iSAX MinDist that prunes sequences based on the range intervals determined for quantizing time series.

The second contribution of this paper is the definition of an approximate lower bounding function used in iSAX. It builds on iSAX MinDist, and its tightness can also be tuned by the users to offer a wide range of response time versus accuracy trade-offs.

The performance of iSAX depends on the one hand on the lower bounding function used, and on the other, it depends on the way the tree used for indexing is built. iSAX and iSAX2.0 quantize time series in a tree of symbols. As discussed in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], the quantization intervals are determined in order to be as balanced as possible assuming the data in time series are normally distributed. This is, in fact, rarely the case.

The third contribution of this paper is the use of a k-means-based quantization step for iSAX. Using k-means avoids to make any strong assumption on the distribution of data in time series. The resulting index better fits data. Some special care is taken to balance the tree of symbols, as this is key to performance.

Approximate Searches

The third family has approaches that have traded accuracy for response time [START_REF] Chu | Iterative deepening dynamic time warping for time series[END_REF][START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF]. In this case, great performance gains are obtained at the cost of some potential false dismissals. What is central to these schemes is the metrics they use for assessing the similarity of time series, metrics that are no longer lower bounding the DTW. In addition to its exact iSAX MinDist, iSAX defines such an approximate metrics. [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF] show the approximate version of iSAX can answer a query in less than a second, while exact iSAX requires about 10 min.

The last contribution of this paper is to show how approximate lower bounds can be used to boost the performance of the approximate version of iSAX.

Discussion

Experiments made with state-of-the-art datasets [START_REF] Keogh | The ucr time series classification/clustering homepage[END_REF] show that these contributions significantly accelerate DTW. When considering small datasets, sequentially scanning the data collection is still appropriate. In this case, using approximate lower bounding functions is an order of magnitude faster than using exact lower bounds while indeed very rarely triggering false dismissals. It is twice as fast as Iterative deepening dynamic time warping (IDDTW) [START_REF] Chu | Iterative deepening dynamic time warping for time series[END_REF] while achieving comparable recall performance.

Performance gains are even more significant when considering large datasets for which indexing is needed. Approximating iSAX MinDist allows to run iSAX two order of magnitude faster than its exact version while still returning one of the 50 actual nearest neighbors at rank 1 in 87% of the cases. Performance gains are about one order of magnitude compared to the approximate version of iSAX while achieving comparable retrieval performance, quality-wise. These approximate methods can deal with applications for which exact searches would be far too long to perform.

This paper is structured as follows. Section 2 gives the necessary background on the DTW and on few popular techniques boosting its performance. Section 3 details how these techniques can be improved thanks to approximate lower bounding functions. Section 4 introduces iSAX + , an indexing tree borrowing from iSAX that, however, relies on k-means clustering for determining quantization intervals. Note these two sections give pseudocodes. Section 5 presents the extensive experiments made using the datasets of [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] and [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF], showing the algorithms we propose here outperform state-of-the art solutions. Section 6 concludes.

Accelerating DTW

This section reviews the background material needed for this paper. We first give a quick description of the DTW process. We then describe four very classical methods that have been designed to accelerate its computation. We then discuss IDDTW [START_REF] Chu | Iterative deepening dynamic time warping for time series[END_REF] which approximates DTW searches.

DTW

It has been demonstrated that similarity searches in time series is best when using DTW. DTW takes two time series, Q and C, of length m and l respectively, where

Q = q 1 q 2 . . . q m ; C = c 1 c 2 . . . c l (1)
To evaluate the similarity between Q and C, a m × l matrix S is computed such that: ∀(i, j) ∈ {1, . . . , m} × {1, . . . , l}, S i,j = d(q i , c j )

(2)

where d(q i , c j ) is the distance between q i and c j . A warping path W is the set of K elements of S forming a path W = w 1 . . . w K where w k = (i, j) ∈ {1, . . . , m} × {1, . . . , l}. W has to meet three conditions: It has to be continuous, monotonic and has to meet boundary conditions (it goes from (q 1 , c 1 ) to (q m , c l )). Dynamic programming techniques are used to efficiently find this path via a recurrence over a cumulative distance γ i,j defined as follows:

∀(i, j) ∈ {1, . . . , m} × {1, . . . , l}, γ i,j = min S i,j + γ i-1,j-1 S i,j + γ i-1,j S i,j + γ i,j-1 , (3)
A direct implementation of DTW has a time and space complexity of O (m × l).

It has been shown that for some applications, restricting the admissible paths around the diagonal could be beneficial in terms of both complexity and relevance of the resulting metrics. The most commonly used constraints are the Sakoe-Chiba band [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] and the Itakura parallelogram [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF]. The former constrains the admissible paths to lay inside a band of fixed width around the diagonal path, as shown in Fig. 1, while the latter uses a parallelogram that has the diagonal path as one of its diagonals.

Accelerator #1: LB Keogh

LB Keogh [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] is probably the best-known lower bound function. LB Keogh(Q, C) first rescales Q and C to the same length n, then builds the bounding envelope of Q by determining all its maximum and minimum values inside a Sakoe-Chiba band [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] (or an Itakura parallelogram [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF]). It then sums the distances from every part of C not falling within the bounding envelope of Q to the nearest point from that envelope having the exact same time stamp. [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] show, for any Q and C rescaled that LB Keogh ≤ DTW. The complexity of LB Keogh is linear in n.

Accelerator #2: LB PAA

LB PAA is another approach further reducing the cost of computing LB Keogh. LB PAA first requires to chop all sequences C stored in the database into N chunks, containing a number of points that depends on the length of considered sequences. This creates sequences downsampled along their timeline. Downsampled versions of C are then compared to a downsampled version of the envelope of Q in a manner that is similar to the one for LB Keogh. This gives the LB PAA function, lower bounding LB Keogh, and, therefore, DTW. Note the complexity of LB PAA is equivalent to the complexity of computing LB Keogh for time series of length N , which is smaller than n and thus computing LB PAA is more efficient than computing LB Keogh.

Accelerator #3: iSAX

Both LB Keogh and LB PAA must be computed for each sequence C from the database to find the one that is the most similar to Q. Researchers have thus worked on finding indexing techniques where only a small subset of relevant sequences from the database would have to be compared to Q, which would in turn prune the search space even more.

Since LB PAA turns sequences into fixed-length descriptions lying in a Ndimensional space, it is possible to index them, typically in a R-Tree, as proposed in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF]. The index eventually contains minimum bounding rectangles (MBR) in which similar sequences are grouped. At query time, the algorithm finds the sequences from the database that are close to Q using the proximity of their respective MBR, quickly determined through the index.

Building on this, an even more efficient DTW indexing scheme called iSAX was proposed in [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF]. iSAX also downsamples sequences, but it additionally quantizes them. iSAX thus gives a discrete representation to continuous sequences. It turns time sequences into a series of symbols that can easily be indexed in a tree. The quantization boundaries giving symbols for the values are determined according to a standard normal distribution. The size of the alphabet is dynamically adjusted during index construction to balance symbols. iSAX was later extended into iSAX2.0 [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], that achieves more efficient index construction thanks to bulk loading and more balanced clusters via a smart symbol selection scheme to use when splitting a node from the tree.

In this case, the iSAX MinDist metrics defined in [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF] that is used to assess similarity between a query and a node of the tree is in the same spirit as the ones presented above, except that the ranges assigned to nodes are used to replace MBRs. The search process then consists of visiting each leaf of the tree in ascending order of their iSAX MinDist values, and the algorithm stops as soon as all remaining leaves in the tree have iSAX MinDist values greater than the current best distance found.

Accelerator #4: Approximate iSAX

Shieh and [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF] also present a slight variant of the exact iSAX version that uses iSAX MinDist. This variant returns approximate results as it does not scan all theoretically possible leaves, but solely the one leaf with the lowest lower bound. That single leaf is entirely scanned, and LB Keogh is determined for each sequence, which possibly prunes some sequences from any further analysis. As usual, full DTW are computed for the remaining sequences from that leaf to build the final result returned to the user. Returning data from that single leaf is of course very fast. Searching the index for time series that are similar to the query thus aborts as soon as the most probable leaf from the tree has been scanned. For the sake of clarity, we call this way of accelerating the searches iSAX Approx.

IDDTW

Iterative deepening dynamic time warping (IDDTW) [START_REF] Chu | Iterative deepening dynamic time warping for time series[END_REF] aims, as previously introduced lower bounds for DTW, at computing DTW only for those candidate sequences in the database that are likely to be similar to the considered query. This likelihood is estimated by computing DTW between downsampled versions of the sequences. The principle is then to start comparing sequences at a very coarse resolution and keep refining the latter until the sequence can be discarded from the candidate list or, if not, until full resolution is reached. To test if a sequence can be discarded at a given resolution, a learned distribution of the estimation error at this resolution is used. Given this error distribution, if the probability that the computed estimation can lead to a lower DTW value than that of the best sequence so far is lower than a user-defined threshold, the sequence is discarded. If not, this process is repeated at a finer resolution.

This algorithm heavily relies on probabilistic estimation. Therefore, the false dismissal assumption is discarded for the sake of lower computational cost.

Improving accelerators via approximations

This section describes how the efficiency of the four DTW accelerators can be improved thanks to approximations. It first describes how the A Keogh and the A PAA approximate lower bound functions can be defined1 from their exact counterparts LB Keogh and LB PAA. This improves the performance of the accelerators #1 and #2. Then, this section moves to defining A iSAX which improves the performance of the accelerator #3, indexing DTW via iSAX. It also defines A iSAX Approx, built from iSAX Approx which enhances the performance of Approximate iSAX. This section ends by presenting the pseudo-codes for two algorithms derived from that of [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] and [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF], one for a sequential-scan-based search algorithm, the other for an index-based search.

A Keogh

The pruning power of LB Keogh is directly linked to its tightness, that is the ratio of lower bound value to actual DTW value. The tighter it is, i.e., the closer to the real value of the DTW the lower bounding values are, the better. Given a time series, computing its LB Keogh value does not help in knowing whether it is tight, i.e., close or far from its actual DTW value. It is not the value of LB Keogh that matters, but the relative difference between that value and the true value of DTW. It has been observed that LB Keogh is often quite smaller than the true value of the DTW, especially with rapidly varying sequences. What is key is to have a way to distinguish the case where (i) the value of the lower bounding function is small because the DTW is also small, from the case where (ii) it is small because its calculation is poorly tight.

To reach that goal, we first explain why applying a single multiplicative factor to the obtained lower bound is not enough. Given that, we describe the design of a cheap-to-compute upper bounding function that is very similar in spirit to LB Keogh. This upper bound is used to get an indication on the tightness of the lower bounding function LB Keogh. By repeatedly computing the lower and upper bounds on a first set of sequences (denoted in the following as the learning set), it is possible to learn and then model the overall tightness. We then introduce a mechanism to determine, given a tunable parameter that is closely linked to the tightness model, how many DTW values are likely to be overestimated. Overall, this yields approximate lower bounding functions. Thanks to the model for overestimations, it is possible to fine tune the pruning power of 1. Observed median of T values and its quartiles when the input parameter is set so as to imitate a magic lower bound 99% tight. Note that for αLB Keogh and A Keogh, learning values are given between parenthesis as their parameters are set so as to get 99% median tightness on the learning set.

such functions: very aggressive functions that are likely to overestimate DTW will have a very good pruning power, but may in turn create many false dismissals, in contrast to milder functions.

αLB Keogh

As explained above, it has been observed that LB Keogh is often quite smaller than the true value of the DTW, especially when dealing with rapidly varying sequences. A first option to approximate existing lower bounds would have been to learn the typical ratio between DTW values and lower bound ones using a large set of diverse sequences and a set of representative queries. It would then be sufficient to compute a multiplication factor α to apply to LB Keogh during the search process in order to increase the tightness of a new approximate lower bounding function that we call αLB Keogh.

We show that this idea leads to poor approximations through experimentation. Since we believe that at large scale, a database of sequences is highly diverse, we merged into a single database all the datasets listed in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] (that will later be referred to as dataset DS 1 ) and set parameter α in order to achieve 99% tightness on the learning set.

The first two rows of Table 1 show that the median value for tightness is around 92% for LB Keogh, meaning that the returned value is, in median, 8% lower than that of the DTW. Note also that αLB Keogh tends to overestimate DTW, its median tightness being around 105%. Moreover, interquartiles intervals are large for these two methods, which implies, for LB Keogh, a high number of sequences for which the DTW is underestimated and, for αLB Keogh, a high number of sequences for which the DTW is overestimated. This shows that a single multiplicative factor is not sufficient to build a reliable approximative lower bound. Indeed, using αLB Keogh would only slightly increase tightness when the ratio between DTW and LB Keogh is large, while it will return a value that would be too high when the ratio tends towards 1. [START_REF] Vlachos | Indexing multidimensional time-series with support for multiple distance measures[END_REF] and [START_REF] Kashyap | Scalable knn search on vertically stored time series[END_REF] propose to use an upper bound to prune sequences that could not match a given query. In this paper, we are using an upper bound not for pruning, but to determine the interval (between the lower and upper bounds) within which the true value of the DTW is. If the interval is small, then it is likely the lower and upper bounds are tight.

UB Keogh: an Envelope-based Upper Bound

One straightforward way to upper bound the DTW is to compute the Man- hattan distance2 between Q and C. A more elegant way is to rely on the envelopes around sequences that are already used for computing lower bounds. We experimentally observed that the intervals around the true values of the DTW are more stable when the bounds are defined according to the same underlying principle (envelopes). In addition, as it will be detailed later in this paper, relying on envelopes facilitates the computation of an upper bound on downsampled data as it is the case for LB PAA. Building on LB Keogh, it is possible to define an envelope-based upper bound, called UB Keogh. UB Keogh is computed by summing the distances between every point in the query and the furthest point having the exact same time stamp on the envelope of each database sequence. The mathematical definition of UB Keogh as well as the proof it truly upper bounds DTW is provided in the Appendix. By definition:

U Q C L (a) LB Keogh(Q, C) L (b) UB Keogh(Q, C)
LB Keogh ≤ DTW ≤ UB Keogh. ( 4
)
The complexity of computing UB Keogh is the same as the one of computing LB Keogh, i.e., it is linear in n. When processing Q and C, the difference between LB Keogh and UB Keogh gives a clear indication on the tightness of these bounds, a small difference suggesting that the bounds are tight.

Modeling Tightness

It is possible to learn what the values for LB Keogh and UB Keogh are, as well as the ones for the true DTW for each pair of sequences in the learning set. A histogram of the distribution of these values, once normalized, can then be computed. Normalized values are defined as follows:

s-DTW(Q,C) = DTW(Q, C) -LB Keogh(Q,C) UB Keogh(Q,C) -LB Keogh(Q,C) . (5)
s-DTW is closely related to the tightness observed on the learning set: A left-shifted curve for the probability density function of s-DTW suggests most LB Keogh values are tight, while most are loose when right-shifted. We validate this idea through experimentation. Since we believe that at large scale, a database of sequences is highly diverse, we merged into a single database all the datasets listed in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] (that will later be referred to as dataset DS 1 ).

Figure 3a shows, for DS 1 , the distribution of s-DTW over all sequences. It also plots distributions observed for two sub-datasets. Few comments are in order. First, all plots are away from the x-axis zero value, meaning that for every pair of sequences considered, LB Keogh is strictly lower bounding DTW. This is basically the proof that there is room from improvement: it is possible to use a greater value than LB Keogh (i.e. shifting decision value to the right) while still experiencing very few (or even no) DTW overestimations (shown on the y-axis of Fig. 3b). Second, there is a significant shift between the curves, which motivates the use of a learnt histogram modeling the s-DTWdistribution rather than a single shift value.

The distribution of s-DTW values is an indicator of the tightness observed on the learning set: We next show how it can be used as a model, to predict the tightness of sequences outside the learning set. Ideal case N learn =10 N learn =30 N learn =1,000 Fig. 4. Setting tightness as an input parameter. Here, a given tightness value is targeted using dichotomy and we study the behavior of A Keogh with respect to the number of learning sequences. The dataset used for this experiment is the one called "Two Patterns" from the larger dataset that will later be denoted as DS 2 .

Determining the Over-Estimation Rate

By computing the cumulative distribution of s-DTW, one can estimate the expected distance over-estimation rate. Assume one wants to speed up the search process at the price of potential false dismissals. That user decides to allow for p =10% of over-estimated distances. Given the learned cumulative distribution, one can find the value of β on the x-axis that gives p on the y-axis (see the dashed lines on Fig. 3b). β is then used to compute the approximate lower bound equal to

A Keogh = LB Keogh + β(UB Keogh -LB Keogh) (6)
that will subsequently be used to prune sequences. p might be meaningless for some users. Instead, they might prefer achieving a target tightness T . Finding the optimal value for p in order to achieve a target tightness T can be performed by using dichotomy, as explained in more details in Sect. 3.1.5. The last row of Table 1 provides indications on the observed tightness of A Keogh when setting p in order to target 99% tightness. The observed median tightness is very close to this target, and its interquartile interval is small, suggesting a robust estimation that backs the use of information from upper bounding functions to get a better estimation of the actual DTW value.

Discussion

In this section, we will discuss some of the salient properties A Keogh has, making it special with respect to other lower bounding functions.

Learning from data To achieve better tightness, A Keogh relies on the distribution of s-DTW values in the [0; 1] interval. An estimate of this distribution is obtained by computing a histogram of such values on a learning set. Note that we do not aim at gathering general knowledge from the data, but rather dataset-specific information. Hence, one can either use a learning set if there is one available or build one by using a subpart of the dataset to be searched.

Cost of off-line learning Building such a histogram requires off-line DTW computations. More precisely, a learning set made of N learn sequences implies

N learn (N learn -1) 2
DTW computations that have to be done once and for all at the learning step. This fixed off-line cost is affordable and helps saving time for the critical phase during which user queries have to be processed in a timely fashion. Experimental results presented in this paper do not refer to this cost as we believe the on-line query cost is the key for such methods. To illustrate that few learning sequences are required for our method to work properly, we show in Fig. 4 the difference between expected and obtained tightness. In this case, a learning set made of 30 sequences is sufficient to achieve a good modeling of tightness since the bias that can be observed in this case is still present when using 1,000 learning sequences, which means it is due to the bias between the learning and test sets. In this case, then, 30×29 2 = 435 off-line DTW computations are sufficient to learn characteristics of a dataset that is made of 1,000 learning sequences.

Moreover, the number of sequences used for learning drives the resolution of parameter p, as the latter is determined by the resolution of the cumulative density function learned on these sequences. A reasonable condition is that the given p value corresponds to at least one observation in the training set:

p • N learn (N learn -1) 2 ≥ 1, (7)
which is satisfied for

N learn ≥ 1 2 + 1 4 + 2 p . (8)
Using a typical value of p = 0.01, this condition becomes N learn ≥ 15, which is satisfied by the learning set used in this section.

Setting an optimal value for p Setting an optimal value for parameter p depends on the kind of problem the user wants to tackle: When dealing with relatively small datasets, she might want to achieve very good accuracy results, as the query cost is not a strong issue in this case, whereas when dealing with very large datasets, setting a large value for p might be the only way to get a result in an acceptable amount of time. Finally, note that, when setting p = 0, A Keogh=LB Keogh.

Setting a constraint on T As stated above, it might be more meaningful for a user to set an expected value for T rather than for p. This can be done by using dichotomy as the function of T giving p is nondecreasing. Indeed, tightness is defined as the ratio of A Keogh to DTW, and is then a nondecreasing function of β (cf Equation 6). As β is a nondecreasing function of p (this comes from the definition of a percentile), we get that tightness is a nondecreasing function of p. Dichotomy consists in recursively splitting the [0,1] p value interval into smaller ones. To do so, at each step, the current interval is divided into two parts, and the middle of the interval is used as a p value. If the tightness computed on the learning set for this p value is lower than the targeted one, then the upper half interval is kept for the next iteration; otherwise, the lower half interval is considered. This process is performed on the learning set so as to learn a lookup table that can be used afterward for the user to give T as an input.

Algorithm 1 Sequential scan based searching algorithm.

Require: Q[1..n], DB[1..n db ][1..n], β (L, U) ← get envelope(Q) for i = 1..n db do lb ← LB Keogh(L, U, DB[i]) ub ← UB Keogh(L, U, DB[i]) A Keogh[i] ← lb + β (ub -lb) end for argsorted A ← argsort(A Keogh) best so far ← ∞ best so far idx ← -1 for i in argsorted A do if best so far ≤ A Keogh[i] then break end if if DTW(Q, DB[i]) < best so far then best so far ← DTW(Q, DB[i]
) best so far idx ← i end if end for return (best so far, best so far idx)

The effectiveness of this process is illustrated in Fig. 4, which shows that this dichotomy process is efficient even if using very few learning sequences.

If the user prefers to set a complexity constraint in terms of the expected amount of DTW computations, it is then possible to set the value for p by numerically inverting the relationship between this amount of computations and p. The proof that the amount of DTW computations is a nondecreasing function of p comes from the fact that this amount is a nondecreasing function of the tightness that has been shown earlier in this section to be a nondecreasing function of p.

A PAA

The previous section defined UB Keogh and explained how to calculate A Keogh. It is possible to obtain UB PAA and A PAA in a very similar manner. Few comments are in order, however. First, LB PAA downsamples all sequences into N chunks. The envelope of Q is therefore defined according to its downsampled version. UB PAA is thus also defined on this downsampled version. Mathematical definition of UB PAA is in the Appendix. Second, LB PAA ≤ DTW ≤ UB PAA. Last, the formula for modeling tightness (see Eq. ( 5)) has to be slightly changed in a trivial manner to use LB PAA and UB PAA. Note also that when the rate of DTW over estimations p is set to 0, then A PAA=LB PAA.

A Keogh and A PAA for Sequential Scan

LB Keogh and LB PAA are defined in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF]. They are both used together with a search algorithm doing a sequential scan of the entire collection. Either lower bounds are computed on every sequence C from the database to quickly get it compared to Q.

This section presents a slight variant of this sequential scan search process that uses either LB Keogh, UB Keogh and A Keogh or LB PAA, UB PAA and A PAA. The pseudo-code in this algorithm focuses on A Keogh. It is straightforward to adapt to A PAA.

Prerequisites n db sequences are stored in a database DB. Each sequence has been rescaled to length n.

Off-Line Learning For each pair of sequences (C i , C j ) from the learning set, DTW, LB Keogh, UB Keogh and s-DTW are computed. Then, the average cumulative distribution for s-DTW is calculated, and a lookup table C keeping track of the quantiles of this distribution is built.

Query Time User Input The user inputs are Q and p. Q is the query, and it gets rescaled to length n. p is the average overestimation rate. p is used to probe C, resulting in the target quantile β.

On-Line Searching The sequential-scan-based searching algorithm is detailed in Algorithm 1. An important point here is that when setting β > 0, the break instruction can be triggered before finding the exact nearest neighbor.

A iSAX and iSAX indexing

Indexing sequences is needed at large scale. The most popular time series indexing scheme is probably iSAX [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF]. For efficiency, iSAX uses a lower bound function called iSAX MinDist. Following the principles defined above, it is possible to create an approximate lower bound function A iSAX built on top of iSAX MinDist.

The definition of iSAX MinDist is slightly more complex compared to LB Keogh or LB PAA as it relies on the intervals defined during the index tree construction by the quantization process. With iSAX it is possible that some intervals associated with nodes in the index tree have infinite bounds. In this case, computing an upper bound is problematic. To facilitate this computation, we had to slightly modify the tree building process by storing, for each node that has at least one infinite bound, the corresponding extremum value for the sequences that are stored in the node or one of its children nodes. Therefore, the upper bound will use the actual extremum value instead of the infinite bound in its formula. That upper bound is called iSAX MaxDist and its mathematical definition is in the Appendix. Then, iSAX MinDist and iSAX MaxDist are used seamlessly in Eq. ( 5), eventually defining A iSAX.

The resulting iSAX index-based searching algorithm is detailed in Algorithm 2. Note this algorithm uses iSAX MinDist, iSAX MaxDist and A iSAX as well as the sequential scan algorithm presented above (and thus is relies on A Keogh).

Prerequisites n db sequences are stored in a database DB. Each has been downsampled to length N .

Off-Line Indexing iSAX 2.0 off-line indexes all sequences as in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF].

Off-Line Learning Same as in Sect. 3.3. Additionally, however, quantiles for A iSAX and A Keogh must be learned as both lower bounding functions are used at query time.

User Input The user inputs Q and p. The two quantiles β LB Keogh and β iSAX MinDist are obtained from p and C.

On-Line Searching The index based searching algorithm is detailed in Algorithm 2. It is built on iSAX (see [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF]), with specifics due to the use of approximate lower bounds. Note that, isax approximate search, called to initialize the process, returns the leaf node having a SAX signature that matches the one of the query3 and seq scan is Algorithm 1.

iSAX Approx indexing

Shieh and [START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF] describe a variant of iSAX that returns approximate results. We have presented this variant in Sect. 2.5 and called it iSAX Approx. iSAX Approx determines the one best iSAX-leaf, computes LB Keogh on all sequences from that leaf for pruning, and then computes full DTW on the remaining sequences. As it uses LB Keogh, it is possible to patch this algorithm such that it uses A Keogh instead.

It is therefore quite easy to turn iSAX Approx into A iSAX Approx. iSAX Approx is already an approximate search scheme, due to its leaf picking strategy. A iSAX Approx somehow piles up another layer of approximation as it uses A Keogh. Using A Keogh at the heart of iSAX Approx typically divides by 7 the number of sequences that go through a full DTW process for comparable recall performances. The pseudo-code for this A iSAX Approx strategy is simply made of the four first lines of Algorithm 2 as it directly returns (idx bsf, dist bsf) once the best iSAX node scanned by Algorithm 1.

iSAX + : Quantizing with k-means

iSAX is a very successful approach because it turns continuous time series into discrete sequences of symbols, for which efficient high-dimensional indexing schemes exist. The cornerstone of iSAX is therefore its quantization process which is clearly defined in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF]. It makes one strong assumption: The data in time series before quantization is assumed to be normally distributed. From that assumption, iSAX determines the quantization intervals Algorithm 2 Index based searching algorithm for iSAX. to balance the probability with which symbols will appear in the quantized sequences. [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF] show it is key for performance to have (roughly) equiprobable symbols.

Require: Q[1..n], β LB Keogh , β iSAX MinDist (L, U) ← get envelope(Q) node bsf ← isax approximate search(Q) (idx
Unfortunately, the data in real-world time series are unlikely to follow a normal distribution. Note the biggest datasets used in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF] are synthetic and have been created so that considered time series have marginal normal distribution. As real-world time series rarely stick to a normal distribution, quantifying them as iSAX does is unlikely to give equiprobable symbols, which, in turn, hurts performance.

It is possible to relax this normal distribution assumption by using a variant of the k-means approach. k-means is a widely used unsupervised classification algorithm that takes as input a set of points x and a number k of clusters to build. It makes no assumption on the distribution of data. It tends to minimize intra-cluster variance defined as follows:

V = k i=1 x∈Ci x -c i 2 (9)
where c i is the centroid of cluster C i .

Algorithms belonging to the k-means family have been used as unstructured quantizers in the context of image and video search [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF][START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF][START_REF] Paulevé | Locality sensitive hashing: A comparison of hash function types and querying mechanisms[END_REF][START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF]. In this context, k-means has proved to nicely fit the data in space. This is the rationale for using k-means for quantizing time series. k-means, however, is also known to create Voronoi cells having rather uneven cardinalities.

The variant we use here not only does k-means but also balances the clusters it creates. We now detail the way k-means can be used for iSAX as well as the cluster balancing strategy.

k-means for Building the iSAX-Tree

The original iSAX tree forms a hierarchy of nodes. We therefore use k-means in a hierarchical way [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] to get a similarly shaped index tree. The root node of the original iSAX index tree has b w children, where b is the base cardinality of symbols and w is the word length, i.e., the number of symbols in a sequence. Therefore, the root node when using the k-means approach has the same number of children nodes and results from determining k = b w centroids. Other levels below in the tree are made by running k-means with k = 2 as for the original iSAX which has also 2 child nodes per parent node.

As it is the case for the original iSAX approach, building an index tree with k-means requires to split a node into two child nodes once the parent node gets fully filled with sequences. This very traditional recursive node splitting process eventually creates leaf nodes containing at most th sequences. Not surprisingly, a node keeps track of its MBR and its associated centroid, that centroid being produced by k-means. Inserting a new sequence in the index tree is performed by going down the tree and, at each level, by determining the centroid that is the closest to the sequence to insert. If that sequence has to be inserted in an entirely filled node, then a split occurs. k-means is applied to the sequences in that node that are then moved to child nodes according to the centroids newly determined.

Balancing the k-means-based iSAX-Tree

It has been shown in [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] that a k-means computed over a large set of high dimensional features produces clusters having quite different cardinalities. This, in turn, increases the response time variance of the algorithms processing the data in clusters, and performance might be particularly bad when using highly populated clusters. Having balanced clusters reduces this variance and improves performance, overall. Note this is also why the original version of iSAX creates a balanced index tree. [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] balance the clusters produced by a k-means as follows. It first projects all the data points as well as the centroids computed by the kmeans into a higher dimensional space. Then, it iteratively enlarges the distances between the centroid of a cluster and the corresponding data points in that cluster in proportion to the cardinality of the cluster. This moves inward the frontiers of highly populated clusters and reduces their cardinality because it assigns some of the points from that cluster to other, less crowded, clusters. [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] can be run until all clusters contain the same number of points, or it can be stopped earlier, after a given number of iterations or when the standard deviation computed over the cardinalities of clusters falls below a threshold.

The mathematical formulas used in [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] can be significantly simplified when the algorithm has to balance the cardinality of k = 2 clusters. In this case, they can be turned into a simple algebraic equation given in the Appendix.

iSAX +

This section presents the pseudo-code for the complete iSAX + indexing scheme (Algorithm 3). It borrows a lot from iSAX. It makes use of k-means and of the procedure described in [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] to create balanced clusters as described above.

The build tree function is in charge of splitting overfilled nodes upon insertion. It uses the th parameter defined for the original iSAX. At splitting time, the function split is invoked. It first quantizes the data using k-means. Then, it balances the children nodes using the mechanisms presented in [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF] and briefly sketched above.

Note the pseudo-code for querying the tree is not presented here as it is identical to Algorithm 2. The only (minor) difference lies in the approximate search function that relies on the proximity of the centroids obtained by the k-means (and not on the set of symbols as for the original iSAX).

Experiments

This section presents the extensive performance evaluation showing the techniques proposed in this paper outperform state-of-the art solutions. The first set of experiments compares LB Keogh and A Keogh, as LB Keogh is the best exact lower bound ever proposed. This comparison involved evaluating their respective tightness and their resulting effectiveness in pruning more sequences. We also evaluate the impact of relaxing the no false dismissal assumption in the case of A Keogh on the quality of the results returned to the user. Note this involves checking the ranking of candidate sequences that will go through a full DTW process as this ranking might differs from the one determined by LB Keogh. Note LB PAA and A PAA are also compared.

The second set of experiments compares A Keogh with IDDTW. IDDTW also approximates DTW, but it uses downsampling instead of an approximate lower bound.

The third set of experiments focuses on the impact of the above techniques on iSAX. It shows its performance is improved when using the approximate lower bounds A iSAX and A iSAX Approx.

Finally, we give the performance of iSAX + that relies on a balanced k-means quantization step.

For the sake of reproducibility of experiments, all evaluations are conducted using publicly available [START_REF] Keogh | The ucr time series classification/clustering homepage[END_REF] and widely used datasets [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF][START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF]. In all experiments, the baseline similarity measure is a DTW constrained by a Sakoe-Chiba band with a size set to 10% of the length of sequences, as suggested in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF].

Algorithm 3 Basic operations for iSAX + tree.

function build tree(DB[1..n db ][1.

.n], th, b, w, α, r) root ← new node() root.insert (DB) while there exists a node v such that card

(v) > th do if is root (v) then split (v, α, r, b w ) else split (v, α, r, 2) end if end while return root end function function split(v, α, r, k) centroids ← k-means(v.stored data, k) (centroids, assign) ← balance (centroids, v, k, α, r) for c in centroids do v c ← new node (centroid = c) v c .stored data ← assign[c] v c .MBR ← MBR (v c .stored data) v.children.add (v c ) end for end function function insert(v, S, α, r) if has children (v) then v 0 = closest child (v, S) insert (v 0 , S) else if card(v) == th then split (v, α, r, 2) v = closest child (v, S) end if v.stored data.

add (S) end if end function

The experiments use two datasets. The first one, referred to as DS 1 , is precisely defined in [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF]: It is made of 32 time series, 50 subsequences of length 256 being randomly extracted from each dataset to build the test set. In order to learn our lookup table C, we built our own learning set per time series using 100 randomly extracted subsequences of length 256.

The second dataset, that we call DS 2 , is available at [START_REF] Keogh | The ucr time series classification/clustering homepage[END_REF] and was used in several publications. It is made of 20 sub-datasets, each split in advance into a learning and a test set. When running experiments involving DS 2 we naturally use these ready-made learning sets to build C. Note that, these sub-datasets also provide classification labels for every sequence. When evaluating classification, we assign to the query the label of the most similar sequence returned by the search algorithm.

For datasets DS 1 and DS 2 and for both A Keogh and A PAA, we report the first-nearest-neighbor (1-NN) retrieval rate, that is, the ratio of queries for which their true nearest neighbor (according to the DTW) is ranked first.

Note that, when necessary, estimator medians are reported together with their inter-quartile interval.

Compared Tightness of Lower Bounds

Figure 5 compares the tightness achieved by LB Keogh and A Keogh on datasets DS 1 and DS 2 . It can be seen that the observed tightness for LB Keogh varies significantly between series: For example, while it is tight for series #3 or #10 with DS 2 , it is very loose for series #8 or #13. In contrast, A Keogh is overall much closer to 1 for both datasets. Note that, when p = 0.05 or p = 0.1, the improved tightness is achieved with high precision, as shown by the inter-quartile intervals that do not cross the 100% tightness limit for any dataset. As predicted, setting higher values for p leads to higher achieved tightness. It is important to see that the results are more stable when using large learning datasets as it is the case of DS 2 , where inter-quartile intervals are so small that they can hardly be seen on figures.

Similarity of Candidate Lists

Algorithm 1 sorts all sequences from the database according to the value of A Keogh. That order depends on the value given to p and might might possibly differ from the one that is determined when sequences are ordered according to their LB Keogh value. Moreover, with A Keogh, when p > 0, the break instruc- tion can be fired before identifying the best match, causing a false dismissal. It is therefore key to check if the order according to which sequences are sorted along their lower bound is somehow similar when considering DTW and LB Keogh on the one hand or DTW and A Keogh on the other. If the latter order is quite similar, then it is likely that the best match will also be returned when using A Keogh.

The Spearman correlation coefficient ρ is a widely used metric for assessing the similarity of two ordered lists. The closer to 1 ρ is, the more similar the lists are. It is defined as:

ρ = 1 - 6 n db i=1 (r 1 i -r 2 i ) 2 n db (n 2 db -1) , (10)
where r 1 i and r 2 i are the ranks of sequence i in both lists. Table 2 gives the median values for ρ computed on datasets DS 1 and DS 2 . In this experiment, the set of r 1 i has been produced by ordering the n db sequences along their DTW value; r 2 i is either LB Keogh, A Keogh, LB PAA or A PAA. In all cases, the values for ρ do not vary much when p ranges from 0 to 0.1. Note also the large overlap between inter-quartile intervals suggests the observed differences for ρ between the three versions of A Keogh (resp. A PAA) are not significant. Finally, it is interesting to notice that A PAA has higher ρ values than LB PAA. In other words, when a lower bound is far from tight, it can be useful to grasp information from the upper bound so as to generate ordered candidate lists that are closer to the ones that DTW would get.

Another important point is that the observed ρ values are much smaller for LB PAA (resp. A PAA) than they are for LB Keogh (resp. A Keogh). This implies that smaller values should be used for p when trying to approximate LB PAA as each DTW overestimation is more likely to induce a false dismissal.

Retrieval performances

We used DS 2 to evaluate the performances of our proposed approximate lower bounds in terms of retrieval performances as well as computational cost. Re- trieval is expressed by both the 1-NN retrieval rate and the correct classification rate. The computational cost is expressed as a ratio of the number of DTW computations induced by the approximate lower bound to the number of DTW computations induced by its corresponding exact lower bound. Note that, this cost does not include our off-line learning step, as it aims at reflecting on-line query cost. Table 3 presents these results for A Keogh and A PAA. The first thing to notice is that in all cases, the computational cost decreases considerably even for small values of p. Moreover, as previously suggested, small values of p should be used when considering LB PAA to reduce the likelihood of false dismissals. The classification rate drops much slower than does the 1-NN retrieval rate, showing that when the true nearest neighbor is not ranked first, another sequence of the same class is, which is key to performance, quality-wise. For more detailed information, refer to Table 4 that gives the cost versus accuracy trade-off for all sub-datasets in DS 2 .

Comparison with IDDTW

Both IDDTW and A Keogh approaches rely on modeling the expected error that is induced by approximation. We therefore compare both approaches in terms of accuracy and computational cost, using DS 1 . Here, the computational cost is determined slightly differently from what was done previously, as IDDTW is an iterative process running DTW computations at different scales. We therefore define this cost as the total number of elementary distance computations that are needed for searching the database. For example, if a database of n db sequences of length n was searched using DTW with no temporal constraint, the cost would be n db × n 2 , since n db DTW would be evaluated, each of which would use n 2 distance computations. This amount is then divided by the complexity of the reference method that is a full scan of the database using DTW constrained to a Sakoe-Chiba band.

Figure 6 shows that A Keogh offers great improvement over IDDTW, as comparable retrieval performances can be obtained at a cost divided by two. This can be explained by the fact that each IDDTW comparison between the query and a candidate sequence uses several estimations to infer how likely that candidate is the best match. Running such estimations multiple times (IDDTW) introduces more errors than what a single test would do (A Keogh). That phenomenon, a well-known problem in statistics, explains why relying on A Keogh better performs than using IDDTW.

iSAX

In this section, we compare the performance of A iSAX and A iSAX Approx with the performance of the original version of iSAX2.0. We have built two indexing trees, one for learning and one for testing. The learning tree is required to determine the parameters needed by A iSAX and A iSAX Approx. Then, the learned parameters are used to construct the other tree, subsequently used for searching. Each tree stores 100,000 sequences that are random walks of length 256. Each node of the tree can store up to 1,000 sequences. The trees are built using the iSAX2.0 algorithm to ensure better balance in the leaves cardinality. Two sets of 1,000 random walk sequences are used as queries: one for learning and the other one for testing.

A groundtruth is computed in terms of DTW, as for previous experiments. Together with the 1-NN retrieval rate, we also report 10-NN (resp. 50-NN) retrieval rate that is the ratio of returned nearest neighbors that are among the true 10 (resp. 50) nearest neighbors of the query point. iSAX MinDist is used as a reference in terms of cost.

We measured average query times of around 10s for iSAX MinDist while approximate methods could perform several orders of magnitude faster: Querying the index took on average 100ms for iSAX Approx, 300ms for A iSAX with p = 0.01 and 10ms for A iSAX Approx with p = 0.01. Note that, these timings are coherent with the cost expressed here (that is related to the number of DTW computations), which is why we will keep presenting this cost in the following tables.

Results presented in Table 5 show that A iSAX outperforms iSAX MinDist and that A iSAX Approx outperforms iSAX Approx. For example, when p = 0.001, iSAX Approx runs more than 6 times more DTW computations than what A iSAX Approx does while achieving similar retrieval rates (this refers to raws #2 and #6 of the table). It is interesting to observe the retrieval 10-NN (resp. 50-NN) rates given by columns 5 and 6. The rates are excellent while the computation cost is significantly reduced.

iSAX +

We compare here the performance of iSAX2.0 and iSAX + that build their index tree using two different approaches (iSAX + uses k-means). We first ran a set of experiments using the very same datasets as the ones presented above, i.e., that are made with normally distributed random walks. It is crucial to note sequences built that are particularly adapted to iSAX2.0. As we did previously, we set the parameters for iSAX + using a first set of sequences. These parameters are then used in the experiments below. Both iSAX2.0 and iSAX + try as much as possible to determine their quantization intervals to get roughly equiprobable symbols. They use different mechanisms, however. To check the effectiveness of these mechanisms, we measured the imbalance of the final intervals determined on the one hand by iSAX2.0 and on the other after having ran the cluster balancing phase that follows the completion of the k-means in the case of iSAX + . Table 6 shows the ratio between the less and the more probable symbols from the trees build by iSAX2.0 and iSAX + when using the learning or the test tree. This table shows the iSAX + balancing phase achieves a much better balance of symbols probabilities than what iSAX2.0 does. This is a very nice result.

We now turn to the retrieval results when using iSAX2.0 or iSAX + for searching. These results are given in Table 7. iSAX2.0 and iSAX + show comparable retrieval performance when p = 0.001-slightly lower in the case of 1-NN and identical otherwise. The retrieval cost is, however, much lower. One claim made in this paper is that iSAX + better performs than iSAX2.0 when non-normally distributed sequences are used. We therefore ran experiments involving random walks using a Beta distribution with parameters α = β = 0.5. There are 100,000 such walks in the database, each has a length of 256 and 1,000 other such walks have been created for querying the system. Table 8 presents the retrieval performances and the computing costs for iSAX + and iSAX2.0 when using these sequences. A iSAX + is much faster than iSAX MinDist and iSAX Approx while achieving better recall compared to the latter.

Comparing Tables 7 and8 enlighten the benefit of relying on a method that does not make use of any assumption on the distribution of the data in space. Figure 7 give more comprehensive results for this non-normal feature setup. It shows that the best trade-offs in this case are obtained either by A iSAX Approx or A iSAX + , depending on the computational cost one can afford.

A case study: indexing DNA datasets

As shown in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], indexing full genomes using time series representations enables to link chromosomes between species that share common ancestors. For example, Rhesus macaques and humans have a common ancestor that lived 25 million years ago, which means their genomes share common attributes. Yet, the mapping between their chromosomes is not straightforward. One possible way to find such a mapping is to turn DNA sequences into time series and use related indexing techniques.

Here, we use the method proposed in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], that is, each DNA symbol is changed into a numerical value equal to the one of the previous element in the sequence plus a value that depends on the symbol. Obtained sequences are then downsampled by a factor of 25 in order to reduce noise. We finally use a sliding window of length 400 and step 10 so as to extract subsequences. We use the iSAX2.0 algorithm with parameters w = 10, b = 2, th = 1, 000 to index the whole Rhesus macaque genome (that is made of 10, 194, 500 such subsequences). The obtained index is queried using 300 sequences randomly picked from each human chromosome and their transposed versions. As in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF], we only consider longest chromosomes for the mapping, and for each of these chromosomes, we keep track of the 10 nearest neighbors among all retrieved neighbors. This case study is typical of problems for which exact methods would be really long to run, while approximate approaches can provide useful information. This in fact motivated the use for iSAX Approx rather than iSAX MinDist in [START_REF] Camerra | iSAX 2.0: Indexing and mining one billion time series[END_REF]. For this reason, we chose to use a high value for our parameter p, that is p = 0.1, so that low query times could be achieved. We present correspondence matrices in Fig. 8 by coloring each cell according to the amount of these nearest neighbors for the considered pair: the darker the cell, the more neighbors. Query times reached by A iSAX Approx are significantly lower than the ones of iSAX Approx, while obtained correspondence matrices are visually similar. Note that, in this case, performing exact search would take around 100days (estimation based on the query times measured for a few of the 12,000 queries).

Conclusion

Reducing the cost of computing the DTW similarity measure between sequences is crucial to many applications. For some applications, indexing is not appropriate. Their performance, however, can be improved by using approximate lower bounding distance functions as described in this paper. The accuracy of these functions can range from rough results with possible false dismissals returned extremely fast to exact results returned more slowly-the users can set this accuracy versus speed trade-off. For other applications, indexing is mandatory, for example to cope with scale. In this case, these approximate lower bounding distance functions can also improve the performance of retrievals, as demonstrated here. Furthermore, this paper also shows that relying on a balanced k-means quantification step significantly improves the behavior of the very popular iSAX indexing scheme. Future work includes the design of more robust indicators than the s-DTW estimate. This includes the proposition of specific lower bounds dedicated to the match of rapidly varying sequences, for which LB Keogh lacks tightness. We will also investigate the specific case of sequences of high-dimensional features, in which relying on assumptions on the distribution of the data is unaffordable. Finally, for many applications, the aim is not to match entire sequences, but rather spot similar subsequences inside long sequences. In this case, novel indexing approaches have to be developed to allow for approximate matching of subparts of time series. 8. Correspondence matrices between human and Rhesus macaque genomes. Y-axis corresponds to human chromosomes, while X-axis represents Rhesus macaque ones. Blue dots correspond to ground-truth matches. In (b), p parameter is set to 10%. Difference map is given to evaluate the similarity between both matrices, where white cells indicate perfectly similar results while black cells stand for opposite results. Cells marked with solid green line are those for which our algorithm is better while cells marked with dashed red lines are the ones for which iSAX Approx performs better. Other cells all have equal values for both techniques. Best viewed in color.

Appendices A. Proofs and mathematical definitions for upper bounds

We prove here that UB Keogh is upper bounding DTW when the latter is restricted to a Sakoe-Chiba band. We also introduce upper bounds related to LB PAA and iSAX MinDist for which we omit the proofs as they follow the exact same principles.

Definition 1. Let UB Keogh be:

UB Keogh(Q, C) = n i=1    (c i -L i ) if c i > U i (U i -c i ) if c i < L i max(U i -c i , c i -L i ) otherwise (11)
Lemma 1. For any two sequences Q and C of length n, the following inequality stands:

L 1 (Q, C) ≥ DTW(Q, C)
where the considered DTW is constrained to a Sakoe-Chiba band of width r.

Proof. Let Q and C be two sequences of length n. Manhattan distance L 1 corresponds to the alignment that follows the diagonal path. Hence, this distance is associated with one of the possible paths considered by the DTW algorithm and is therefore greater than the cost of the minimal path, that is the value returned by DTW, which concludes the proof for Lemma 1.

Proposition 1. For any two sequences Q and C of length n, the following inequality stands:

UB Keogh(Q, C) ≥ DTW(Q, C)
where the considered DTW is constrained to a Sakoe-Chiba band of width r.

Proof. It is important to notice that each term in the sum that occurs in the definition of UB Keogh is related to exactly one term in the computation of the L 1 distance. The only difference is that for UB Keogh, the i-th term corresponds to the distance between the i-th point in the candidate sequence and its furthest corresponding point in the envelope of the query, while for L 1 the same term is equal to the distance between the i-th point in the candidate sequence and one of its possible corresponding points in the envelope of the query. The latter distance is then, by definition, smaller than the former, and the following inequality is then straightforward, coming from Lemma 1:

UB Keogh(Q, C) ≥ L 1 (Q, C) ≥ DTW(Q, C). C 1 C 2 C ′ 0 C 0 ∆ 1 h C ′ 1 M ∆ 2

Fig. 1 .

 1 Fig. 1. Examples of constrained DTW paths. Admissible paths are restricted to the gray areas.

Fig. 2 .

 2 Fig. 2. Example of LB Keogh and UB Keogh. U and L are defined as in Keogh and Ratanamahatana (2005). The right hand side of each figure shows the overall area accounted for when computing each bound.

FigureFig. 3 .

 3 Fig.3. Experimental probability density function (a) and cumulative density function (b) for the whole merged dataset DS 1 and for its sub-datasets "Robot Arm" (labeled #20 later) and "Spot Exrates" (labeled #24 later).

Fig. 5 .

 5 Fig. 5. Median of tightness and its quartiles for A Keogh compared to the tightness of LB Keogh.

Fig. 6 .

 6 Fig.6. Comparison of IDDTW and A Keogh in terms of cost versus accuracy. The cost is defined as the number of elementary operations involved in the search process and accuracy is expressed in terms of 1-NN retrieval rate.

  Fig. 8. Correspondence matrices between human and Rhesus macaque genomes. Y-axis corresponds to human chromosomes, while X-axis represents Rhesus macaque ones. Blue dots correspond to ground-truth matches. In (b), p parameter is set to 10%. Difference map is given to evaluate the similarity between both matrices, where white cells indicate perfectly similar results while black cells stand for opposite results. Cells marked with solid green line are those for which our algorithm is better while cells marked with dashed red lines are the ones for which iSAX Approx performs better. Other cells all have equal values for both techniques. Best viewed in color.

Fig. 9 .

 9 Fig. 9. Balancing k-means for the k = 2 case.

Table 2 .

 2 Spearman correlation coefficient. Here, the median of ρ values is computed for all sub-datasets and the corresponding quartiles are reported.

	Dataset	Method	Median of ρ	Quartiles of ρ
		LB Keogh	0.966	[0.874;0.996]
		A Keogh p = 0.05	0.967	[0.865;0.997]
	DS 1	A Keogh p = 0.1 LB PAA	0.967 0.917	[0.861;0.997] [0.501;0.991]
		A PAA p = 0.05	0.927	[0.653;0.993]
		A PAA p = 0.1	0.926	[0.653;0.993]
		LB Keogh	0.914	[0.820;0.982]
		A Keogh p = 0.05	0.914	[0.819;0.982]
	DS 2	A Keogh p = 0.1 LB PAA	0.915 0.614	[0.819;0.982] [0.310;0.875]
		A PAA p = 0.05	0.661	[0.194;0.874]
		A PAA p = 0.1	0.663	[0.181;0.867]

Table 3 .

 3 Compared accuracy and computational cost of A Keogh and A PAA. In all cases, the cost of the exact lower bound is used as a reference.

		Method	p	Cost	1-NN	Classification
	LB Keogh	-	1.000	1.000	0.918
			0.01	0.165	0.715	0.913
	A Keogh	0.05	0.062	0.410	0.879
			0.10	0.042	0.315	0.853
	LB PAA	-	1.000	1.000	0.918
			0.001 0.770	0.942	0.918
		A PAA	0.01	0.370	0.603	0.906
			0.05	0.132	0.243	0.861
		Dataset	T	Cost	1-NN	Classi-fication
	1	50words	0.402 0.437	0.802	0.765
	2	Adiac	0.177 0.182	0.673	0.519
	3		Beef	0.581 0.526	1.000	0.500
	4		CBF	0.268 0.236	0.953	0.994
	5	Coffee	0.357 0.372	0.929	0.821
	6	ECG200	0.309 0.335	0.970	0.850
	7	FISH	0.138 0.354	0.874	0.851
	8	FaceAll	0.121 0.234	0.679	0.766
	9	FaceFour	0.114 0.728	0.966	0.841
	10	Gun Point	0.589 0.428	0.967	0.920
	11	Lighting2	0.342 0.599	0.984	0.852
	12	Lighting7	0.349 0.566	0.945	0.808
	13	OSULeaf	0.277 0.211	0.839	0.612
	14	OliveOil	0.030 0.698	0.867	0.867
	15	SwedishLeaf	0.220 0.382	0.782	0.789
	16	Trace	0.600 0.597	0.890	0.990
	17	Two Patterns	0.264 0.036	0.463	1.000
	18 synthetic control 0.086 0.551	0.860	0.987
	19	wafer	0.541 0.294	0.748	0.986
	20		yoga	0.398 0.342	0.861	0.839

Table 4 .

 4 Per-dataset accuracy and computational cost. The cost of LB Keogh is used as a reference and A Keogh is used with p set to 0.01.

Table 5 .

 5 Retrieval performance, together with the ratio of DTW computations (denoted as cost), for iSAX.

	p	Method	Cost	1-NN 10-NN 50-NN
	-	iSAX MinDist	1.0000	1.00	1.00	1.00
	-	iSAX Approx	0.0095	0.12	0.59	0.92
	0.001	A iSAX	0.0722	0.39	0.84	0.99
	0.010	A iSAX	0.0324	0.18	0.77	0.99
	0.050	A iSAX	0.0034	0.03	0.26	0.77
	0.001 A iSAX Approx 0.0014	0.12	0.57	0.92
	0.010 A iSAX Approx 0.0003	0.05	0.49	0.90
	0.050 A iSAX Approx 0.0001	0.02	0.17	0.58

Table 6 .

 6 Compared balance of iSAX and iSAX + trees.

	Tree		Imbalance factor γ Test tree Learning tree
	iSAX 2.0	1.926	1.927
	iSAX +	1.11	1.10
	p		Search type	Cost	1-NN 10-NN 50-NN
	-	iSAX MinDist	1.000	1.00	1.00	1.00
	-	iSAX + MinDist	1.7997	1.00	1.00	1.00
	-	iSAX + Approx	0.0141	0.12	0.47	0.81
	0.001		A iSAX +	0.6756	0.89	1.00	1.00
	0.010		A iSAX +	0.0400	0.46	0.99	1.00
	0.050		A iSAX +	0.0004	0.05	0.37	0.88
	0.001 A iSAX + Approx 0.0022	0.11	0.47	0.80
	0.010 A iSAX + Approx 0.0004	0.11	0.34	0.77
	0.050 A iSAX + Approx 0.0001	0.02	0.15	0.47

Table 7 .

 7 Compared performances for approximate search using iSAX + and iSAX2.0. Sequences used here are drawn from normal distribution.

  Comparison of several indexing methods in the case of non-normal features.

		1
		0.8
	10-NN retrieval rate	0.4 0.6 0.2			iSAX_MinDist iSAX_Approx A_iSAX A_iSAX_Approx iSAX+_Approx
					A_iSAX+
				A_iSAX+_Approx
		0
		0.0001	0.001	0.01	0.1	1
				Computational cost
	Fig. 7.
	p	Search type	Cost	1-NN 10-NN 50-NN
	-	iSAX MinDist	1.000	1.00	1.00	1.00
	-	iSAX Approx	0.0129	0.20	0.81	0.95
	-	iSAX + MinDist 1.9314	1.00	1.00	1.00
	-	iSAX + Approx	0.0159	0.14	0.53	0.76
	0.001	A iSAX		0.2363	0.68	0.96	1.00
	0.001	A iSAX +		0.2914	0.95	1.00	1.00
	0.010	A iSAX +		0.0125	0.24	1.00	1.00
	Table 8. Compared performances for approximate search using iSAX + and iSAX2.0: the case
	of non-normal features

The A prefix stand for approximate.

Note that, in this paper, we chose to rely on the DTW formulation that is used in the work of[START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF], which leads to the Manhattan distance being an upper bound while, when using the same formulation as in[START_REF] Keogh | Exact indexing of dynamic time warping[END_REF], DTW is upper-bounded by Euclidean distance.

There exists exactly one such node. isax approximate search is in fact the method defined in[START_REF] Shieh | iSAX: Indexing and mining terabyte sized time series[END_REF] and used in the next section describing iSAX Approx indexing.

Definition 2. Let us define UB PAA as:

Lemma 2. For any two sequences Q and C of length n, the following inequality stands:

Proposition 2. For any two sequences Q and C of length n, the following inequality stands:

where the considered DTW is constrained to a Sakoe-Chiba band of width r.

Proof. It is straightforward that proving Lemma 2 is sufficient to prove, using Proposition 1, that Proposition 2 holds. Let Q and C be two sequences of length n and W be the minimum cost path used for DTW computation with a Sakoe-Chiba band of with r between Q and C. So as to prove:

it is sufficient to prove that, for all i ∈ {1, . . . , N } :

we get :

which concludes the proof.

Definition 3. Let us define iSAX MaxDist as:

where

B. Balancing k-means

When k = 2, balancing k-means does not require any iterative process as proposed in [START_REF] Tavenard | Balancing clusters to reduce response time variability in large scale image search[END_REF]. It is possible to derive elevation h that the most populated clusters' centroid will get in order for both clusters to finally get equal populations without resorting to any iterative process. Let us assume, without loss of generality, that k-means produced two centroids C 1 and C 2 and that cluster C 1 is more populated than C 2 . Using notations introduced in Fig. 9, intersection between the line (C 1 , C 2 ) and the boundary between classes C 1 and C 2 is then C 0 , middle of the line segment [C 1 , C 2 ]. We aim at evaluating elevation h such that this point moves to C 0 that is the median of projected data points. It is straightforward that if one builds a new boundary that is parallel to the original one and passes through C 0 , both clusters will be equally populated. After solving the related system of equations, one gets:

where x 1 and x 2 are known from the k-means and x 0 is the median of projected data points.

Author Biographies

Romain Correspondence and offprint requests to: Romain Tavenard. Université de Rennes 2, Rennes, France. Email: romain.tavenard@univ-rennes2.fr


Richard Mark Soley
email: soley@omg.org

Modeling Enterprise Interoperability: Taming the Information Explosion

Keywords: Enterprise Interoperability, information explosion, Business Process Modeling, Model Driven Architecture

The problems of enterprise interoperability, portability, maintenance and integration are not exactly new. From the first time code was stored in memory, the problems of legacy integration with new users and new uses of computing systems began. The explosion of computing vendors and tools hasn't exactly made the problem any easier. In fact, it's the explosion of information in general that is causing the problem. We expect information at our fingertips, but somehow we expect that to come about magically, despite different developers, different development styles, different coding languages, operating systems, instruction set architectures --and a general lack of planning (or indeed, reading of the literature). The most important problem is the enormous explosion of information available in the world, and the increasing demands for globalized, mobile, agile, connected business processes across newly digital value chains. The resulting complexity makes integration even harder than it was before -and it was near impossible before. There is some hope, however. That hope is formal modeling, with associated metrics and continuous improvement of processes based on customer and supplier feedback. That is much easier to achieve, however, when those business models are "live" --that is, rather than simply documenting the business process, they in fact are the business process. This requires not only formal models, but formal models with well-defined semantics. The combination of Business Process Modeling (BPM) with Model Driven Architecture (MDA) promises just that. This keynote will discuss the driving factors for BPM and MDA, and the standards that support the approach.

Brief Biography

Richard Mark Soley, Ph.D is Chairman and CEO of the Object Management Group, Inc. (OMG). Dr. Soley was instrumental in the founding of OMG in 1989 (as the founding Chief Technology Officer), and since 1997 has led the organization as Chief Executive Officer.

As CEO, Dr. Soley is responsible for the vision and direction of the world's largest consortium of its type. After joining the nascent organization in 1989, Dr. Soley led the development of OMG's world-leading standardization process and the original CORBA® specification. In 1996, he led the effort to move into vertical market standards (starting with healthcare, finance, telecommunications and manufacturing) and modeling, leading first to the Unified Modeling Language TM (UML®) and later the Model Driven Architecture® (MDA®). He also led the effort to establish the SOA Consortium in January 2007, leading to the launch of the Business Ecology Initiative (BEI) in 2009. The Initiative focused on the management imperative to make business more responsive, effective, sustainable and secure in a complex, networked world, through practice areas including Business Design, Business Process Excellence, Intelligent Business, Sustainable Business and Secure Business. In addition, Dr. Soley is the Executive Director of the Cloud Standards Customer Council, helping end-users transition to cloud computing and direct requirements and priorities for cloud standards throughout the industry. He was also directly involved in the creation of both the Eclipse Foundation and Open Health Tools.

Since 1989, OMG has become one of the world leaders in the creation of international software industry standardization, building worldwide communities to deliver software standards that have changed the face of the software development field. Hundreds of member companies (approximately half software vendors and half software users, with government agencies & research institutions as well), volunteering thousands of product management and software development professionals, are currently developing some hundred software standards in fields as diverse as business processing modeling languages & methodologies, systems modeling languages, cloud computing, software modernization & real-time & embedded systems, but also vertically-oriented standards in financial services, insurance, healthcare, manufacturing, life sciences, military command & control, military & civil communications, civil government management and new areas like smart energy grids and systems safety assurance for consumer devices from automobiles to consumer electronics.

In more than two decades at OMG, Dr. Soley has led the development of a leadership community on the OMG Board of Directors that has included high-level executives from IT vendors including IBM, Hewlett Packard, Oracle and Microsoft but also key vertical-market leaders like THALES, Citigroup, HSBC, Lockheed Martin, Northrop Grumman, John Deere & Company and many others, including many small and medium-sized firms as well. Through relationships around the world, Dr. Soley has become quite well known in the speaking circuit, giving hundreds of speeches in support of OMG's best-practices & technology programs in North & South America, all over Asia, Europe, Africa and the Middle East, on issues relevant to standards, the adoption of new technology and creating successful companies. Dr. Soley also serves on numerous industrial, technical and academic conference program committees all over the world. Before helping to found OMG, Dr. Soley was a successful entrepreneur, participating in and leading startups in the computer software and hardware arena. Start-up companies included A.I. Architects (the first leader in acceleration hardware for personal computers, and originator of the DOS extender), Symbolics (the most successful start-up in artificial intelligence hardware & software), PictureTel (the early leader in video telephony, now part of Polycom) and many others. Dr. Soley continues this early entrepreneurial focus with an active participation in venture investment & advising, with companies including United Villages (the leader in rural delivery of fast-moving consumer goods in India, with a recent partnership with Oxigen), rollApp (a Ukrainian/US leader in rapid application virtualization for cloud delivery), Polymita (a Spanish company which was focused on being the world leader in business process automation, recently sold to Red Hat), Strategic Security Air (delivering lightweight, long-endurance, low-cost surveillance aircraft for civil and light military usage), ProcessUnity (a world leader in business process modeling for regulatory compliance and other innovative risk solutions), SpaceCurve (the leader in real-time geospatial-temporal databases for location services), ClinicalBox (the leader in automated surgical coordination), and so forth.

A

  native of Baltimore, Maryland, U.S.A., Dr. Soley's academic background includes the SB, SM and PhD degrees in Computer Science & Engineering from MIT, with undergraduate minor in Spanish and graduate concentration in business (with coursework from MIT's Sloan School and Harvard Law School). Dr. Soley has continued his academic & research relationships worldwide with associations as an advisor, Visiting Professor or Adjunct Professor at the Curtin University in Perth, Australia, Colorado State University, Wuhan University in the People's Republic of China, Bentley College near Boston, the Software Engineering Methods and Tools (SEMAT) international project, Fraunhofer Institut für Software-und Systemtechnik in Berlin, Goethe Universität Frankfurt am Mein and others. Dr. Soley is a Life Member of the national engineering honor fraternity Tau Beta Pi; a Member of the national electical engineering honor fraternity Eta Kappa Nu, and a Life Member of the scientific research honor fraternity Sigma Xi, as well as a member of AAAI, ACM and IEEE; he is also a Founding Member of the Boston CTO Club, and a Member of the New York City CTO Club.