Chao Zhang
email: chao.zhang@uca.fr

Reza Akbarinia
email: reza.akbarinia@inria.fr

Farouk Toumani
email: farouk.toumani@uca.fr

Zhang Chao

Limos

Efficient Incremental Computation of Aggregations over Sliding Windows

Keywords: Information systems → Data streaming Data Stream, Streaming Algorithm, Sliding Window Aggregation ACM Reference Format:

Computing aggregation over sliding windows, i.e., finite subsets of an unbounded stream, is a core operation in streaming analytics. We propose PBA (Parallel Boundary Aggregator), a novel parallel algorithm that groups continuous slices of streaming values into chunks and exploits two buffers, cumulative slice aggregations and left cumulative slice aggregations, to compute sliding window aggregations efficiently. PBA runs in 𝑂 (1) time, performing at most 3 merging operations per slide while consuming 𝑂 (𝑛) space for windows with 𝑛 partial aggregations. Our empirical experiments demonstrate that PBA can improve throughput up to 4× while reducing latency, compared to state-of-the-art algorithms.

INTRODUCTION

Nowadays, we are witnessing the production of large volumes of continuous or real-time data in many application domains like traffic monitoring, medical monitoring, social networks, weather forecasting, network monitoring, etc. For example, every day around one trillion messages are processed through Uber data analytics infrastructure 1 while more than 500 million tweets are posted on Twitter [START_REF] Gupta | Real-Time Tweet Analytics Using Hybrid Hashtags on Twitter Big Data Streams[END_REF]. Efficient streaming algorithms are needed for analyzing data streams in such applications. In particular, aggregations [START_REF] Mesiar | Aggregation Functions[END_REF], having the inherent property of summarizing information from data, constitute a fundamental operator to compute real-time statistics in this context. In the streaming setting, aggregations are typically computed over finite subsets of a stream, called windows. In particular, Sliding-Window Aggregation (SWAG) [START_REF] Li | Pane, No Gain: Efficient Evaluation of Sliding-Window Aggregates over Data Streams[END_REF][START_REF] Li | Semantics and Evaluation Techniques for Window Aggregates in Data Streams[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF] continuously computes a summary of the most recent data items in a given range 𝑟 (aka window size) and using a given slide 𝑠. If the range and slide parameters are given in time units (e.g., seconds), then the sliding window is time-based, otherwise, i.e., if these parameters are given as the number of values, it is count-based. Fig. 1 presents an example of computing sum over the count-based sliding window with a range of 10 values and a slide of 2 values.

Stream processing systems (SPSs) [START_REF] Akidau | MillWheel: Fault-Tolerant Stream Processing at Internet Scale[END_REF][START_REF] Carbone | Apache Flink: Stream and Batch Processing in a Single Engine[END_REF][START_REF] Kulkarni | Twitter Heron: Stream Processing at Scale[END_REF][START_REF] Toshniwal | Storm@twitter[END_REF][START_REF] Zaharia | Discretized Streams: Fault-Tolerant Streaming Computation at Scale[END_REF] are ubiquitous for analyzing continuously unbounded data. However, one of the challenges faced by SPSs is to efficiently compute aggregations over sliding windows. This requires the ability of such systems to incrementally aggregate moving data, i.e., inserting new data items and evicting old data items when a window is sliding without recomputing the aggregation from scratch. High throughput and low latency are essential requirements as SPSs are typically designed for real-time applications [START_REF] Hirzel | Sliding-Window Aggregation Algorithms: Tutorial[END_REF].

Two orthogonal techniques have been proposed to meet such requirements: slicing (aka partial aggregation) [START_REF] Carbone | Cutty: Aggregate Sharing for User-Defined Windows[END_REF][START_REF] Krishnamurthy | On-the-Fly Sharing for Streamed Aggregation[END_REF][START_REF] Li | Pane, No Gain: Efficient Evaluation of Sliding-Window Aggregates over Data Streams[END_REF][START_REF] Traub | Efficient Window Aggregation with General Stream Slicing[END_REF], and merging (aka incremental final aggregation) [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF]. Slicing techniques focus on slicing the windows into smaller subsets, called slices, and then computing partial aggregation over slices, called slice aggregations. The benefit of the slicing technique is that slice aggregations (i.e., partial aggregations over slices) can be shared by different window instances. For example, in Fig. 1, a slice aggregation over a slice of 2 values can be shared by 5 windows. On the basis of a slicing technique, final aggregations over sliding windows are computed by merging slice aggregations, e.g., in Fig. 1, SWAGs are computed by merging 5 slice aggregations. During the merging phase, each insertion or eviction is processed over a slice aggregation rather than a value. In modern SPSs, practical sliding windows can be very large [START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF], thereby making the merging phase non-trivial. To efficiently merge slice aggregation, incremental computation is necessary because an approach that recalculate slice aggregations from scratch (hereafter called Recal) is very inefficient [START_REF] Li | Pane, No Gain: Efficient Evaluation of Sliding-Window Aggregates over Data Streams[END_REF][START_REF] Li | Semantics and Evaluation Techniques for Window Aggregates in Data Streams[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF].

The difficulty of incremental merging depends on the considered class of aggregations: invertible or non-invertible. An aggregation is invertible if the merging function for its partial aggregations has an inverse function. For example, 𝑠𝑢𝑚 is clearly invertible, because it has the arithmetical subtraction as the inverse function. Similarly, 𝑎𝑣𝑔 and 𝑠𝑡𝑑 are invertible because they use addition as merging function. Partial aggregations of invertible aggregations can be efficiently merged using the subtract-on-evict algorithm [START_REF] Hirzel | Sliding-Window Aggregation Algorithms: Tutorial[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF], i.e., maintaining a running sum over a sliding window by subtracting evicted values and adding inserted values. However, this is not the case for non-invertible aggregations, e.g., max, min, and bloom filter, where the subtract-on-evict algorithm cannot be applied at the merging phase. Incremental merging of continuous data in the context of non-invertible aggregations is challenging and requires more sophisticated algorithms (e.g., see [START_REF] Arasu | Resource Sharing in Continuous Sliding-Window Aggregates[END_REF][START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF]).

This paper focuses on merging slice aggregations for computing non-invertible aggregations over FIFO sliding windows with an arbitrary range and an arbitrary slide. We propose PBA (Parallel Boundary Aggregator), a novel algorithm that computes incremental aggregations in parallel. PBA groups continuous slices into chunks, and maintains two buffers for each chunk containing, respectively, the cumulative slice aggregations (denoted as 𝑐𝑠𝑎) and the left cumulative slice aggregations (denoted as 𝑙𝑐𝑠) of the chunk's slices. Using such a model, SWAGs can be computed in constant time bounded by 3 for both amortized and worst-case time. Interestingly, the required 𝑐𝑠𝑎 and 𝑙𝑐𝑠 for each window aggregation are completely independent from each other, so the incremental computation of 𝑐𝑠𝑎 and 𝑙𝑐𝑠 can be achieved in parallel. These salient features put PBA ahead of state-of-the-art algorithms in terms of throughput and latency.

In this paper, we make the following main contributions:

• We propose a novel algorithm, PBA, to efficiently compute SWAGs. PBA uses chunks to divide final aggregations of window instances into sub-aggregations that are elements of 𝑐𝑠𝑎 and 𝑙𝑐𝑠 buffers. These two buffers can be computed incrementally and in parallel.

• We analyze the latency caused by SWAG computations with different chunk sizes in PBA, and propose an approach to optimize the chunk size leading to the minimum latency. • We conduct extensive empirical experiments, using both synthetic and real-world datasets. Our experiments show that PBA behaves very well for average and large sliding windows (e.g., with sizes higher than 1024 values), improving throughput up to 4× against state-of-the-art algorithms while reducing latency. For small-size windows, the results show the superiority of the non-parallel version of PBA (denoted as SBA) that outperforms other algorithms in terms of throughput. • To show the benefit of our approach in modern streamprocessing frameworks, we implemented PBA on top of Apache Flink [START_REF] Flink | [END_REF][START_REF] Carbone | Apache Flink: Stream and Batch Processing in a Single Engine[END_REF], called FPBA. Our empirical evaluation shows that FPBA scales well as increasing the parallelism of Flink in both local and cluster modes.

The rest of this paper is organized as follows. In Section 2, we present the background related to SWAG. We discuss state-of-theart algorithms for merging slice aggregations, and parallelization in SPSs in Section 3. PBA is presented in Section 4. Our experimental evaluation is reported in Section 5. We conclude in Section 6. All related proofs and pseudo-codes are included in our technical report, which is available online [1].

BACKGROUND 2.1 Sliding Window Aggregations

.).

We assume that the sliding window respects the FIFO semantic, where values that are first inserted into the window will be first evicted from it. In this case, SWAG is computed by aggregating streaming values according to their arrival order.

Two stream-processing models are widely used by SPSs. One is the micro-batch processing model, where a micro batch containing a finite subset of values is processed at a time, e.g., Apache Spark [4]. The other is the continuous processing model, where streaming values are processed one by one, e.g., Apache Flink [START_REF] Flink | [END_REF]. We consider the continuous processing model in this paper. In this case, the computation of SWAGs is triggered immediately at the boundary of window instances.

Stream Slicing

To share partial aggregations, a stream is sliced to generate disjoint subsets of streaming values, called slices. SWAGs can be computed by merging partial aggregations of slices. The state-of-the-art slicing technique is Cutty [START_REF] Carbone | Cutty: Aggregate Sharing for User-Defined Windows[END_REF], in which a partial aggregate is incrementally computed for each slice. Given a sliding window 𝑊 [𝑟, 𝑠], the size of each slice determined by Cutty is 𝑠 which equals to the slide size. For example, for the sliding window 𝑊 [START_REF] Gray | Data Cube: A Relational Aggregation Operator Generalizing Group-By, Cross-Tab, and Sub-Totals[END_REF][START_REF] Flink | [END_REF], each slice has 3 streaming values, and an instance of 𝑊 [START_REF] Gray | Data Cube: A Relational Aggregation Operator Generalizing Group-By, Cross-Tab, and Sub-Totals[END_REF][START_REF] Flink | [END_REF] will cover 3 slices and 1 streaming value. We consider Cutty as the underlying slicing technique in this paper.

Aggregation Functions

