tsmp: An R Package for Time Series with Matrix Profile

This article describes tsmp, an R package that implements the matrix profile concept for time series. The tsmp package is a toolkit that allows all-pairs similarity joins, motif, discords and chains discovery, semantic segmentation, etc. Here we describe how the tsmp package may be used by showing some of the use-cases from the original articles and evaluate the algorithm speed in the R environment. This package can be downloaded at https://CRAN.R-project.org/package=tsmp.


The matrix profile
The reader may be aware of what a DM is. It is widely used in TS for clustering, classification, motif search, etc. But, even for modestly sized datasets, the algorithms can take months to compute even with speed-up techniques such as indexing (Shieh and Keogh, 2008;Fu et al., 2008), lower-bounding (Keogh and Ratanamahatana, 2005), data discretization (Lin et al., 2003) and early abandoning (Faloutsos et al., 1994). At best, they can be one or two orders of magnitude faster.
The MP is an ordered vector that stores the Euclidean distance between each pair within a similarity join set. One (inefficient) way would be to use the full DM of every iteration of a sliding window join and retrieve just the smallest (non-diagonal) value of each row. The MP also has a companion vector called PI, that gives us the position of the nearest neighbour of each subsequence. This method has a host of interesting and exploitable properties. For example, the highest point on the MP corresponds to the TS discord, the (tied) lowest points correspond to the locations of the best TS motif pair, and the variance can be seen as a measure of the TS complexity. Moreover, the histogram of the values in the MP is the exact answer to the TS density estimation. Particularly, it has implications for TS motif discovery, TS joins, shapelet discovery (classification), density estimation, semantic segmentation, visualisation, rule discovery, clustering, etc. (Yeh et al., 2017b).
Some of the advantages/features of this method: • It is exact, providing no false positives or false dismissals.
• It is simple and parameter-free. In contrast, the more general metric space APSS algorithms require building and tuning spatial access methods and/or hash functions.
• It requires an inconsequential space overhead, just O(n) with a small constant factor.
• It is extremely scalable, and for massive datasets, we can compute the results in an anytime fashion, allowing ultra-fast approximate solutions.
• Having computed the similarity join for a dataset, we can incrementally update it very efficiently. In many domains, this means we can effectively maintain exact joins on streaming data forever.
• It provides full joins, eliminating the need to specify a similarity threshold, which is a nearimpossible task in this domain.
• It is parallelizable, both on multicore processors and in distributed systems (Zhu et al., 2016).

The tsmp package
The tsmp package provides several functions that allow for an easy workflow using the MP concept for TS mining. The package is available from the CRAN at https://CRAN.R-project.org/package=tsmp. In Section 2.3.1 we explain how to install this package. In Section 2.3.2 we describe the syntax for the main functions in tsmp, giving an example of a particular model. In Section 2.3.3 we will further explain the available algorithms for MP computation and its current use. In Section 2.3.4 we show some examples of MP application for data mining.