Modern frameworks provide a uniform interface to capture general aggregations. In the streaming context, aggregation functions can be decomposed into three functions: lift, combine, and lower [START_REF] Hirzel | Sliding-Window Aggregation Algorithms: Tutorial[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF], called hereafter the LCL framework. We use the geometric mean, i.e., 𝑔𝑚(𝑋 ) = ( 𝑛 𝑖=1 𝑥 𝑖 ) 1/𝑛 , as an example to illustrate them: -lift maps a streaming value to a tuple of one or several values, e.g., 𝑔𝑚.𝑙𝑖 𝑓 𝑡 (𝑣) = (𝑣, 1); -combine continuously merges two outputs of lift, e.g., 𝑔𝑚.𝑐𝑜𝑚𝑏𝑖𝑛𝑒 ((𝑣 1 , 𝑛 1 ), (𝑣 2 , 𝑛 2 )) = ((𝑣 1 × 𝑣 2 ), (𝑛 1 + 𝑛 2 )); -lower maps a tuple to a final aggregation result, e.g., 𝑔𝑚.𝑙𝑜𝑤𝑒𝑟 (𝑣, 𝑛) = 𝑣 1/𝑛 .

The combine function is also known as merging function and its outputs are called partial aggregations. We use ⊕ to denote a combine function through this paper.

An aggregation is invertible if its combine function is invertible. Let 𝑥, 𝑦 be any values from the domain of a combine function ⊕. Then ⊕ is invertible if there exists an operation ⊖, such that 𝑥 ⊕ 𝑦 ⊖ 𝑦 = 𝑥, e.g., the combine function of geometric mean is invertible as its inverse function is a pair of division and subtraction.

Aggregations can also be classified in distributive, algebraic and holistic [START_REF] Gray | Data Cube: A Relational Aggregation Operator Generalizing Group-By, Cross-Tab, and Sub-Totals[END_REF]. Distributive aggregations, e.g., sum and max, have an identical constant size for both partial and final aggregation. Algebraic aggregations have a bounded size for partial aggregation, e.g., average and geometric mean. The other aggregations are holistic, which have an unbounded size for partial aggregations, e.g., median and percentile. In this paper, we consider non-invertible aggregations that are distributive or algebraic. Note that, such consideration is consistent with state-of-the-art works [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF].

STATE OF THE ART

This section reviews state-of-the-art algorithms devoted to merging slice aggregations, and discusses slicing techniques and parallelization in SPSs. A more broad vision on SWAG can be found in [START_REF] Hirzel | Sliding-Window Aggregation Algorithms: Tutorial[END_REF] while existing surveys review various facets of SPSs [START_REF] Cugola | Processing Flows of Information: From Data Stream to Complex Event Processing[END_REF][START_REF] Heinze | Cloud-Based Data Stream Processing[END_REF][START_REF] Hirzel | A Catalog of Stream Processing Optimizations[END_REF].

Incremental Merging of Slice Aggregations

To merge slice aggregations for non-invertible aggregations, efficient algorithms are needed. The performance of such algorithms can impact the overall aggregation time, because latency caused by computing SWAG has the nature of propagation towards computations over the following ones. To our best knowledge, TwoStack [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF], DABA [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF], FlatFIT [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF], SlickDeque (Non-Inv) [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF] and SlideSide (Non-Inv) [START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF] are the most efficient state-of-the-art algorithms for computing non-invertible aggregations over FIFO windows. Consider each instance of a sliding window that requires merging 𝑛 slice (partial) aggregations. In Table 1, we present the algorithmic complexities of the state-of-the-art algorithms computed originally in [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF], and also our algorithm PBA (see Section 4.4 for explanation). TwoStack pushes pairs of a partial aggregation and a streaming value into a back stack and pops such pairs from a front stack. Computing an aggregation over a sliding window instance only needs to peek and merge partial aggregations from the front and back stack. However, when the front stack is empty, TwoStack flips all items in the back stack onto the front stack requiring 𝑂 (𝑛) time. DABA [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF] leverages TwoStack to gradually perform the flipping operation of TwoStack which amortizes its time to 𝑂 (1). As a consequence, DABA is in 𝑂 (1) worst-case time.

FlatFIT [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF] only recalculates window aggregation one time per 𝑛 + 1 slides. Intermediate aggregation results are maintained in an appropriate index structure and reused to avoid redundant computations. Hence, on average, each of 𝑛 + 1 slides can be computed by 3 merging steps. SlickDeque (Non-Inv) uses a deque to store streaming values, and not every streaming value will be added into the deque. For example, when computing max over sliding windows, SlickDeque continuously replaces the end of a deque 𝑣 ′ using a new value 𝑣 if 𝑣 > 𝑣 ′ . Adding a value 𝑣 will depend on the rank of 𝑣 w.r.t the number of values in the deque which could be 𝑛. The worst time complexity of SlickDeque is 𝑂 (𝑛) while its amortized time is less than 2 because a new value will never be compared more than twice. SlideSide (Non-Inv) uses two stacks to

Algorithm Time (amortized) Time (worst) Space

TwoStack 3 𝑛 2𝑛 DABA 5 8 2𝑛 FlatFIT 3 𝑛 2𝑛 SlickDeque < 2 𝑛 2𝑛 SlideSide 3 𝑛 2𝑛 PBA < 3 3 3𝑛+13 2
store partial aggregations, like TwoStack, but shares values between them. While, an operation similar to the flipping in TwoStack is still required, resulting in 𝑂 (𝑛) worst-case time.

As shown in Table 1, which contains complexity analysis of PBA described in Section 4.4, PBA is the best for time and space complexity, except the comparison with SlickDeque in amortized time. However, unlike SlickDeque, PBA runs in constant time, which is also experimentally demonstrated in Section 5.1.1. Such a feature of PBA guarantees its performance in the worst case. Note that amortized and worst-case time determine throughput and latency, respectively.

In addition to its low time and space complexities, one of the main advantages of PBA compared to the state-of-the-art algorithms, e.g., TwoStack [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF], DABA [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF], FlatFIT [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF], SlickDeque (Non-Inv) [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF], and SlideSide (Non-Inv) [START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF], is that the latter algorithms cannot use additional threads to get performance gains, i.e., computing intermediate results in parallel. For instance, in TwoStack, the two stacks are interdependent. In SlickDeque, the deque needs to be traversed in the worst case. However, the 𝑙𝑐𝑠 and 𝑐𝑠𝑎 buffers used for computing SWAGs in PBA are completely independent of each other, allowing their incremental computation to be done in parallel.

Combining Slicing and Merging

Various slicing techniques have been proposed in the literature, among which we could mention Panes [START_REF] Li | Pane, No Gain: Efficient Evaluation of Sliding-Window Aggregates over Data Streams[END_REF], Pairs [START_REF] Krishnamurthy | On-the-Fly Sharing for Streamed Aggregation[END_REF], Cutty [START_REF] Carbone | Cutty: Aggregate Sharing for User-Defined Windows[END_REF], and Scotty [START_REF] Traub | Efficient Window Aggregation with General Stream Slicing[END_REF]. The main goal is to determine sizes of intervals over the stream called slices which are then used to split each window into several slices in order to allow sharing the aggregations computed over slices. Slicing techniques can improve throughput by amortizing the total time of merging evenly over each stream value in each slice. However, having only a slicing technique cannot reduce latency caused by merging slice aggregations, which is especially important for a continuous stream-processing model.

A merging technique on top of a slicing one can deal with the latency issues, and additionally improves the throughput. First of all, latency can be reduced, as demonstrated in Cutty and Scotty through combining their slicing techniques with the FlatFAT [START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF] merging algorithm. Eager slicing of Scotty (using FlatFAT for merging) shows up to two orders of magnitudes reduction of latency, compared to lazy slicing (with Recal as the merging algorithm). In addition, an advanced merging algorithm can also improve throughput. Specifically, [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF] shows that TwoStack outperforms Recal, except in the case of very small window sizes, and TwoStack is always faster, up to 6×, than FlatFAT. In this paper, we adopt the Cutty slicing technique and design a constant-time merging algorithm to reduce latency and improve throughput. It is worth noting that our merging approach remains independent from slicing techniques.

Parallelization in SPSs

Parallelization in SPSs is an active research area. A large piece of work in this field focuses on inter-operator parallelism with a clear dominance of key-based splitting techniques [START_REF] Röger | A Comprehensive Survey on Parallelization and Elasticity in Stream Processing[END_REF] in shared nothing cluster architecture. Regarding intra-operator parallelism, [START_REF] Mencagli | Parallel Patterns for Window-Based Stateful Operators on Data Streams: An Algorithmic Skeleton Approach[END_REF] presents parallel patterns to handle generic stateful operators (e.g., sorting, join, grouping, . . . ). PBA fits into the category of exploiting intra-operator parallelism to incrementally compute SWAGs.

Novel SPSs are designed to exploit multi-core processors to process streams in parallel [START_REF] Koliousis | SABER: Window-Based Hybrid Stream Processing for Heterogeneous Architectures[END_REF][START_REF] Zhang | Revisiting the Design of Data Stream Processing Systems on Multi-Core Processors[END_REF]. In this line of work, SABER [START_REF] Koliousis | SABER: Window-Based Hybrid Stream Processing for Heterogeneous Architectures[END_REF] processes fixed-size batches in parallel using a hybrid architecture of CPUs and GPGPUs (i.e., modern GPUs), and merges local results of batches to produce the final result. The parallelism of SABER is suitable for the micro-batch model, giving the opportunity to process each micro batch by separate tasks (one task for each batch [START_REF] Koliousis | SABER: Window-Based Hybrid Stream Processing for Heterogeneous Architectures[END_REF]). PBA is designed for the continuous processing model. In addition, PBA merges at most 3 sub-aggregations to produce the final aggregation for arbitrary cases, i.e., requiring only constant time for merging, which is achieved by using the optimal chunk size computed based on the range and slide parameters in a continuous aggregate query.

Modern distributed SPSs [2-4, 19, 30, 32] compute SWAGs over a stream of key-value pairs in parallel. For example, in Apache Flink [START_REF] Flink | [END_REF][START_REF] Carbone | Apache Flink: Stream and Batch Processing in a Single Engine[END_REF], one can transform a DataStream into a KeyedStream with disjoint partitions, where a key will be assigned to a specific partition. SWAGs over different partitions can be computed in parallel using multiple cores of nodes in a cluster. The parallel computation in PBA is orthogonal to such distributed SPSs, since PBA computes SWAGs over values with the same key in parallel, which makes PBA suitable to extend partition-oriented parallelism in distributed SPSs.

PARALLEL BOUNDARY AGGREGATOR

In PBA, a stream is considered as a sequence of chunks having an identical number of slices. Non-invertible aggregations over window instances are computed by merging two buffers of cumulative aggregations over slice aggregations in chunks. To efficiently compute SWAG, buffers are computed in parallel using two threads. As SPSs are usually running on modern CPUs [START_REF] Zeuch | Analyzing Efficient Stream Processing on Modern Hardware[END_REF], in PBA we assume that the system is equipped with at least two cores allowing to run two tasks in parallel.

Main Ideas

Let us consider a query with the non-invertible aggregation 𝑚𝑎𝑥 over a range of 10 values and a slide of 2 values, i.e., 𝑊 [START_REF] Gray | Data Cube: A Relational Aggregation Operator Generalizing Group-By, Cross-Tab, and Sub-Totals[END_REF][START_REF] Akidau | MillWheel: Fault-Tolerant Stream Processing at Internet Scale[END_REF]. Slices of 2 values can be created, and thus each window instance needs to merge 5 slice aggregations. In general, the final aggregation needs to merge ⌊𝑟 /𝑠⌋ slice aggregations for each instance of 𝑊 [𝑟, 𝑠].

PBA is an 𝑂 (1) time solution for 𝑊 [𝑟, 𝑠], which is achieved by maintaining two buffers: 𝑐𝑠𝑎 (cumulative slice aggregations) and 𝑙𝑐𝑠 (left cumulative slice aggregations). Generally, a 𝑐𝑠𝑎 buffer is computed through accumulating slice aggregations from left to right inside a chunk (an array of slices), and a 𝑙𝑐𝑠 buffer from right to Example review of NBA with the aggregation max over a stream using PBA with a chunk size of 5 slices.

𝒘 𝟏 𝒘 𝟐 𝒘 𝟑 𝒘 𝟒 𝒘 𝟓 𝑚𝑎𝑥 𝑤 2 = 𝑚𝑎𝑥(𝑙𝑐𝑠 0 2 , 𝑐𝑠𝑎 1 [1]) 𝑚𝑎𝑥 𝑤 3 = 𝑚𝑎𝑥(𝑙𝑐𝑠 0 3 , 𝑐𝑠𝑎 1 [2]) 𝒍𝒄𝒔 𝟎 𝒄𝒔𝒂 𝟎 𝒄 𝟎 𝑚𝑎𝑥 𝑤 0 = 𝑐𝑎 0 𝑚𝑎𝑥 𝑤 1 = 𝑚𝑎𝑥(𝑙𝑐𝑠 0 1 , 𝑐𝑠𝑎 1 [0]) 5
left. Then SWAG can be computed by merging 2 elements respectively from 𝑙𝑐𝑠 and 𝑐𝑠𝑎 for this example. The two buffers of PBA for computing 𝑚𝑎𝑥 over 𝑊 [START_REF] Gray | Data Cube: A Relational Aggregation Operator Generalizing Group-By, Cross-Tab, and Sub-Totals[END_REF][START_REF] Akidau | MillWheel: Fault-Tolerant Stream Processing at Internet Scale[END_REF] are shown in Fig. 2, e.g., 𝑐𝑠𝑎 0 and 𝑙𝑐𝑠 0 for chunk 𝑐 0 . Let us explain in detail how the two buffers are computed. Streaming values are partitioned into chunks having an identical number of slices, e.g., 𝑐 0 and 𝑐 1 in Fig. 2 have 5 slices, and chunks are separated by boundaries, e.g., 𝑏 1 between 𝑐 0 and 𝑐 1 . Slice aggregations inside a chunk are sequentially stored in a 𝑠𝑎𝑎 (slice aggregation array), e.g., 𝑠𝑎𝑎 0 and 𝑠𝑎𝑎 1 in Fig. 2. Then, 𝑐𝑠𝑎 and 𝑙𝑐𝑠 are computed over 𝑠𝑎𝑎, which are two arrays having the same length as a 𝑠𝑎𝑎 array. For example, 𝑐𝑠𝑎 1 and 𝑙𝑐𝑠 0 are computed by memorization of each intermediate result during accumulating elements of 𝑠𝑎𝑎 1 and 𝑠𝑎𝑎 0 , illustrated as follows.

-𝑐𝑠𝑎 1 is obtained by accumulating from

𝑠𝑎𝑎 1 [0] to 𝑠𝑎𝑎 1 [4], -𝑙𝑐𝑠 0 is obtained by accumulating from 𝑠𝑎𝑎 0 [4] to 𝑠𝑎𝑎 0 [0].
The last element in a 𝑐𝑠𝑎 array is denoted as a 𝑐𝑎 (chunk aggregation), e.g., 𝑐𝑎 1 = 𝑐𝑠𝑎 1 [4] = 9. Aggregation over a window instance is computed using elements in 𝑙𝑐𝑠 and 𝑐𝑠𝑎 arrays, e.g., in Fig. 2,

𝑚𝑎𝑥 (𝑤 3 ) = 𝑚𝑎𝑥 (𝑙𝑐𝑠 0 [3], 𝑐𝑠𝑎 1 [2]
). Consequently, 𝑚𝑎𝑥 over each window instance, e.g., each of (𝑤 0 , ..., 𝑤 4 ), only requires at most 1 𝑚𝑎𝑥 operation.

Challenges in PBA: In order to obtain the constant-time solution described above, the two buffers, 𝑐𝑠𝑎 and 𝑙𝑐𝑠, must be efficiently computed as streaming values keep arriving. 𝑐𝑠𝑎 can be computed by accumulating streaming values and slice aggregations. The main difficulty is how to compute 𝑙𝑐𝑠 efficiently, since computing a 𝑙𝑐𝑠 array requires traversing a 𝑠𝑎𝑎 array inversely, which means having to stop computing 𝑐𝑠𝑎, resulting in an eventual high delay. We use a parallel strategy to deal with this issue. Our key observation is that SWAG never uses 𝑐𝑠𝑎 and 𝑙𝑐𝑠 from the same chunk, so the inputs required to compute 𝑐𝑠𝑎 and 𝑙𝑐𝑠 are completely independent. For example, in Fig. 2, computing 𝑚𝑎𝑥 (𝑤 0 ) requires 𝑙𝑐𝑠 0 and 𝑐𝑠𝑎 1 . 𝑙𝑐𝑠 0 is computed with 𝑠𝑎𝑎 0 as input, which is complete at the boundary 𝑏 1 . While 𝑐𝑠𝑎 1 is computed with 𝑠𝑎𝑎 1 as input, which is independent of 𝑠𝑎𝑎 0 . The independence between the inputs of 𝑙𝑐𝑠 0 and 𝑐𝑠𝑎 1 allows the task of computing 𝑙𝑐𝑠 0 to be carried out when streaming values reach the boundary 𝑏 1 . This task runs in parallel and simultaneously with the task of computing 𝑐𝑠𝑎 1 . Thus, PBA can keep receiving new streaming values after 𝑏 1 and compute 𝑐𝑠𝑎 1 . Another issue is whether the task of computing 𝑙𝑐𝑠 0 can be finished before the end of 𝑤 1 , i.e., whether PBA needs to wait for the result of 𝑙𝑐𝑠 0 [1] to compute 𝑚𝑎𝑥 (𝑤 1 ). Such an issue is dealt with in Section 4.3, where we identify the optimal chunk size minimize the waiting time.

PBA Model

In PBA, a stream is considered as a sequence of non-overlapping chunks of slices. Each chunk has an identical number of slices and starts at a specific time denoted as boundary.

-𝑏 𝑖 (boundary): the start of the chunks, e.g., 𝑏 1 in Fig. 2; -𝑐 𝑖 (chunk): a sequence of slices from 𝑏 𝑖 to 𝑏 𝑖+1 , e.g., 𝑐 0 in Fig. 2; -𝑏 (chunk size): the number of slices in a chunk, e.g., 5 in Fig. 2.

We incrementally compute partial aggregations for every chunk to compute SWAGs. PBA applies two kinds of incremental aggregations: (i) computing 𝑐𝑠𝑎 and (ii) computing 𝑙𝑐𝑠. To compute 𝑐𝑠𝑎, PBA adopts a two-level incremental computation inside each chunk: (i.1) accumulating streaming values inside each slice to have a slice aggregation, and (i.2) accumulating slice aggregations inside each chunk to have a chunk aggregation. To compute 𝑙𝑐𝑠, PBA maintains intermediate results during accumulating slice aggregations from right to left in a slice aggregation array 𝑠𝑎𝑎 𝑖 .

Dividing SWAG into Sub-Aggregations: As the window instances, whose starts coincide with boundaries in every chunk, are trivial cases requiring only 𝑐𝑠𝑎 elements e.g., 𝑤 0 in Fig. 2, we focus on the other cases in the sequel. To illustrate the PBA model, we first define 𝑘 = ⌊ ⌊𝑟 /𝑠 ⌋ 𝑏 ⌋, which denotes the number of chunks a window instance can fully cover, e.g., 𝑘 = ⌊ ⌊10/2⌋ 5 ⌋ = 1 in Fig. 2. In the PBA model, we use boundaries to divide a window instance 𝑤 of 𝑊 [𝑟, 𝑠] into several sub-parts and merge partial aggregations of sub-parts to have a final aggregation of 𝑤. The number of sub-parts can be 𝑘 + 1 or 𝑘 + 2, depending on 𝑊 [𝑟, 𝑠] and the used chunk size. We denote a partial aggregation of each sub-part as a sub-aggregation, and all sub-aggregations of 𝑤 as (𝑆𝐴 0 , 𝑆𝐴 1 , ...), i.e., 𝑆𝐴 0 is the first one (see Fig. 3 for an example). Then the first sub-aggregation is an element in 𝑙𝑐𝑠, and the last one in 𝑐𝑠𝑎. Each of the others is a chunk aggregation, i.e., the last element in a 𝑐𝑠𝑎 buffer. An example of 𝑘 = 1 is shown in Fig. 2, where each window instance is divided into 2 sub-parts and SWAGs are computed by merging 2 sub-aggregations. To illustrate the case of 𝑘 > 1, Fig. 3 presents the example of computing a count-based 𝑊 [16, 1] with a chunk size of 4, where we have 𝑘 = 4 and each window instance is divided into 5 sub-parts, e.g., 𝑚𝑎𝑥

(𝑤 1 ) = 𝑚𝑎𝑥 (𝑙𝑐𝑠 0 [1], 𝑐𝑎 1 , 𝑐𝑎 2 , 𝑐𝑎 3 , 𝑐𝑠𝑎 4 [0]).
Efficient Computation of 𝑙𝑐𝑠: In the final aggregation stage, since 𝑙𝑐𝑠 and 𝑐𝑠𝑎 buffers used for merging are computed over completely independent 𝑠𝑎𝑎 arrays, we compute 𝑙𝑐𝑠 and 𝑐𝑠𝑎 buffers in parallel, i.e., two tasks run simultaneously and each task focus on computing a single buffer. Specifically, we maintain two threads in PBA. The main thread is used to compute 𝑐𝑠𝑎 and the final aggregation, and the other for processing the task of computing 𝑙𝑐𝑠

Optimal Chunk Size in PBA