Input arguments and example
The tsmp has a simple and intuitive workflow. First, you must compute the MP of the desired TS. Depending on the task, the user might want to follow one of three paths: univariate self-join, AB-join or multivariate self-join. One exception is the SiMPle algorithm that is a multivariate AB-join and will be explained in Section 2.3.3.
The main function is tsmp(), which has the following usage: tsmp(..., window_size, exclusion_zone = 1/2, mode = c("stomp", "stamp", "simple", "mstomp", "scrimp"), verbose = 2, s_size = Inf, must_dim = NULL, exc_dim = NULL, n_workers = 1, .keep_data = TRUE) The first argument ellipsis (the three dots) receives one or two TS. For self-joins, the user must input just one TS; two for AB-joins. Multivariate TS may be input as a matrix where each column represents one dimension. Alternatively, the user may input the Multivariate TS as a list of vectors. The second argument window_size is the size of the sliding window. These are the most basic parameters you need to set.
Further parameters are: • exclusion_zone, is an important parameter for self-joins. This is used to avoid trivial matches and is a modifier of the window_size, i.e., for an exclusion_zone of 1/2, and window_size of 50, internally the result will be 25.
• mode, here the user may choose the algorithm used for the MP calculation. stomp, stamp and scrimp return equal results, although differing in some practical attributes, and they will be further explained in Section 2.3.3. mstomp is designed for Multivariate TS self-join only. simple is designed for Multivariate TS for self-join and AB-join, which will also be further explained in Section 2.3.3.
• verbose, controls the verbosity of the function. 0 means no feedback, 1 means text messages only, 2 (the default) means text messages and progress bar, and 3 also plays a sound when finished.
• s_size, controls the anytime algorithms. This is just a way to end the algorithm in a controlled manner because the anytime algorithms can be stopped anytime and the result will be returned.
• must_dim, is an optional argument for the mstomp algorithm. See next item.
• exc_dim, as must_dim, is an optional argument for the mstomp algorithm. These arguments control which dimensions must be included and which must be excluded from the multidimensional MP.
• n_workers, controls how many threads will be used for the stamp, stomp, and mstomp. Note that for small datasets, multiple threads add an overhead that makes it slower than just one thread.
• .keep_data, TRUE by default, keeps the input data inside the output object. This is useful for chained commands.

Example data
We think that the best and simple example to demonstrate the tsmp package is the motif search.
The tsmp package imports the %>% (pipe) operator from the magrittr package that makes the tsmp workflow easier.
The following code snippet shows an example of the workflow for motif search: R> data <-mp_fluss_data$walkjogrun$data R> motifs <-tsmp(data, window_size = 80, exclusion_zone = 1/2) %>% + find_motif(n_motifs = 3, radius = 10, exclusion_zone = 20) %T>% plot() The find_motif() function is an S3 class that can receive as the first argument the output of tsmp() function as a univariate or multivariate MP. This allows us to use the pipe operator easily. The plot() function is also an S3 class extension for plotting objects from the tsmp package and works seamlessly.

Computational methods
There are several methods to compute the MP. The reason for that is the unquenchable need for speed of the UCR's researchers. Before starting, let's clarify that the time complexity of a brute force algorithm has a time complexity of O(n 2 m), for n being the length of the reference TS and m the length of the sliding window (query) that is domain dependent.

STAMP
This was the first algorithm used to compute the MP. It uses the MASS (Mueen et al., 2015) as the core algorithm for calculating the similarity between the query and the reference TS, called the DP. The ultimate MP comes from merging the element-wise minimum from all possible DP. This algorithm has the time complexity of O(n 2 log n) and space complexity of O(n) (Yeh et al., 2017b). The anytime property is achieved using a random approach where the best-so-far MP is computed using the DP that have been already calculated.

STOMP
This was the second algorithm used to compute the MP. It also uses the MASS to calculate the DP but only for the first iteration of each batch. The researchers noticed that they could reuse the values calculated of the first DP to make a faster calculation in the next iterations. This results on a time complexity of O(n 2 ), keeping the same space complexity of O(n). This algorithm is also suitable for a GPU framework (although this was not yet implemented in tsmp package) (Zhu et al., 2016). The main drawback of STOMP compared with STAMP is the lack of the anytime property. In scenarios where a fast convergence is needed (e.g., finding the top-k motifs) it may be required only 5% of the MP computation to provide a very accurate approximation of the final result.

SCRIMP
The SCRIMP algorithm is still experimental at the time of this article. It combines the best features of STOMP and STAMP, having a time complexity of O(n 2 ) and the anytime property (UCR, 2016).

SiMPle
The SiMPle algorithm is a variation designed for music analysis and exploration (Silva et al., 2018). Internally it uses STOMP for MP computation and allows multidimensional self-joins and AB-joins. The resulting MP is computed using all dimensions. One major difference is that it doesn't apply any z-normalization on the data, since for music domain this would result in spurious similarities.