The chunk size is crucial for PBA. Specifically, a carelessly selected chunk size can cause delay to compute final aggregations. For instance, consider the example in Fig. 2, and let 1 𝑇 𝑢 be the time required to execute ⊕, i.e., a merging function or a combine function. Although the two tasks for computing 𝑙𝑐𝑠 0 and 𝑐𝑠𝑎 1 are simultaneously started at the boundary 𝑏 1 , 𝑙𝑐𝑠 0 [1] cannot be obtained earlier than 𝑐𝑠𝑎 1 [0] because computing 𝑙𝑐𝑠 0 [1] (3 𝑇 𝑢 ) requires two more 𝑇 𝑢 than computing 𝑐𝑠𝑎 1 [0] (1 𝑇 𝑢 ). Thus, PBA needs to wait for 2 𝑇 𝑢 to use 𝑙𝑐𝑠 0 [1] and 𝑐𝑠𝑎 1 [0] to compute 𝑚𝑎𝑥 (𝑤 1 ).

The following two theorems identify the optimal chunk size, denoted as 𝑏 * , for count-based windows and time-based window, respectively. The corresponding proofs are provided in our technical report, which is available online [1].

Theorem 1. Given a count-based sliding window 𝑊 [𝑟, 𝑠], the optimal chunk size 𝑏 * in PBA can be obtained as follows: -if ⌊𝑟 /𝑠⌋ ⩽ 𝑠 + (𝑟 mod 𝑠) + 1, then 𝑏 * = ⌊𝑟 /𝑠⌋;

-otherwise 𝑏 * = ⌊ 𝑠 (⌊𝑟 /𝑠⌋ + 1) + 𝑟 mod 𝑠 + 1 𝑠 + 1 ⌋.
With 𝑏 * , we have 𝑘 = ⌊ ⌊𝑟 /𝑠 ⌋ 𝑏 * ⌋ = 1. Theorem 2. Given a time-based sliding window 𝑊 [𝑟, 𝑠], the optimal chunk size 𝑏 * in PBA can be obtained as follows:

-if ⌊𝑟 /𝑠⌋ ⩽ 𝑠+2𝑇 𝑢 +𝐶 𝑇 𝑢 , then 𝑏 * = ⌊𝑟 /𝑠⌋; -otherwise 𝑏 * = ⌊ ⌊𝑟 /𝑠 ⌋ (𝑠+𝑇 𝑢 )+𝑠+2𝑇 𝑢 +𝐶 𝑠+2𝑇 𝑢 ⌋,
where 𝐶 = 0 if 𝑟 mod 𝑠 = 0, otherwise 𝐶 = 𝑠 (𝑟 mod 𝑠) + 1 𝑇 𝑢 . With 𝑏 * , we have 𝑘 = ⌊ ⌊𝑟 /𝑠 ⌋ 𝑏 * ⌋ = 1. The optimal chunk size ensures that (i) each 𝑙𝑐𝑠 buffer can be obtained without waiting, and (ii) 𝑘 = 1 which leads to the minimum number of merging operations to compute final aggregations.

Complexity Analysis

As previous works [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF][START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF] present complexities by computing a count-based 𝑊 [𝑛, 1] with a range of 𝑛 and a slide of 1, i.e., merging 𝑛 slice aggregations for each window instance, we also rest on such a scenario to discuss complexities of PBA.

Given a sliding window 𝑊 [𝑛, 1], as PBA is used with the optimal chunk size 𝑏 * , we have 𝑏 * = ⌊ 𝑛+2 2 ⌋ and 𝑘 = ⌊ ⌊𝑟 /𝑠 ⌋ 𝑏 * ⌋ = 1 according to Theorem 1. Thus, a window instance will be divided into either 𝑘 + 1 = 2 or 𝑘 + 2 = 3 sub-parts, i.e., spanning 2 or 3 chunks. Based on this, we discuss the time and space complexity of PBA below.

Time Complexity: In the final aggregation stage, merging 2 or 3 sub-aggregations of sub-parts needs 1 or 2 merging operations. Considering that the last streaming value of a window will be merged with 𝑐𝑠𝑎, the total number of merging operations by PBA is 2 or 3. Therefore, the worst-case time is 3, and the amortized time is less than 3.

Space Complexity: As a window instance can span at most 3 chunks, w.l.o.g, we denote the three chunks as 𝑐 𝑖-2 , 𝑐 𝑖-1 , 𝑐 𝑖 , and consider the current streaming value is arriving at a slice of 𝑐 𝑖 . Thus, 3 𝑠𝑎𝑎 arrays of size 𝑏 * are needed. Note that the 𝑙𝑐𝑠 arrays of the three chunks will be computed over 𝑠𝑎𝑎 arrays and do not require extra space. In addition, PBA only needs the 𝑐𝑠𝑎 array of the current chunk 𝑐 𝑖 , i.e., 𝑐𝑠𝑎 𝑖 , and maintains one element in 𝑐𝑠𝑎 𝑖 , the current one, rather than the entire 𝑐𝑠𝑎 𝑖 array. For windows spanning 3 chunks, the chunk aggregation 𝑐𝑎 𝑖-1 of chunk 𝑐 𝑖-1 is also needed. Thus, the total space required by PBA is:

3𝑏 * + 2 = 3⌊ 𝑛+2 2 ⌋ + 2 ≤ 3𝑛+13 2 .

EXPERIMENTAL EVALUATION

In this section, we evaluate the performance of PBA against the state-of-the-art algorithms for computing non-invertible aggregations over sliding windows with FIFO semantics. Moreover, we integrate PBA into Apache Flink and present the corresponding improvement in local and cluster modes, respectively.

PBA Compared to Alternatives

State-of-the-Art Algorithms: In our experimental evaluation we compare PBA to the two best state-of-the-art algorithms: SlickDeque [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF] and TwoStack [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF]. SlickDeque shows the best performance in terms of throughput and latency [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF] while TwoStack is slightly worse than SlickDeque but behaves better than the others [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF][START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF]. Note that, although DABA [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF] has a constant time complexity in the worst case, it does not show a better performance result compared to SlickDeque in terms of throughput and latency [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF], or TwoStack in terms of throughput [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF]. We also consider the naive algorithm algorithm and call it Recal (recalculate). As the sliding window moves, it calculates an aggregation from scratch. To study the performance of our algorithm in the case of single-core machines, we implemented a sequential version of PBA, denoted as SBA (sequential boundary aggregator), which uses only one thread to compute SWAG. SBA follows the same procedure as PBA, except that the task of computing 𝑙𝑐𝑠 in SBA is done in the same thread as the one that computes 𝑐𝑠𝑎 buffer (SBA uses only one thread). In order to slice a data stream, we set each slice size to be equal to the slide size (as described in Section 2.2 for details). Notice that this is the underlying slice solution for all algorithms tested in our evaluation. We implemented PBA and SBA using Java 8, and the state-of-the-art algorithms based on the pseudo-codes in [START_REF] Tangwongsan | Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time[END_REF] and [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF]. The source codes of PBA, SBA, and the implementation of state-of-the-art algorithms are available online [1]. Datasets: We use two datasets, a real-world and a synthetic one, for throughput experiment. The real-world dataset is the DEBS12 Grand Challenge Manufacturing Equipment Dataset [START_REF] Jerzak | The DEBS 2012 Grand Challenge[END_REF], which is also used by SlickDeque [START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF]. This dataset contains energy consumption recorded by sensors of manufacturing equipment with the frequency of 100 records per second, and each tuple is associated with a timestamp, which we will use in our time-based experiments. The original dataset contains 32 million tuples, which we made into 130 million tuples, i.e., 4 copies. The synthetic dataset contains 200 million values generated using the uniform distribution. Note that in our experiments we observed the same tendency for PBA with different datasets, this is why we have not used more than two test datasets. For the latency experiments, we used the synthetic dataset and tested latency results for 1 million window instances.

In our experiments, as aggregation function we use max, commonly used in related works [START_REF] Anatoli | FlatFIT: Accelerated Incremental Sliding-Window Aggregation For Real-Time Analytics[END_REF][START_REF] Anatoli | SlickDeque: High Throughput and Low Latency Incremental Sliding-Window Aggregation[END_REF][START_REF] Tangwongsan | General Incremental Sliding-Window Aggregation[END_REF][START_REF] Theodorakis | SlideSide: A Fast Incremental Stream Processing Algorithm for Multiple Queries[END_REF] as the representative of non-invertible aggregations. We first evaluate throughput for count-based windows with a fixed slide while varying the ranges. Then we evaluate the case with a fixed range while varying the slides. Similarly, we evaluate the throughput for time-based windows in both settings (i.e., with a fixed slide and varying ranges, and then with a fixed range and varying slides). We test latency for both count-based and time-based windows. To further study the latency of different algorithms, we also test a more complex aggregation, MomentSketch = (𝑚𝑖𝑛, 𝑚𝑎𝑥, 𝑠𝑢𝑚, 𝑐𝑜𝑢𝑛𝑡, 𝑥 𝑖 , ... 𝑥 𝑙 𝑖 , 𝑙𝑜𝑔𝑥 𝑖 , ... (𝑙𝑜𝑔𝑥 𝑖 ) 𝑙 ) [START_REF] Gan | Moment-Based Quantile Sketches for Efficient High Cardinality Aggregation Queries[END_REF] with 𝑙 = 10, which can approximate percentiles.

Experiment Setup: We run the throughput and latency experiments at a server with Ubuntu 16.01.3 LTS, 6 virtual CPUs of Intel(R) Xeon(R) 2.40GHz, and 16GB main memory, where the heap size of JVM is configured to 14GB.

Throughput Experiments.

The workload for evaluating max over count-based windows is as follows: 1) range from 2 3 to 2 20 with the slide size set to one; 2) slide from 1 to 10 with a range of 2 17 . For the time-based windows the workload is as follows: 1) range from 6 to 60 minutes with a slide of 10 milliseconds; 2) slide from 10 milliseconds to 100 milliseconds with a range of 1 hour. We do not test for larger slide sizes as the overall throughput will be amortized by the time of computing slice aggregations, such that the differences between merging algorithms cannot be observed (see Section 3.2 for detailed discussion). We run the workload 20 times and report the median of all evaluation results. We stop reporting the throughput of the Recal algorithm when the corresponding throughput is too small. Our experimental results are presented in Fig. 4, Fig. 5, and Fig 6 . We discuss the results in the sequel.

Count-based windows with a fixed slide. Our throughput results for this case are shown in Fig. 4. PBA shows the best throughput with large window sizes, and SBA (the sequential version of PBA) shows a higher throughput than TwoStack and SlickDeque, even for small window sizes. More precisely, the throughput of PBA increases w.r.t to the increase of window size 𝑟 when 𝑟 < 2 15 , and keeps consistent when 𝑟 ⩾ 2 15 . When 𝑟 ⩾ 2 15 , the throughput of PBA is 2.5× higher than that of TwoStack and 4× higher than that SlickDeque. Note that, the throughput of PBA increases w. window size 𝑟 when 𝑟 < 2 15 , because a smaller window size will lead to a smaller optimal chunk size in PBA. Thus, the parallel task in PBA, i.e., computing 𝑙𝑐𝑠, will consequently require less time but will be launched more frequently, compared to a case with a larger window size. In such a case, the constant overhead for launching a parallel task will be relatively more obvious. For example, when 𝑟 = 2 5 , the optimal chunk size of PBA is 17. Then, the corresponding parallel computation for each chunk in PBA requires applying only 15 times of the 𝑚𝑎𝑥 operation over 2 elements, and this cost is much less than the cost of launching a parallel task. But, if 𝑟 = 2 20 , the optimal chunk size is 5242884, which requires 5242882 times of 𝑚𝑎𝑥 over 2 elements. In the latter case, launching a parallel task can bring significant benefits.

Count-based windows with various slide sizes. Our throughput results for this case are shown in Fig. 5. We observe that, with a larger slide size, combining slicing with merging can improve throughput. PBA shows better throughput with a small slide shown in Fig. 5, because 𝑟𝑎𝑛𝑔𝑒 𝑠𝑙𝑖𝑑𝑒 is still large in such cases, which means there are still a large number of slice aggregation that need to be merged. Then, as the slide size is getting larger, the throughput of different algorithms are getting closer. There are mainly two reasons that explain such a behavior. The first reason comes from the fact that the ratio 𝑟𝑎𝑛𝑔𝑒 𝑠𝑙𝑖𝑑𝑒 is smaller in these cases, which leads to less slice aggregations to be merged for each window instance. Thus, the required computations are not enough to exploit the full capacity of PBA. The second reason is that the overall difference between various algorithms will be amortized by the cost of computing partial aggregations over each slice (this is due to the large number of ⊕ operations performed in each slice). Time-based windows with various range and slide sizes. Our throughput results of this case are shown in Fig. 6. PBA shows the overall best throughput with large window sizes, and SBA (the sequential version of PBA) shows a higher throughput than TwoStack and SlickDeque. The corresponding reasons are identical to the case of count-based windows. Note that, in Fig. 6 (a), PBA outperforms the others because 𝑟𝑎𝑛𝑔𝑒 𝑠𝑙𝑖𝑑𝑒 is large even for the case of 𝑟𝑎𝑛𝑔𝑒 = 6 minutes.

Latency

Experiments. The workload consists in: (i) computing the aggregation Max over count-based windows with a slide of 1 and ranges of 2 13 and 2 14 respectively, and (ii) computing the aggregation MS (MomentSketch) over time-based windows with a slide of 1 second (containing around 100 streaming values) and ranges of 1 and 2 hours respectively. In this experiment, we compare TwoStack, SlickDeque, SBA and PBA. We capture the computation time in nanoseconds, corresponding to the time required by an algorithm to compute the aggregation over each window when the ends of windows are reached. We run the experiments for 1 million windows, and report the latency for the last 970 thousand ones.

Our latency results are presented in Fig. 7. In Table 2 and Table 3, we present 7 statistical summaries for the latency experiments, i.e., Min (minimum), Max (maximum), Avg (average), Std (standard deviation), LQ (lower quartile, i.e., 25% percentile), Med (median), HQ (higher quartile, i.e., 75% percentile). On the whole, we observe that both SlickDeque and PBA do not have latency spikes, and PBA performs better than SlickDeque. We discuss detailed results of all tested algorithms in the sequel.

As it can be observed, PBA has the overall best latency results in terms of all statistical summaries, especially for the Max and Std metrics shown in Table 2 and3 SlickDeque has almost the same latency as TwoStack for Avg, but much better than TwoStack for Max and Std. As shown in Fig. 7, TwoStack has latency spikes periodically, and thus TwoStack does not perform well in terms of Max and Std. Such latency spikes become larger when the window size increases. TwoStack requires moving all elements of one stack to another one periodically, and the corresponding cost is linearly increasing w.r.t to the window size. This is why the periodical red spike of TwoStack for the case of 16384 is higher than the corresponding one for the case of 8192, and similar for time-based window with a range of 2 hours compared to a range of 1 hour.

We also observe that, as aggregation getting more complex, the latency is also getting higher, which can be observed in Table 2 and Table 3. As an aggregation defined in the LCL framework can have an arbitrarily complex ⊕ operation, the difference in the latency caused by using a different merging algorithm can be more obvious.

Apache Flink Integration

We integrate PBA into Apache Flink [START_REF] Flink | [END_REF][START_REF] Carbone | Apache Flink: Stream and Batch Processing in a Single Engine[END_REF], denoted as FPBA. PBA is integrated into Flink through KeyedProcessFunction, i.e., a lowlevel operator in Flink. To support parallel tasks of PBA, a thread pool for each task slot in Flink is launched, which will be used to compute the 𝑙𝑐𝑠 buffers for all key-value pairs associated with the corresponding task slot. Our goal in this subsection is to show that can be easily integrated into modern distributed and parallel SPSs, and to evaluate FPBA in both local and cluster modes. In this experiment, we consider 2 10 as the representative of large windows. We compute 𝑚𝑎𝑥 over count-based sliding windows with a slide of 1 and a window size of 2 10 , and capture the throughput of Flink and FPBA by using a different degree of parallelism in local mode, and a different number of nodes in cluster mode.

Datasets: We use the DEBS14 Grand Challenge Smart Home Dataset [START_REF] Jerzak | The DEBS 2014 Grand Challenge[END_REF], which contains recordings from smart plugs in private households. Each tuple of this dataset includes a power consumption and a plugin ID. The first 10 million tuples contain 14 plugin IDs with a slight data skew, which we made into different sizes with the identical distribution in our experiments, i.e., 60 million tuples with 84 plugin IDs for local execution, and 240 million tuples with 336 plugin IDs for cluster execution.

Experiment Setup: We run the cluster experiment at a Flink cluster having 1 master node and 6 worker nodes running Ubuntu server 16.01.3 LTS. The master node has 6 virtual CPU cores of Intel XEON 3GHz, and 16GB main memory, and the job manager of Flink is configured with 6GB of memory. Every worker node has 4 virtual CPU cores of Intel XEON 3GHz, and 8GB main memory, and the task manager of Flink is configured with 6GB of memory. We run the local experiment only at the master node of the Flink cluster. We run both local and cluster experiments 10 times, where we consider the first 3 runs as warm-ups, and we repeat the whole execution 3 times. We report the median throughput for Flink and FPBA. We use Flink 1.10.0 and Java 8 in this experiment.

We present the results in Fig. 8. In Fig. 8 (a), we observe the throughput of FPBA or Flink linearly increases up to 3 parallel instances (called parallelism in Flink). Each parallel instance of the execution pipeline of FPBA, or Flink, contains 2 tasks, i.e., reading data and computing SWAG. Therefore, up to 3 parallel instances, the task for computing SWAG can run on a dedicated CPU core at a machine with 6 cores. For more parallel instances, the throughput growth becomes less linear compared to the first phase. In Fig. 8 (b), we observe the throughput of FPBA and Flink linearly scale up w.r.t more nodes used in a cluster.

CONCLUSION

In this paper, we presented PBA, a novel algorithm that calculates incremental aggregations in parallel for efficiently computing SWAGs, and also its sequential version SBA. We also proposed an approach to optimize the chunk size, which guarantees the minimum latency for PBA. Our experimental results demonstrate that SBA has higher throughput than state-of-the-art algorithms irrespective of window sizes, and PBA shows best performance for large windows (i.e., improving throughput up to 4× while reducing latency). Thus, a simple strategy for obtaining high performance can be using SBA for small windows and PBA for larger ones. In addition, FPBA, an integration of PBA into Apache Flink, shows significant improvement in both local mode and cluster mode.

As future works, we plan to extend PBA to support both multiquery processing and out-of-order streams. The idea is to replace the array structure of 𝑙𝑐𝑠 buffers with a tree structure, where each leaf node stores a slice aggregation and any intermediate node stores partial aggregation of all its children (e.g., in the spirit of the array-based binary heap implementation of FlatFAT). In such a scenario, the major benefit of using PBA model is that any updates caused by late values can be processed only over 𝑙𝑐𝑠 buffers in parallel, such that the computation for accumulating streaming values will not be delayed. One technical challenge underlying the development of such an approach lies in optimizing the chunk size that ensures the overall minimum latency, while taking into account watermarks and multiple combinations of ranges and slides.

Figure 1 :

 1 Figure 1: Example of computing sum over a count-based sliding window with a range of 10 values and a slide of 2 values.

A

  sliding window 𝑊 [𝑟, 𝑠], 𝑟 > 𝑠, with a range 𝑟 and a slide 𝑠 is an infinite sequence of window instances (𝑤 0 , 𝑤 1 , 𝑤 2 , ...). Sliding window aggregation means computing an aggregate function 𝛼 over each window instance of a sliding window𝑊 [𝑟, 𝑠], i.e., 𝛼 (𝑊 [𝑟, 𝑠]) = (𝛼 (𝑤 0 ), 𝛼 (𝑤 1 ), 𝛼 (𝑤 2 ), ..

Figure 2 :

 2 Figure 2: Example of computing the count-based 𝑊 [10, 2]with the aggregation max over a stream using PBA with a chunk size of 5 slices.

Figure 4 :Figure 5 :

 45 Figure 4: Throughput for count-based windows by varying the window size and with a fixed slide size.

Figure 6 :

 6 Figure 6: Throughput for time-based windows.

Figure 7 :

 7 Figure 7: Latency experiments.

Figure 8 :

 8 Figure 8: Throughput on top of Apache Flink.

Table 1 :

 1 Algorithmic Complexities (Non-Invertible SWAG).

  r.t to the

		Recal	TwoStack	SlickDeque	SBA	PBA
	Million items/s	Window size 2 3 2 4 2 5 2 6 2 7 2 8 2 9 2 10 2 11 2 12 2 13 2 14 2 15 2 16 2 17 2 18 2 19 2 20 0 50 100 150 200 250 300 50 350 10 20 30 40 Million events/s Window size 0 2 3 2 4 2 5 2 6 2 7 2 8 2 9 2 10 2 11 2 12 2 13 2 14 2 15 2 16 2 17 2 18 2 19 2 20
		(a) DEBS'12 dataset [15]

Table 2 :

 2 latency summaries of Max over count-based windows.

	Algo. Range Min	Max	Avg	Std	LQ, Med, HQ
	TS	2 13	57	74543 69.21	754.3	60, 60, 61
	SD	2 13	57	667	81.77	22.12	64, 79, 88
	SBA	2 13	43	15160 50.39 233.64	45, 46, 47
	PBA	2 13	40	559	53.93	10.41	45, 48, 65
	TS	2 14	63	150024 77.55 1077.86	67, 68, 70
	SD	2 14	61	647	83.07	21.86	65, 80, 90
	SBA	2 14	43	31031 50.06 331.25	45, 46, 47
	PBA	2 14	42	537	46.06	5.24	44, 45, 47

Table 3 :

 3 latency summaries of MS over time-based windows. PBA shown in Table1. Moreover, PBA shows better latency results with a larger window size, shown in Fig.7, and Table2and 3 (the reason has been discussed in Section 5.1.1). SBA shows similar results as PBA except for Max and Std in Table2 and 3. SBA computes a 𝑙𝑐𝑠 buffer at each boundary, which postpones the computation of a 𝑐𝑠𝑎 buffer until the computation of 𝑙𝑐𝑠 has been completed. Therefore, SBA has a high latency periodically, e.g., in Fig.2, at boundary 𝑏 1 , SBA needs to firstly compute 𝑙𝑐𝑠 0 , and then 𝑐𝑠𝑎 1 , and this incurs high latency for 𝑤 1 . Such a high latency of SBA can be observed in Fig.7. Note that, PBA does not have high latency due to computing 𝑙𝑐𝑠 and 𝑐𝑠𝑎 in parallel.

	Algo. Range Min	Max	Avg	Std	LQ, Med, HQ
	TS	1 h	203 365582 337.65 5348.91 223, 230, 258
	SD	1 h	267	13638 385.28	101.5	336, 370, 415
	SBA	1 h	131 103716 194.85 1566.49 160, 163, 166
	PBA	1 h	129 12837 164.24 73.35	141, 145, 168
	TS	2 h	207 934212 381.1 9999.85 226, 235, 290
	SD	2 h	266	13175 392.72	84.86	342, 380, 427
	SBA	2 h	132 226099 201.77 2471.56 162, 164, 168
	PBA	2 h	135 12971 171.75 64.59	145, 150, 191

ACKNOWLEDGMENTS

This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25).