mSTOMP
The mSTOMP algorithm was designed to motif search in multidimensional data (Yeh et al., 2017a). Performing motif search on all dimensions is almost guaranteed to produce meaningless results, so this algorithm, differently from SiMPle, doesn't compute the MP using all dimensions naïvely, but the d-dimensional MP for every possible setting of d, simultaneously, in O(dn 2 log d) time and O(dn) space. The resulting MP allow motif search in multiple dimensions and also to identify which dimensions are relevant for the motifs founded.

Motif search
In Section 2.3.2 we have shown a basic example of the workflow for motif search. Let's take a look at the result of that code: As we can see, this is a summary that tsmp package automatically generates from the resulting object. One nice property is that the object always holds the original MP and by default also holds the input data so that you can keep mining information from it. If the dataset is too big or you are concerned about privacy, you may set the argument .keep_data = FALSE.

R> motifs
In addition to this summary, you can see the results using plot() in Figure  This dataset is the WalkJogRun PAMAP's dataset (Reiss and Stricker, 2012). It contains the recording of human movements in three states, walking, jogging and running. As we can see, the plot shows the motifs of each state. Experienced readers might say that this is not the purpose of motif search, and we agree. The result shown here was achieved using a large radius and exclusion_zone to force the algorithm to look for distant motifs. Semantic segmentation is the proper algorithm for this task, and we will show this in the next section.

Semantic segmentation
As previously explained, the resulting object holds the original data and MP. So let's save some time and use the resulting object from the last section to try to find where the human subject started to jog and to run: We can see that this object now holds information of the FLUSS algorithm (Gharghabi et al., 2017), but the motif information is still there and can be retrieved using as.motif(). In Figure 2 we can see the graphic result of the segmentation. R> plot(segments, type = "data")  is 3800 and 6800). The middle graphic shows the data. The lower graphic shows the normalised arc counts with correction for the "edge-effect" (Gharghabi et al., 2017).

Time series chains
As a final example of practical application, let's search for a new kind of primitive: time series chains (Zhu et al., 2018a). This algorithm looks for patterns that are not just similar but evolve through time. The dataset used in this example is a record of the Y-axis of a mobile phone accelerometer while placing it on a walking subject's pocket (Hoang et al., 2015). The authors of this dataset wanted to analyse the stability of the mobile phone as it slowly settles in the pocket. This is a good example of a pattern that changes through time. Let's start with the workflow for this example: Here we see that the algorithm found 58 chains. Id est, it found 58 evolving patterns with at least three elements, and the best one is presented in the last line, a chain with six elements. Figure 3 shows the patterns discovered.

Speed
While this new method for TS data mining is extremelly fast, we have to take into consideration that the R environment is not as fast as a low-level implementation such as C/C++. In Table 1 we present the comparison to the MATLAB version that is available at the UCR. Yeh et al. (2017b) shows that the slowest algorithm (STAMP) can be hundreds of times faster than the MK algorithm (the fastest known exact algorithm for computing TS motifs) (Yoon et al., 2015), while the R implementation is just 1.65 to 8.04 times slower than MATLAB's, which is not a problem for an R researcher.

Conclusion
The examples in Section 2.3.4 show how straightforward the usage of tsmp package is. Regardless, these examples are just a glimpse of the potential of the MP. Several new algorithms based on MP are being developed and will be gradually implemented in the tsmp package (Linardi et al., 2018;Zhu et al., 2018b;Gharghabi et al., 2018;. Yeh et al. (2017a) for example, have developed an algorithm to allow MDS visualisation of motifs. Gharghabi et al. (2018) have developed a new distance measure that better suits repetitive patterns .
The MP has the potential to revolutionise the TS data mining due to its generality, versatility, simplicity and scalability (UCR, 2016). All existing algorithms for MP have been proven to be flexible to be used in several domains using very few parameters and they are also robust, showing good performance with dimensionality reduced data and noisy data. In addition, a yet to be published article shows a fantastic score of > 10 18 pairwise comparisons a day using GPU for motif discovery (Zimmerman et al., 2018).
The tsmp package is the first known MP toolkit available on any statistical language, and we hope it can help researchers to better mining TS and also to develop new methods based on MP.