Heraldo Borges

Reza Akbarinia
email: reza.akbarinia@inria.fr

Florent Masseglia
email: florent.masseglia@inria.fr

Anomaly Detection in Time Series

Keywords:

Data mining has become an important task for researchers in the past few years, including detecting anomalies that may represent events of interest. The problem of anomaly detection refers to discovering the patterns that do not conform to expected behavior. This paper analyzes recent studies on the detection of anomalies in time series. The goal is to provide an introduction to anomaly detection and a survey of recent research and challenges. The article is divided into three main parts. First, the main concepts are presented. Then, the anomaly detection task is defined. Afterward, the main approaches and strategies to solve the problem are presented.

Introduction

In many real-world applications including load demand forecasting [START_REF] Qiu | Empirical mode decomposition based ensemble deep learning for load demand time series forecasting[END_REF], geography [START_REF] Feyrer | Trade and income-exploiting time series in geography[END_REF], human activity recognition [START_REF] Ignatov | Real-time human activity recognition from accelerometer data using convolutional neural networks[END_REF], stock return [START_REF] Koijen | The cross-section and time series of stock and bond returns[END_REF] and others [START_REF] Chandola | Anomaly detection: A survey[END_REF], the data is collected in the form of time series. Anomaly detection in this type of data refers to discovering any abnormal behavior within the data encountered in a specific time interval. Anomaly detection has been widely used in several application areas. For instance, cardiologists are interested in identifying anomalous parts of ECG signals to diagnose heart disorders [START_REF] Pourbabaee | Deep convolutional neural networks and learning ECG features for screening paroxysmal atrial fibrillation patients[END_REF]. Economists are interested in anomalous parts of share prices to analyze and build economic models [START_REF] Lahmiri | A variational mode decompoisition approach for analysis and forecasting of economic and financial time series[END_REF]. Meteorologists are interested in anomalous parts of weather data to predict future consequences [START_REF] Soares | Ensemble of evolving data clouds and fuzzy models for weather time series prediction[END_REF].

Consistent with [START_REF] Chandola | Comparative evaluation of anomaly detection techniques for sequence data[END_REF], we can define an anomaly as a point in time where the system's behavior is unusual and significantly different from previous normal behavior. An anomaly can mean a change with positive consequences, such as an increase in the number of sales on a sales website, and a change with negative consequences, such as a change in the rotation frequency of a jet engine's turbine. In either case, whether positive or negative, the similarity is that anomalies must often be considered and carefully taken into account.

Anomalies can be temporal if the temporal sequence of data is relevant; i.e., the data maybe anomalous only in a specific temporal context. Temporal anomalies are often subtle and hard to detect in real applications. Detecting temporal anomalies in applications is valuable as they can allow reactions (e.g. serve as an early warning for problems with the underlying system, or decision support for adjustments in order to benefit from the anomaly). This article aims to provide an organized overview of existing research to detect anomalies in time series. The objective is to provide an understanding of the problem of detecting anomalies and how existing techniques are related. The next section provides some definitions. In Section 3, a brief overview of time series anomaly detection techniques is reported.

Definitions

A time series is an ordered list of observations X =< x 1 , ..., x n >, where each value x p is an observation collected at time p. Let S p = x p , x p+1 + ... + x p+m-1 be the subsequence of size m, which starts at the position p of the time series X. Figure 1 presents an example of the time series, where O1 is an example of a subsequence.

Time-series patterns

We briefly describe some main properties of time series important for anomaly detection methods.

We can say that a time series has a trend if the observed mean µ is not constant, but varies over time. The trend can have a linear behavior. The time series in Figure 2 has a positive trend over the years.

The seasonality is the recurrence of oscillations periodically. A time series is seasonal if it is influenced by seasonal aspects such as year, period, or year. Thus, there is always a fixed-term, where the oscillation reappears. Figure 3 shows a seasonal time series with annual house sale values.The seasonality can be observed, since the real estate market is generally not active at the beginning of the year and sales generally increase in the middle of the year, decreasing again at the end of the year. A stationary time series is one whose properties do not depend on the instant at which the series is observed [START_REF] Hyndman | Forecasting : principles and practice[END_REF]. Thus, time series that have trends, or seasonality, are not stationary, as these properties affect the value of the time series at different times. On the other hand, a series of white noise is stationary because it doesn't matter when it's observed, it looks the same at all times.

The given definition assumes that a stationary series will have the following characteristics:

1. The verification of constant mean, shown that there is no trend in the series. 2. Throughout the series period, there is a constant autocorrelation 3. The series does not show periodic movements, indicating that there is no seasonality

Anomalies

Several studies have attempted to define the condition of anomalous data. In [START_REF] Hawkins | Identification of Outliers[END_REF], Hawkins defined anomaly as: "an observation that deviates so much from other observations that it raises the suspicion that a different mechanism generated it." In [START_REF] Pincus | Outliers in statistical data[END_REF], Barnett and Lewis present another interpretation of this definition: "an outlier is an observation (or subset of observations) that appears to be inconsistent with the rest of this dataset." We can then generalize the definition of anomalies as described in [START_REF] Chandola | Anomaly detection: A survey[END_REF] "patterns in data that do not conform to a welldefined notion of normal behavior".

Activities such as a system failure, error in the capture of information by a sensor, or even earthquakes are examples of reasons for anomalies in the data. Even in different scenarios, all of these anomalies produce important information for analysts. The importance of analyzing and understanding these occurrences is an important point for detecting anomalies.

Anomaly types

We can understand the problem of detecting anomalies for time series to identify outliers concerning some standard or usual signal. While there are many types of anomalies, let's focus on the most important types [START_REF] Chandola | Anomaly detection: A survey[END_REF].

1. Point anomaly: A point anomaly happens when an observation or a set of several individual observations diverges from the other set of observations. One can observe a punctual discontinuity in a period of time outside the normality of the values. For example, a purchase with a very high amount within a set of credit card transactions. This type is closest to the concept presented in [START_REF] Hawkins | Identification of Outliers[END_REF]. 2. Contextual anomaly: An apparently normal observation may diverge when analyzed within a specific context, such as a period of time or a certain location. For example, low temperatures are normal in the winter period but are contextual anomalous when they occur in the summer, where high temperatures are expected. 3. Collective anomalies: It occurs when no individual observation is an anomaly in itself, but a group of these observations, when combined, exhibit behavior that diverges from the rest of the data. In the stock market, for example, the fall in the price of an asset does not deviate significantly from the normal range, but the combination of successive falls indicates a collective anomaly.

The contextualization and understanding of the type of anomaly to be identified can help choose the best detection model.Each approach is aimed at certain types of anomalies, having certain advantages and disadvantages according to the subjective definition of normal observations and anomalies in a given context.

Anomaly Detection Several works provide an overview of techniques for detecting outliers [START_REF] Aggarwal | Outlier Analysis[END_REF]. Anomaly detection is also known as outlier detection, event detection, novelty detection, drift detection, change point detection, failure detection, or misuse detection [START_REF] Chandola | Anomaly detection: A survey[END_REF]. These different approaches and terms converge towards the same purpose: detecting abnormal patterns that deviate from the rest of the data, called anomalies or outliers.

Anomaly detection is a broad field that has been studied in the context of a large number of application domains, such as intrusion detection, fraud detection, failure detection in industry 3.0, system health monitoring, event detection in network sensors, and detection of disturbances in the ecosystem [START_REF] Thudumu | A comprehensive survey of anomaly detection techniques for high dimensional big data[END_REF]. Due to the great diversity of types and techniques, it is necessary to understand the properties of anomalies:

1. Temporality: anomalies can be temporal when associated with some temporal information. These are usually anomalies found in time series as well as streaming data and medical data. 2. Labeling: Some datasets have a set of labeled anomalies, where the location of each instance in the dataset is indicated. In general, this set of labeled anomalies represents only a subset of the anomalies present. Annotated anomalies are used by supervised learning methods, where the method learns through examples to detect new anomalies. 3. Dimensionality: anomalies can present univariate data when represented in only one dimension or multivariate when presenting a set of variables. Multivariate data is generated, for example, on sensors.

Anomaly Detection in Time Series Detecting anomalies in temporal data differs from detecting non-temporal data. In non-temporal data, as in spatial data, it is possible to detect an anomaly by its location in less dense or more distant regions. Another way is to calculate the deviation from the anomalous observations to the rest of the data. In this domain, it is understood that the observations are independent of each other.

When we look at the temporal data sets, we can see that this premise is not true, as the observations are not completely independent of each other. Previous observations may influence new observations. Thus, the variation between values tends to be smooth and gradual, without major variations.

Take, for example, a list of ten values with the pricing of the value of an asset in the stock market, measured every one minute: $12, $15, $16, $18, $17, $42, $43, $18, $19, $18. Since these points are dependent on each other, we can see that the sudden increase to the value of $42 could be an anomaly.

Anomaly detection methods for time series can be divided into two main categories [START_REF] Aggarwal | Outlier Analysis[END_REF]. There are methods that are based on time series prediction, which represent most of the statistical methods [START_REF] Thudumu | A comprehensive survey of anomaly detection techniques for high dimensional big data[END_REF]. Another category are methods that are based on unusual forms of the time series. We will present some of these approaches in the next section.

Anomaly Detection Approaches For Time Series

In this section, we discuss various anomaly detection applications. We organize the approaches into three categories: statistics-based, cluster-based, and matrix profile-based.

Statistical Based Approaches

We present some stationary linear autoregressive models. Most of the methods used in time series are linked to linearity and stationarity. A process is stationary when its mean, variance and autocovariance are invariant with respect to time, such as, for example, the Auto-Regressive (AR) and Moving Averages (MA) and Auto-Regressive Moving Averages (ARMA) models.

Autoregressive Model (AR) In an autoregressive model, denoted by AR, we project the variable of interest based on the linear combination of a finite set of values before to the variable (independent variables) and an error value [START_REF] Hyndman | Forecasting : principles and practice[END_REF].

Thus, an autoregressive model of order p, AR(p), can be written as:

Z t = ϕ 1 Z t-1 + ϕ 2 Z t-2 + . . . + ϕ p Z t-p + a t (1)
where a t is white noise and the terms Z t-1 , Z t-2 , . . . , Z t-p are terms that are independent of a t . This model assumes that the current value of the series is a linear combination of the past p values of the series and a white noise a t . We refer to Equation 1 as an AR(p) model, an autoregressive model of order p.

Autoregression-based techniques were extended to detect contextual anomalies in time series. Initially, the AR(p) model is fitted to the data. Then, for each Z t instance, the residual of the instance is determined for calculating the anomaly score. This residual is the instance value that falls outside the regression model. Thus, the anomaly is scored as the difference between the estimated value and the residual value [START_REF] Chandola | Anomaly detection: A survey[END_REF].

Moving Average Model (MA) Instead of using, as in Autoregressive Model (AR), the p observations prior to the variable of interest, a Moving Average Model (MA) considers the last q errors t , t-1 , . . . , t-q prediction in a regression model. We refer to this model as an MA (q), a moving average model of order q [13]:

Z t = c + t + θ 1 t-1 + θ 2 t-2 + . . . + θ q t-q (2)
where the θ 1 , θ 2 , . . . , θ q are the parameters of the model, c is the expectation of Z t (often assumed to equal 0), and the t , t-1 , . . . , t-q are the white noise error terms.

It is important to note that this is a moving average model used to predict future values, unlike smoothing models, where average values are used to eliminate some of the randomnesses of the data. Therefore, we consider that each Z t comes to be considered a weighted moving average of the latest forecast errors.

In the model of Equation 2, the moving average of the previous q observations, also called the MA of order q, is considered to calculate the current estimated value. The next step includes checking that the estimated value is within the predefined confidence band. Confidence band is the interval defined as a multiple of the standard deviation of the moving average of the previous period. If the value is greater than the maximum expected value (higher confidence band), it is then flagged as an anomaly.

Autoregressive Moving Average Model (ARMA) Autoregressive Moving Average (ARMA) models play a key role in time series modeling, being widely used for time series analysis and prediction. The ARMA model is a combination of the Autoregressive (AR) model, which describes the analytical component of the signal, and the Moving Average (MA) model, which describes the noise component of the signal. Its linear structure also presents a substantial simplification of linear prediction. Compared to pure AR or MA models, ARMA models provide the most efficient linear model in stationary time series, given the ability to model the unknown process with minimal parameterization [START_REF] Wang | Real-time prediction and anomaly detection of electrical load in a residential community[END_REF].

An autoregressive and moving average (ARMA) model denoted by ARMA(p, q), combining the AR(p) and MA(q) models, can be written in a single equation in the form:

Z t = c + t + θ 1 t-1 + θ 2 t-2 + . . . + θ q t-q + ϕ 1 Z t-1 + ϕ 2 Z t-2 + . . . + ϕ p Z t-p (3)
Let Z t be the sign Z at time t. We assume that Z t linearly depends on the previous values Z t-1 , . . . , Z t-p where p is the order of the autoregression. Thus, we have from AR that ϕ 1 , . . . , ϕ p are auto-regression parameters that can be learned from historical data and used to predict or find similar time series, c is a constant, and t N (0, σ 2 ) is Gaussian noise. From MA we have θ 1 , . . . , θ q are learnable parameters of the model, µ is the expectation of Z t and ti N (0, σ 2 ) are the terms of Gaussian noise.

Using historical data, we can select p and q in ARMA(p,q) model and learn the model coefficients θ and φ based on which we can make a future prediction. A substantial deviation from the prediction result in a time-series anomaly. We can define substantial as two standard deviation from a moving average of Z. The model parameteters θ and φ are learned via Maximum Likelihood Estimation (MLE). Given the simplicity of the model it may be sufficient for many applications.

The ARMA model has the advantage of having a clear mathematical and statistical basis. However, there are some disadvantage, like the selection of the best values for p and q in order to find a better model for detecting anomalies [START_REF] Braei | Anomaly detection in univariate time-series: A survey on the state-of-the-art[END_REF]. If high values are used, the generated model will find a large number of false-negatives, thus, a low number of anomalies. In contrast, using small values for p and q, the generated model will identify a high number of false positives, that is, label a large number of observations as anomalies, which actually are not. There are some ways to adjust the model, such as the use of correlograms, the use of cross-validation [START_REF] Aggarwal | Outlier Analysis[END_REF] and the Box-Jenkins method [START_REF] Box | Time series analysis : forecasting and control[END_REF].

Clustering-Based Approaches

Clustering-based approaches are used to group similar datasets into clusters. Although clustering and anomaly detection tasks are different in nature, sev-eral clustering-based anomaly detection techniques have been developed. During the clustering process, these approaches consider the identified observations and clusters and then detect the anomalies [START_REF] Chandola | Anomaly detection: A survey[END_REF].

We can group these techniques into two groups based on different assumptions. In first, the dataset is grouped into similar data clusters. Instances that do not belong to any cluster are marked as anomalies. Several algorithms do not force all instances to belong to a cluster; thus, identified clusters are removed from the dataset, and residual instances are noted as anomalies. A disadvantage of these techniques is that they are not optimized to identify anomalies, as their initial objective is to identify clusters. Among these techniques, we have DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] as the most used.

In the second group of techniques, initially, the dataset into clusters, using a clustering algorithm. Then, the distance of each instance to the centroid of the nearest cluster is calculated. Instances that are further away from your centroid are marked as anomalies. A series of anomaly detection techniques that follow this approach have been proposed using different clustering algorithms, among these, the most notorious K-Means clustering [START_REF] Likas | The global k-means clustering algorithm[END_REF].

Density-Based Spatial Clustering of Applications with Noise (DB-SCAN) Density-based spatial clustering of applications with noise (DBSCAN)

is a data clustering non-parametric algorithm [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. The method is significantly effective in identifying clusters of arbitrary shape and sizes, identifying and separating noise from the data, and detecting "natural" clusters and their arrangements within the data space, without any preliminary information about the groups.

Given a set of points in some space, DBSCAN groups points that are intimately very close together (points with many nearby neighbors), marking them as outliers (anomalies) points that are isolated in low-density regions (whose closest neighbors are very far away) [START_REF] Celik | Anomaly detection in temperature data using DBSCAN algorithm[END_REF]. Two important user-defined parameters are required: neighborhood distance epsilon (eps) and a minimum number of points minpts. For a given point, the points in the eps distance are called neighbors of that point. If the total number of points neighboring a point is greater than minpts, this group is called a cluster.

DBSCAN labels the data points in three categories: 1) core: the points that have at least minpts number of points in the eps distance; 2) border: those that are not core points but are the neighbors of core points; 3)outlier (anomalous): those that do not belong to any cluster. Emadi et al. [START_REF] Emadi | A novel anomaly detection algorithm using DB-SCAN and SVM in wireless sensor networks[END_REF] propose an algorithm based on DBSCAN to detect anomalies in Wireless Sensor Networks.

K-Means Clustering

One of the most explored techniques for grouping data is the K-means [START_REF] Likas | The global k-means clustering algorithm[END_REF]. It is a cluster based on centroids that partitions the dataset into k clusters of similar instances. We can define the k-means algorithm in a sequence of steps. Initially, the number K of clusters is defined. Then, the k centroids are initialized, which can be done by arbitrary choice. Then, for each object, the distance to the centroid of all clusters and connected to the nearest centroid is calculated. Then the centroids of the modified clusters are recalculated. Finally, the distance from the objects to the centroids is recalculated, updating the connections with the nearest ones. This last step is repeated until there are no more updates [START_REF] Kumari | Anomaly detection in network traffic using k-mean clustering[END_REF].

Originally, the k-means approach was not defined to work with time series. One of the possibilities to use k-means as anomaly detection in time series is through the use of the sliding window [START_REF] Li | Clustering-based anomaly detection in multivariate time series data[END_REF] approach. With this approach, a set of subsequences of equal lengths of the same time series is generated. In this set, the k-means algorithm is applied until it converges on the searched k clusters. To perform the anomaly search, the distance between each subsequence and its associated centroid is computed. This subsequence is marked as an anomaly if this distance value is greater than a defined threshold value δ. One of the biggest challenges of this approach is the correct parameterization of the k amounts of clusters.

Matrix Profile Technique

In 2016, Yeh et al. [START_REF] Yeh | Matrix profile i: All pairs similarity joins for time series: A unifying view that includes motifs, discords and shapelets[END_REF] published a novel technique to perform all-pairs-similaritysearch on two time-series, producing two new series: the Matrix Profile and the Matrix Profile Index. The Matrix Profile is defined as a data structure containing the z-normalized Euclidean distances between each subsequence of the first series and its closest corresponding subsequence of the second time series. The Array Profile Index contains the index of the closest matching substring in the second series for each substring. By itself, the Matrix Profile can be used to detect anomalies in contexts where anomalies are defined by unique behavior [START_REF] Paepe | A generalized matrix profile framework with support for contextual series analysis[END_REF]. In fact, in the Matrix Profile vector the anomalies can be detected in the points with high values, because the distance of the subsequences represented by these points to their closest matching subsequence is high.

In general, given two series of n real values, S1 ∈ R n and S2 ∈ R n and a subsequence length m, the Matrix Profile M ∈ R n-m+1 and a Matrix Profile Index I ∈ R n-m+1 are new series such that for each i ∈ [0, n-m], I i contains the index of the start of the subsequence of S2 with length m that best matches S1 i,m and M i contains the corresponding distance. In the case a self-join is performed where S1 = S2, an additional constraint is added to prevent trivial matches, where subsequences match themselves or nearby subsequences, called exclusion zone. The default distance measure used is the z-normalized Euclidean distance, which removes the effect of a changing data offset over time and thus focuses more on shape instead of amplitude. Typical causes of a changing offset are wandering baselines in sensors or natural phenomena (e.g., the gradual change in temperature throughout seasons) [START_REF] Paepe | A generalized matrix profile framework with support for contextual series analysis[END_REF]. Figure 4 shows an example of the application of the matrix profile in a time series for the detection of anomalies. In this figure, the top image represents a time series, and the bottom image its matrix profile. The top 3 anomalies are marked in red, i.e., the three points where the matrix profile has high values.

Matrix profile is usually able to detect subsequences with unusual shapes in the data. It can be used for detecting point anomalies and also collective anomalies (i.e., a sequence of abnormal points). In the rest of this section, we present the main algorithms for calculating the Matrix Profile.

STAMP The Matrix Profile was originally published together with the STAMP (Scalable Time Series Anytime Matrix Profile) [START_REF] Yeh | Matrix profile i: All pairs similarity joins for time series: A unifying view that includes motifs, discords and shapelets[END_REF], an anytime algorithm to calculate the Matrix Profile over a time series and the corresponding Index. Internally, STAMP uses a similarity search algorithm called MASS [START_REF] Mueen | The fastest similarity search algorithm for time series subsequences under euclidean distance[END_REF] that under z-normalized Euclidean iteratively calculates the distance profile of each subsequence, which is the distance of the subsequence to every subsequence, by using the Fast Fourier Transform (FFT).

The STAMP is outlined in Algorithm 1. In line 2, the length of T B is extracted. In line 3, the matrix profile P AB and matrix profile index I AB are initialized. From lines 4 to line 6, the distance profiles D are calculated, using each subsequence B[idx] in the time series T B and T A . The pairwise minimum for each element in D is performed with the paired element in P AB (i.e., min(D[i], P AB [i]) for i = 0 to length(D) -1). Then, as the minimum pair operations are performed, I AB [i] is updated with idx (when D[i] ≤ P AB [i]). Finally, the results, i.e., P AB and I AB , are returned in line 7. In this format, STAMP computes the matrix profile for the general similarity join. It is possible to change the algorithm to compute the self-similarity join matrix profile of a time series T A , just by replacing T B in line 2 with T A , replace B with A in line 5, and ignore trivial matches in D when performing ElementWiseMin in line 6.

The overall complexity of the algorithm is O(n 2 log n) where n is the length of the time series. Since all subsequences are compared using the MASS algorithm, the n log n factor comes from the FFT subroutine. STOMP The STOMP algorithm is similar to STAMP [15] in that it can be seen as highly optimized nested loop searches, with the repeated calculation of distance profiles as the inner loop. However, while STAMP must evaluate the distance profiles in random order (to allow its anytime behavior), STOMP performs an ordered search. It is by exploiting the locality of these searches that STOMP can reduce the time complexity by a factor of O(log n). STOMP uses the z-normalized Euclidean distance d i,j , as shown below, of two time series subsequences T i,m and T j,m using their dot product, QT i,j :

d i,j = 2m 1 - QT i,j -mµ i µ j mσ i σ j (4)
Here m is the subsequence length, µ i is the mean of T i,m , µ j is the mean of T ( j, m), σ i is the standard deviation of T i,m , and σ j is the standard deviation of T j,m . Note that QT i,j can be decomposed as:

QT i,j = m-1 k=0 T i+k T j+k (5)
The time required to compute d i,j depends only on the time required to compute QT i,j . To solve this problem, STOMP pre-computes and stores the means and standard deviation in O(n) space and time, thus, it takes O(1) to compute d i,j [START_REF] Zhu | Matrix profile ii: Exploiting a novel algorithm and gpus to break the one hundred million barrier for time series motifs and joins[END_REF].

The pseudo-code of STOMP algorithm is shown in Algorithm 2. It begins in line 1 by computing the matrix profile length l. In line 2, it calculates the mean and standard deviation of every subsequence in T . Line 3 calculates the first dot product vector QT with the algorithm in TABLE I. Line 5 initializes the matrix profile P and matrix profile index I. The loop in lines 6-13 calculates the distance profile of every subsequence of T in sequential order, with lines 7-9 updating QT according to [START_REF] Celik | Anomaly detection in temperature data using DBSCAN algorithm[END_REF]. Then update QT [START_REF] Aggarwal | Outlier Analysis[END_REF] in line 10 is done with the pre-computed QT f irst in line 3. Line 11 calculates distance profile D according to Equation (4). Finally, line 12 compares every element of P with D: if D[j] < P [j], then P

QT [j] ← QT [j -1] -T [j -1] × T [i -1] + T [j + m -1] × T [i + m -1] 10 QT [1] ← QT f irst [i]
[j] = D[j], I[j] = i.
The time complexity of STOMP is O(n 2 ). Thus, it can achieve a O(log n) factor speedup over STAMP [15]. The O(log n) speedup makes little difference for small datasets, however, when considering the datasets with millions of data points, this O(log n) factor begins to produce a significant performance gain. SCRIMP++ An extension of STOMP is proposed in [START_REF] Zhu | Matrix profile xi: Scrimp++: Time series motif discovery at interactive speeds[END_REF]. SCRIMP is an anytime algorithm that computes the matrix profile algorithm combining the anytime component of STAMP with the speed of STOMP. The optimization of SCRIMP++ is performed using the incremental calculation of the D diagonals of the scalar product of Equation 4, as follows:

Q i,j = Q i-1,j-1 -t i-1 t j-1 + t i+m-1 t j+m-1 (6)
Equation 6 presents the incremental approach, where the values of the diagonal cells can be calculated using the cell value previously calculated. This approach reduces the number of operations required for the new calculation.

Algorithm 3 shows the pseudo-code of SCRIMP++. Line 2 precomputes the means and standard deviations of all subsequences in T. In line 4, matrix profile P and matrix profile index I are initialized. In lines 6-16, the diagonals of the returnP, I distance matrix are iteratively evaluated, being chosen in random order. Figure 5 shows an example of this process. Diagonal distance values, such as d 1,k , d 2,k , . . . , d n-m+2-k,n-m+1 are calculated one by one. If d i,i+k-1 , for any i < n-m+1, referenced by d in line 12, is less than P i (line 13) or P i+k-1 (line 15), then the associated matrix profile and index values are updated. The iterative algorithm can be interrupted by the user to analyze the matrix profile and index values.

AAMP In many applications it is preferred to not normalize the time-series data because the anomalies can be deteted based on the point values, and not the shapes. AAMP [START_REF] Akbarinia | Efficient matrix profile computation using different distance functions[END_REF] has been designed for such applications. It is an efficient algorithm for computing matrix profile with the pure (non-normalized) Euclidean distance. AAMP is executed in a set of iterations, such that in each iteration the distance of subsequences is computed incrementally. The time complexity of AAMP is O(n × (n -m)) with small constants, where n is the time series length and m the subsequence length. The experiments reported in [START_REF] Akbarinia | Efficient matrix profile computation using different distance functions[END_REF] show that the performance of AAMP is significantly better than that of STAMP and SCRIMP++ (an improved version of STOMP).

The main idea behind AAMP is that for computing the distance between subsequences it uses diagonal sliding windows, such that in each sliding window, the Euclidean distance is incrementally computed only between the subsequences that have a precise difference in their start position. Let T i = t i , t i+1 , . . . , t i+m-1 and T j = t j , t j+1 , . . . , t j+m-1 be two subsequences. The sliding windows in Fig. 5. A SCRIMP++ iteration evaluates a randomly selected diagonal, thus updating the matrix profile AAMP allow to use Equation 7for incremental computation of the distance between subsequences T i and T j (denoted by D i,j ) by using the yet computed distance between subsequences T i-1 and T j-1 (denoted as D i-1,j-1 ):

[H]D i,j = D 2 i-1,j-1 -(t i-1 -t j-1 ) 2 + (t i+m-1 -t j+m-1 ) 2 Algorithm 4 shows the pseudo-code of AAMP. Initially, the algorithm sets all the values of the matrix profile to infinity (i.e., maximum distance). Then, it performs n -m -1 iterations using a variable k (1 ≤ k ≤ n -m -1). In each iteration k, the algorithm compares each subsequence T i,m with the subsequence that is k positions far from it, i.e., T i,m+k . To do this, AAMP firstly computes the Euclidean distance of the first subsequence of the time series, i.e., T 1,m , with the one that starts at position k, i.e., T k,m . This first distance computation is done using the normal formula of Euclidean distance. Then, in a sliding window, the algorithm incrementally computes the distance of other subsequences with the subsequences that are k position far from them, and this is done by using Equation 7in O(1). If the computed distance is smaller than the previous minimum distance that is kept in the matrix profile P , then it is updated with this new lower value.

Conclusion

In this paper, we presented a survey of anomaly detection methods in time series datasets. We firstly presented the main concepts related to anomalous data in different applications, and then defined the anomaly detection task. Afterwards, we described the important approaches for anomaly detection in three main categories: statistical based, clustering based, and matrix profile based.

Fig. 1 .

 1 Fig. 1. Time series

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sample time-series showing a constant trend

Fig. 4 .

 4 Fig. 4. Matrix Profile Anomaly Detection Example

Algorithm 1 : 2 nB ← Length(TB) 3 PAB 4 for each idx in idxes do 5 D

 12345 STAMP (T A , T B , m) Input: Two time series, TA and TB Interested subsequence length m Output: A matrix profile PAB and associated matrix profile index IAB of TA join TB 1 begin ← infs, IAB ← zeros, idxes ← 1 : nB -m + 1 ← MASS(B[idx], TA) 6 PAB, IAB ← ElementWiseMin(PAB, IAB, D, idx) 7 return PAB, IAB

Algorithm 2 : 2 n 4 QT 5 D 6 P 7 for i=2 to l do 8 for j=l downto 2 do 9

 22456789 STOMP (T, m)Input: A time series T and a subsequence length m Output: Matrix profile P and the associated matrix profile indexI of T 1 begin ← Length(T ), l ← n -m + 1 3 µ, σ ← ComputeMeanStd (T, m) ← SlidingDotProduct(T [1 : m], T ), QT f irst ← QT ← CalculateDistanceProfile (QT, µ, σ) ← D, I ← ones // initialization

11 D

 11 ← CalculateDistanceProfile (QT, µ, σ, i) 12 P, I ← ElementWiseMin (P, I, D, i) 13 return P, I

Algorithm 4 : 2 for i=1 to n do 3 P

 423 AAMP algorithmInput: T : time series; n: length of time series; m: subsequence length Output: P : Matrix profile;1 begin [i] = ∞ ; 4 for k=1 to n-m-1 do 5 dist = Euc Distance(T1,m, T k,m ) 6 if dist < P[1] then 7 P[1] = dist; 8 if dist < P[k] then 9 P[k] = dist; 10 for i=2 to n -m + 1 -k do 11 dist = (dist 2 -(ti-1 -t i-1+k ) 2 + (ti+m-1 -t i+m+k-1 ) 2 12 if dist < P[i] then 13 P[i] = dist;14 if dist < P[i+k] then 15 P[i+k] = dist;