An efficient aggregation method for the symbolic representation of temporal data

ABBA is a symbolic representation of time series based on an adaptive polygonal chain approximation and a mean-based clustering algorithm. Formally, we consider the problem of aggregating a time series $T=[t_{0},t_{1},\ldots,t_{N}]\in\mathbb{R}^{N+1}$ into a symbolic representation $S=[s_{1},s_{2},\ldots,s_{n}]\in\mathbb{A}^{n}$, $n\ll N$, where each $s_{j}$ is an element of an alphabet $\mathbb{A}=\{a_{1},a_{2},\ldots,a_{k}\}$ of $k$ symbols. The notation is also summarized in Table 1. The ABBA method consists of two steps, compression and digitization, which we review below.

The ABBA compression step computes an adaptive piecewise linear continuous approximation of $T$. Given a compression tolerance tol, the method adaptively selects $n+1$ indices $i_{0}=0<i_{1}<\cdots<i_{n}=N$ so that the time series $T=[t_{0},t_{1},\ldots,t_{N}]$ is approximated by a polygonal chain going through the points $(i_{j},t_{i_{j}})$ for $j=0,1,\ldots,n$. This gives rise to a partition of $T$ into $n$ pieces $P_{j}=[t_{i_{j-1}},t_{i_{j-1}+1},\ldots,t_{i_{j}}]$, each of integer length $\texttt{len}_{j}:=i_{j}-i_{j-1}\geq 1$ in the time direction. The criterion for the partitioning is that the squared Euclidean distance of the values in $P_{j}$ from the straight polygonal line is bounded by $(\texttt{len}_{j}-1)\cdot\texttt{tol}^{2}$. More precisely, starting with $i_{0}=0$ and given an index $i_{j-1}$, we find the largest possible $i_{j}$ such that $i_{j-1}<i_{j}\leq N$ and

$\displaystyle\sum_{i=i_{j-1}}^{i_{j}}\Big{(}\,\underbrace{t_{i_{j-1}}+(t_{i_{j}}-t_{i_{j-1}})\cdot\frac{i-i_{j-1}}{i_{j}-i_{j-1}}}_{\text{straight line approximation}}\ \ -\underbrace{t_{i}}_{\text{actual value}}\!\!\!\!\!\!\Big{)}^{2}\leq(i_{j}-i_{j-1}-1)\cdot\texttt{tol}^{2}.$

(1)

Each linear piece $P_{j}$ of the resulting polygonal chain $\widetilde{T}$ is described by a tuple $(\texttt{len}_{j},\texttt{inc}_{j})$, where $\texttt{inc}_{j}=t_{i_{j}}-t_{i_{j-1}}$ is the increment in value. The whole polygonal chain can be recovered exactly from the first value $t_{0}$ and the tuple sequence

$$(\texttt{len}_{1},\texttt{inc}_{1}),\ldots,(\texttt{len}_{n},\texttt{inc}_{n})\in\mathbb{R}^{2}.$$

(2)

The ABBA digitization step assigns the tuples in (2) into $k$ clusters $S_{1},S_{2},\ldots,S_{k}$. Before doing so, the tuple lengths and increments are separately normalized by their standard deviations $\sigma_{\texttt{len}}$ and $\sigma_{\texttt{inc}}$, respectively. A further scaling parameter scl is used to assign different weight (“importance”) to the length of each piece in relation to its increment value. Hence, one effectively clusters the scaled tuples

$$p_{1}=\left(\texttt{scl}\frac{\texttt{len}_{1}}{\sigma_{\texttt{len}}},\frac{\texttt{inc}_{1}}{\sigma_{\texttt{inc}}}\right),\ldots,p_{n}=\left(\texttt{scl}\frac{\texttt{len}_{n}}{\sigma_{\texttt{len}}},\frac{\texttt{inc}_{n}}{\sigma_{\texttt{inc}}}\right).$$

(3)

If $\texttt{scl}=0$, then clustering is performed on the increments alone, while if $\texttt{scl}=1$, the lengths and increments are clustered with equal weight.

In the original ABBA method [15] the cluster assignment is performed using k-means [34], by iteratively aiming to minimize the within-cluster-sum-of-squares

$$\texttt{WCSS}=\sum_{i=1}^{k}\sum_{(\texttt{len},\texttt{inc})\in S_{i}}\Bigg{\|}\left(\texttt{scl}\frac{\texttt{len}}{\sigma_{\texttt{len}}},\frac{\texttt{inc}}{\sigma_{\texttt{inc}}}\right)-\overline{\mu}_{i}\Bigg{\|}^{2},$$

(4)

with each 2d cluster center $\overline{\mu}_{i}=(\overline{\mu}_{i}^{\texttt{len}},\overline{\mu}_{i}^{\texttt{inc}})$ corresponding to the mean of the scaled tuples associated with the cluster $S_{i}$.

Given a partitioning of the $n$ tuples into clusters $S_{1},\ldots,S_{k}$, the unscaled cluster centers $\mu_{i}$

$$\mu_{i}=(\mu_{i}^{\texttt{len}},\mu_{i}^{\texttt{inc}})=\frac{1}{|S_{i}|}\sum_{(\texttt{len},\texttt{inc})\in S_{i}}(\texttt{len},\texttt{inc})$$

are then used to define the maximal cluster variances in the length and increment directions as

	$\displaystyle\mathrm{Var}_{\texttt{len}}$	$\displaystyle=$	$\displaystyle\max_{i=1,\ldots,k}\frac{1}{|S_{i}|}\sum_{(\texttt{len},\texttt{inc})\in S_{i}}\left|\texttt{len}-\mu_{i}^{\texttt{len}}\right|^{2},$
	$\displaystyle\mathrm{Var}_{\texttt{inc}}$	$\displaystyle=$	$\displaystyle\max_{i=1,\ldots,k}\frac{1}{|S_{i}|}\sum_{(\texttt{len},\texttt{inc})\in S_{i}}\left|\texttt{inc}-\mu_{i}^{\texttt{inc}}\right|^{2},$

respectively. Here, $|S_{i}|$ denotes the number of tuples in cluster $S_{i}$. The aim is to seek the smallest number of clusters $k$ such that

$$\max(\texttt{scl}\cdot\mathrm{Var}_{\texttt{len}},\mathrm{Var}_{\texttt{inc}})\leq\texttt{tol}_{s}^{2}$$

(5)

with a tolerance $\texttt{tol}_{s}$. Once the optimal $k$ has been found, each cluster $S_{1},\ldots,S_{k}$ is assigned a symbol $a_{1},\ldots,a_{k}\in\mathbb{A}$, respectively. Finally, each tuple in the sequence (2) is replaced by the symbol of the cluster it belongs to, resulting in the symbolic representation $S=[s_{1},s_{2},\ldots,s_{n}]\in\mathbb{A}^{n}$.

The final step of the ABBA digitization procedure just outlined, namely the determination of the required number of clusters $k$, requires repeated runs of the k-means algorithm. This algorithm partitions the set of $n$ time series pieces into the specified number of $k$ clusters by minimizing the sum-of-squares error (SSE) between each point and its nearest cluster center. Choosing the right number of clusters and performing repeated SSE calculations is found to be a major computational bottleneck of ABBA, in particular when the number of pieces $n$ and the final parameter $k$ turn out to be large (say, with $n$ in the hundreds and $k=10$ or above). While there are other approaches to estimate the “right” number of clusters, such as silhouette coefficients [41], these still require a clustering to be computed first and hence are of high computational complexity.

While there are a large number of alternative clustering methods available (such as STING [48], BIRCH [54], CLIQUE [1], WaveCluster [44], DenClue [23], spectral clustering [50], Gaussian mixture model [12, 20, 7], DBSCAN [19] and HDBSCAN [8], to name just a few of the more commonly used), many of these are inappropriate for time series data and/or require difficult parameter choices. For example, the density-based DBSCAN algorithm is not appropriate for time series pieces because there is no reason to expect that similar pieces $p_{i}$ would aggregate in regions of high density.

We start off with $n$ data points $p_{1},\ldots,p_{n}$ in $\mathbb{R}^{2}$, obtained as in (3) but sorted in one of the following two ways.

• •

  Binned and lexicographically sorted: if $p_{i}(1)=p_{j}(1)$ and $p_{i}(2)\leq p_{j}(2)$ then we require that $i\leq j$. Here, $p_{i}(1)$ and $p_{i}(2)$ refer to the first and second coordinate of $p_{i}$, respectively. This form of lexicographic sorting is meaningful if the first components $p_{i}(1)$ ($i=1,\ldots,n$) take on only a small number of distinct values. This can be achieved by assigning the $p_{i}$’s into a small number of bins depending on $p_{i}(1)$. Within ABBA this condition is naturally satisfied as each $p_{i}(1)$ corresponds to the scaled length of a piece $P_{i}$, which is an integer.

• •

  Sorted by norm: if $\|p_{i}\|\leq\|p_{j}\|$ then necessarily $i\leq j$, where $\|\cdot\|$ corresponds to a norm. In all our numerical tests we assume that this is the $\ell_{p}$ norm $\|x\|_{p}=(|x_{1}|^{p}+|x_{2}|^{p})^{1/p}$ for $p=2$ (Euclidean norm) or $p=1$ (Manhattan norm).

These sorting conditions are essential for the efficiency of our aggregation algorithm in various ways. First, it is clear that two data points $p_{i}$ and $p_{j}$ with a small distance

$$\operatorname*{dist}(p_{i},p_{j}):=\|p_{i}-p_{j}\|$$

(6)

must also have similar norm values since

$$\big{|}\|p_{i}\|-\|p_{j}\|\big{|}\leq\operatorname*{dist}(p_{i},p_{j}).$$

(7)

Hence, a small distance between norm values is a necessary condition for the closeness of data points. While it is not a sufficient condition, norm-sorted points with adjacent indices are more likely to be close in distance than without sorting, facilitating the efficient assignment of nearby points into clusters. Moreover, the imposed sorting can be used to stop searching for data points near a given point $p_{i}$ prematurely without traversing all of the remaining points.

Indeed with lexicographic sorting, if $p_{i}(1)=p_{i+1}(1)$ and $\operatorname*{dist}(p_{i},p_{i+1})<\min_{j}|p_{i}(1)-p_{j}(1)|$ (the minimum bin width around $p_{i}$), then we know that $p_{i+1}$ is closest to $p_{i}$ among all points $p_{i+1},\ldots,p_{n}$. On the other hand, if $p_{k}(1)>p_{i}(1)+\alpha$ for some $k>i$ and $\alpha>0$, then $\operatorname*{dist}(p_{i},p_{j})>\alpha$ for all $j\geq k$. As a consequence, we can stop searching through the bins for points $p_{j}$ near $p_{i}$ as soon as the first component $p_{j}(1)$ is large enough. Also within a bin, if $p_{k}(1)=p_{i}(1)$ but $p_{k}(2)>p_{i}(2)+\alpha$, we know that $\operatorname*{dist}(p_{i},p_{j})>\alpha$ for all the following indices $j\geq k$.

Similarly, when the data is sorted by norm, then from (7) we know that when $\|p_{k}\|>\|p_{i}\|+\alpha$ for some index $k>i$, we have $\operatorname*{dist}(p_{i},p_{j})>\alpha$ for all $j\geq k$. Hence, we can terminate the search for points near $p_{i}$ as soon as $j$ reaches this index $k$.

With these termination conditions for the search in place, we now propose our aggregation approach in Algorithm 1. It assumes the data points $p_{1},\ldots,p_{n}$ to be sorted in one of the two ways described above, and requires a hyperparameter $\alpha>0$ for controlling the group size. Hereforth, we term “group” is used instead of “cluster” in order to clearly distinguish this approach from other clustering methods. We will discuss the role of the parameter $\alpha$ in Section 4.1 and in the numerical experiments, Section 5.5. Starting with the first sorted point, which is considered the starting point of the first group, the algorithm traverses points in the sorted order and assigns them to the current group if the distance to the starting point does not exceed $\alpha$. A visual illustration of Algorithm 1 is shown in Figure 1.

We emphasize that the starting points are fixed while the groups are being formed, hence they do not necessarily correspond to the mean centers of each group (as would be the case, e.g., in k-means). This is a crucial property of this algorithm as the mean center of a group is only known once the group is fully formed. By using this concept of starting points we avoid repeated center calculations, resulting in significantly reduced computational complexity (see the analysis in Section 4.2). As the points are traversed in sorted order, we use the early stopping conditions for the search whenever applicable. Only when all $n$ points have been assigned to a group, the centers are computed as the mean of all data points associated with a group (Step 6 in Algorithm 1).

The fABBA method for computing a symbolic representation of a time series $T=[t_{0},t_{1},\ldots,t_{N}]$ follows the same compression phase as the original ABBA method described in Section 2. That is, given a tolerance tol, we compute the same polygonal chain approximation as in (1), resulting in the same tuple sequence (2), and the same sequence of scaled points $p_{1},\ldots,p_{n}$ as in (3).

The fABBA digitization then proceeds using Algorithm 1 described in Section 3 to assign the points $p_{1},p_{2},\ldots,p_{n}$ to groups $S_{1},\ldots,S_{k}$, each represented by a symbol in $\mathbb{A}=\{a_{1},\ldots,a_{k}\}$, resulting in the symbolic representation

$$S=[s_{1},s_{2},\ldots,s_{n}]\in\mathbb{A}^{n}.$$

Each group $S_{i}$ has an associated group center $\mu_{i}$ corresponding to the mean of all points associated with that group. These mean centers can be used to reconstruct an approximation of the time series $T$ from its symbolic representation.

Throughout this section we assume that the points $p_{1},p_{2},\ldots,p_{n}$ are already scaled and sorted by their Euclidean norm. Recall that in Algorithm 1 we effectively use a starting point $p_{i}$ as a “temporary” group center when computing the distances $\operatorname*{dist}(p_{i},p_{j})$ in Step 3. By construction, all points in a group, say, $S_{i}=\{p_{i_{1}},p_{i_{2}},\ldots,p_{i_{\ell}}\}$, satisfy

$$\sum_{j=1}^{\ell}\operatorname*{dist}(p_{i_{1}},p_{i_{j}})^{2}\leq\ell\alpha^{2}.$$

(8)

So, upon replacing the temporary group center $p_{i_{1}}$ by the true group mean in Step 6 of the algorithm,

$$\mu_{i}=\frac{1}{\ell}\sum_{j=1}^{\ell}p_{i_{j}},$$

and using that the mean minimizes the sum of squares on the left-hand side of (8), we obtain

$$\mathrm{Var}_{i}:=\frac{1}{\ell}\sum_{j=1}^{\ell}\operatorname*{dist}(\mu_{i},p_{i_{j}})^{2}\leq\alpha^{2}.$$

(9)

This is a variance-based cluster condition similar to the one used in the ABBA method; see also (5).

We now analyze how the use of starting points as “temporary” group centers affects the WCSS defined in (4). The WCSS is a measure for the coherence of the groups. The smaller its value, the more tightly packed are points within each group. In order to make the analysis tractable, we assume that both grouping approaches—measuring distance from a starting point versus measuring distance from the true group mean—lead to the same groupings. For clarity we denote by $\texttt{WCSS}_{\mu}$ the within-cluster-sum-of-squares as defined in (4) using the group means ${\mu}_{i}$ as the centers, while $\texttt{WCSS}_{\mathrm{sp}}$ uses the starting points. It turns out that both of these quantities cannot be arbitrarily large. Indeed, let $p_{1},p_{2},\ldots,p_{n}\in\mathbb{R}^{2}$ be aggregated into $S_{1},\ldots,S_{k}$ with associated starting points $sp_{1},\ldots,sp_{k}$. Then

	$\displaystyle\texttt{WCSS}_{\mathrm{sp}}$	$\displaystyle=\sum_{i=1}^{k}\sum_{p\in S_{i}}\|p-sp_{i}\|^{2}=\sum_{i=1}^{k}\sum_{p\in S_{i},p\neq sp_{i}}\|p-sp_{i}\|^{2}$
		$\displaystyle\leq\sum_{i=1}^{k}(|S_{i}|-1)\alpha^{2}=\alpha^{2}\left(\sum_{i=1}^{k}|S_{i}|-k\right)$
		$\displaystyle=\alpha^{2}(n-k).$

Similarly, it can easily be shown that

$$\texttt{WCSS}_{\mu}\leq\alpha^{2}n:$$

Consider each mean $\mu_{i}$ as an additional element of a the group $S_{i}$, which is also the starting point of that group $S_{i}$. Hence, we are effectively grouping $n+k$ points into $k$ groups and so the $\texttt{WCSS}_{\mu}$ bound follows from the bound $\alpha^{2}(n-k)$ for $\texttt{WCSS}_{\mathrm{sp}}$ when $n$ is replaced by $n+k$.

Both of these bounds are only attained with worst-case point distributions. In practice, the measured $\texttt{WCSS}_{\mathrm{sp}}$ and $\texttt{WCSS}_{\mu}$ are usually smaller than these bounds suggest. This has to do with the fact that, in expectation, each term $\sum_{p\in S_{i}}\|p-sp_{i}\|^{2}$ in $\texttt{WCSS}_{\mathrm{sp}}$ will usually be smaller than the worst-case value $(|S_{i}|-1)\alpha^{2}$. Indeed, let us assume that the points $p$ in a group $S_{i}$ follow a uniform random distribution in the disk $\{(x,y):x^{2}+y^{2}\leq\alpha^{2}\}$ with center 0. (The latter assumption is without loss of generality as the distances within a group are shift-independent.) Then the distance of a random point from the origin is a probability variable, say $R$. Since the area of the disk is $A(\alpha)=\pi\alpha^{2}$, the distance to the origin has a cumulative distribution function

$$F_{R}(r)=\mathbb{P}(R\leq r)=\frac{A(r)}{A(\alpha)}=\frac{r^{2}}{\alpha^{2}},\quad 0\leq r\leq\alpha,$$

and a probability density $f_{R}(r)=F_{R}^{\prime}(r)=2r/\alpha^{2}$. Hence, the expected $\texttt{WCSS}_{\mathrm{sp}}$ is

$$\mathbb{E}[\texttt{WCSS}_{\mathrm{sp}}]=(n-k)\cdot\mathbb{E}[R^{2}]=(n-k)\int_{0}^{\alpha}r^{2}\cdot f_{R}(r)\mathrm{d}r=\frac{\alpha^{2}(n-k)}{2},$$

which is a factor of two smaller than the worst-case.

Let us illustrate the above with an experiment on synthetic data using $n=1000$ data points $p_{i}\in\mathbb{R}^{2}$ chosen as 10 separated isotropic Gaussian blobs each with 100 points. We run Algorithm 1 on this data for various values of $\alpha$, measuring the obtained $\texttt{WCSS}_{\mathrm{sp}}$ as a function of the number of groups $k$. For each grouping we then compute the mean centers of each group and compute the corresponding $\texttt{WCSS}_{\mu}$ as well. Both curves are shown in Figure 2, together with the probabilistic estimate just derived. We find that the $\texttt{WCSS}_{\mathrm{sp}}$ matches the probabilistic estimate very well, while $\texttt{WCSS}_{\mu}$ is smaller. Note that $\texttt{WCSS}_{\mu}$ is the actual WCSS achieved by Algorithm 1 after Step 6 has been completed and all starting points are replaced by group means. We also show the WCSS achieved by k-means++ [3] with the cluster parameter $k$ chosen equal to the number of groups. (K-means++ is a derivative of k-means with an improved initialization using randomized seeding.) As expected, k-means++ achieves an even smaller WCSS as its objective is the minimization of exactly that quantity. However, it is striking to see that all three WCSS values are rather close to each other, with the aggregation-based algorithm achieving $\texttt{WCSS}_{\mu}$ values which are just a moderate multiple of k-means++, while not actually performing any WCSS calculations. (Note that the WCSS is a “squared” quantity. So if two WCSS values differ by a factor of, say 2, in the units of distance this would only be a factor of $\sqrt{2}$.) We will demonstrate in Section 5 that the tiny worsening of WCSS incurred with Algorithm 1 is barely noticeable in the reconstruction of compressed time series.

Let us turn our attention to the computational complexity of Algorithm 1, which is dominated by two operations involving the data: (i) the initial sorting of $n$ data points in Step 1 and (ii) the distance calculations in Step 3. For the sorting operation we can assume that an algorithm with worst-case and average runtime complexity $\mathcal{O}(n\log n)$ has been used, such as introsort [38].

The number of total distance calculations can be counted as follows. For each unassigned point $p_{i}$, Step 4 is guaranteed to include at least this point itself in a group, so the outer loop over all unassigned pieces (beginning with Step 2) will end after at most $n$ iterations. (The trivial distance $\operatorname*{dist}(p_{i},p_{i})=0$ does not need to be computed.) In the worst case, with the early termination condition never invoked and no groups of size larger than one being formed, each $p_{i}$ will be compared to all $p_{j}$ with indices $j>i$. In total, $(n-1)+\cdots+1=n(n-1)/2=\mathcal{O}(n^{2})$ distance calculations are required in this case. On the other hand, if the sorting is perfect so that points belonging to the same group have consecutive indices, and if early termination conditions are always invoked exactly before a new group is formed, every point $p_{j}$ in a group will be involved in exactly one nontrivial distance calculation $\operatorname*{dist}(p_{i},p_{i+1})$, leading to linear complexity $\mathcal{O}(n)$.

In Section 5 we will demonstrate in a comparison involving a large number of time series, that for points $p_{i}$ arising from the ABBA compression step the average number of distance calculations per point is indeed very small, in some cases as small as just 2 and rarely above 10, in particular, when Euclidean norm sorting is used. (See Figure 10 which shows an almost bounded number of the average distance calculations as the number of pieces $n$ increases.) In practice this means that Algorithm 1, and hence the digitization phase of fABBA, is dominated by the $\mathcal{O}(n\log n)$ cost for the initial sorting in Step 1. As a consequence, fABBA outperforms the original ABBA implementation based on k-means++ significantly.

One might wonder whether the use of other clustering algorithms might be more appropriate. We argue that in the context of fABBA the WCSS criterion is indeed the correct objective function, advocating distance-based clustering methods like k-means and Algorithm 1. This is justified by the observation that the error of the ABBA reconstruction can be modelled as a Brownian bridge (see the analysis in [15, Section 5], which we do not repeat here) and the overall error is controlled by the total cluster variance, which coincides with the WCSS up to a constant factor.

Nevertheless, we provide here for completeness a comparison of the digitization performance obtained with other clustering algorithms, namely k-means++, spectral clustering, Gaussian mixture model (GMM), DBSCAN, and HDBSCAN. To this end we sample 100 time series of length 5000 from a standard normal distribution and then measure the mean squared error (MSE) and dynamic time warping (DTW) distance [43] of the reconstruction, digitization runtime, and the number of distinct symbols used. For the computation of DTW distances we use the implementation in the tslearn library [46] with all default options. For HDBSCAN we use the open-source software [37, 36]. The remaining algorithms are implemented in scikit-learn [39]. All algorithms are executed through the Cython compiler for a fair runtime comparison.

The number of distinct symbols $k$ for the symbolic representation is not a user-specified parameter for Algorithm 1, DBSCAN and HDBSCAN, but we must specify $k$ for the other clustering algorithms (k-means++, spectral clustering, GMM). Therefore, we first run fABBA with $\text{{tol}}=0.15$ (digitization) and $\alpha=0.5$ (aggregation) and use the number of generated groups as the cluster parameter for k-means++, spectral clustering, and GMM (weights are initialized using k-means). For DBSCAN, we use all default options for the parameters except for setting $\epsilon$ identical to $\text{{tol}}=0.15$ and specifying $\texttt{min\_samples}=2$. For HDBSCAN we use $\texttt{min\_samples}=2$, $\texttt{min\_cluster\_size}=12$ and $\texttt{alpha}=0.5$ to achieve roughly the same number of symbols $k$ as used by the other algorithms.

The results are given in Figure 3 and Table LABEL:table:com_other_clusterings. In terms of average runtime, our new aggregation Algorithm 1 is the fasted (51.215 ms), closely followed only by DBSCAN (59.057 ms). The box-and-whiskers plots in Figure 3 depict the distributions of the four metrics, with each box showing the first quartile $Q_{1}$, median, and third quartile $Q_{3}$, while the whiskers extend to points that lie within the 1.5 interquartile range of $Q_{1}$ and $Q_{3}$. Observations that fall outside this range are displayed separately. We see that k-means++ is the only method that outperforms Algorithm 1 in terms of MSE and DTW error, but only by a small margin and at the cost of being almost 8 times slower. Therefore, Algorithm 1 indeed seems to be the most appropriate method to use in this test.

To further the findings in Section 4.3 we perform a direct comparison of Algorithm 1 with k-means++ [3] on a different dataset. We use the Shape sets, a dataset containing eight benchmarks corresponding to different shapes that are generally difficult to cluster. The benchmarks are referred to as (1) Aggregation, (2) Compound, (3) Pathbased, (4) Spiral, (5) D31, (6) R15, (7) Jain, and (8) Flame.

The test is conducted by first normalizing each dataset and then running Algorithm 1 using the tolerance $\alpha=1.5$. Afterwards k-means++ is run using the same number of clusters as determined by Algorithm 1. The number of data points in each benchmark and the number of (predicted) labels are listed in Table LABEL:table:clusteringbenchmark. Note that this comparison is not primarily about the accuracy of identifying clusters “correctly”. The grouping in Algorithm 1 is controlled by the parameter $\alpha$ and the variance-based criterion (9), not by a criterion for cluster quality. Nevertheless, when comparing the adjusted Rand indices [24] of the obtained clusters in Figure 4 (a), Algorithm 1 is surprisingly competitive with k-means++: on four of the eight benchmarks Algorithm 1 outperforms k-means++ and on two other benchmarks the adjusted Rand indices are almost identical. The main take-away from this test, however, is that the runtime of Algorithm 1 with the same number of clusters is significantly reduced, by factors well above 50, in comparison to k-means++. See Figure 4 (b).

Performance profiles are an effective and widely used tool for evaluating the performance of a set of algorithms on a large set of test problems [14]. They can be used to compare quantities like runtime, number of function evaluations, number of iterations, or the memory used.

Let $\mathcal{P}$ represent a set of problems and let $\mathcal{S}$ be a set of algorithms. Let $s_{i,j}\geq 0$ be a performance measurement of algorithm $i\in\mathcal{S}$ run on problem $j\in\mathcal{P}$, assuming that a smaller value of $s_{i,j}$ is better. For any problem $j\in\mathcal{P}$, let $\widehat{s}_{j}=\min\{s_{i,j}:i\in\mathcal{S}\}$ be the best performance of an algorithm on problem $j$ and define the performance ratio $r_{i,j}$ as

$$\small r_{i,j}=\frac{s_{i,j}}{\widehat{s}_{j}}.$$

(10)

We set $r_{i,j}=0$ if solver $i$ fails on problem $j$. We define

$$\small k(r_{i,j},\theta)=\left\{\begin{array}[]{rcl}1&&\text{if}\hskip 5.0pt{r_{i,j}<\theta}\\ 0&&\text{otherwise}\\ \end{array}\right.$$

(11)

where $\theta\geq 1$. The performance profile of solver $i$ is

$$\rho_{i}(\theta)=\frac{\sum_{j\in\mathcal{P}}k(r_{i,j},\theta)}{|\mathcal{P}|},\quad\theta\geq 1.$$

(12)

In other words, $\rho_{i}(\theta)$ is the empirical probability that the performance ratio $r_{i,j}$ for solver $i\in\mathcal{S}$ and each problem $j\in\mathcal{P}$ is within a factor $\theta$ of the best possible ratio. The function $\rho_{i}$ is the cumulative distribution function of the performance ratio. Its graph, the performance profile, reveals several performance characteristics in a convenient way. In particular, solvers $i$ with fast increasing performance profile are preferable, provided the problem set $\mathcal{P}$ is suitably large and representative of problems that are likely to occur in applications.

We perform experiments running the SAX, 1d-SAX, ABBA, and fABBA methods on the 201,161 time series included in the UCR 2018 Archive [11]. All of the four methods fall in the class of numerical transforms, are based on piecewise polynomial approximation, and hence are rather easily comparable. Our selection of symbolic representations is guided by the public availability of Python implementations, e.g., within the tslearn library [46]. We follow the experimental design in [15, Section 6] to allow for a direct comparison. Each original time series $T=[t_{0},t_{1},\ldots,t_{N}]$ is first processed by the fABBA compression with the initial tolerance ${\texttt{tol}}=0.05$, returning the number of pieces $n$ required to approximate $T$ to that accuracy. The compression rate $\tau_{c}$ is defined as

$$\small\tau_{c}=\frac{n}{N}\in(0,1].$$

(13)

If the compression rate $\tau_{c}$ turns out to be larger than 0.2, we increase the tolerance tol in steps of 0.05 and rerun the fABBA compression until $\tau_{c}\leq 0.2$. If this cannot be achieved even with a tolerance of ${\texttt{tol}}=0.5$, we consider the time series as too noisy for symbolic compression and exclude it from the test. (In other words, if the compression is almost as long as the original time series, there is no point using any symbolic representation.) The numbers of time series finally included in the test are listed in Table 4.

After the compression phase we run the fABBA digitization using 2-norm sorting with the scaling parameter ${\texttt{scl}}=1$, corresponding to unweighted 2d clustering. The digitization rate $\tau_{d}$ is

$$\tau_{d}=\frac{k}{n}\in(0,1]$$

(14)

where $k$ is the number of symbols assigned by fABBA. This is then followed by running SAX, 1d-SAX, and ABBA (${\texttt{scl}}=1$) with the same number of pieces $n$ and symbols $k$. Finally, we reconstruct the time series from their symbolic representations and measure the distance to the original time series $T$ in the Euclidean norm and DTW distance. We also measure the Euclidean and DTW distance of the differenced time series (effectively ignoring errors in representing linear trends).

The performance profiles for the parameter choices $\alpha=0.1$ and $\alpha=0.5$ are shown in Figure 5 and Figure 6, respectively. The mean values of the compression rate $\tau_{d}$ for $\alpha=0.1$ and 0.5 are approximately 0.896 and 0.615, respectively. The performance profiles reveal that ABBA and fABBA consistently achieve similar reconstruction accuracy in all four distance metrics and lead to significantly more accurate reconstructions than SAX and 1d-SAX. The only exception is $\alpha=0.5$ with the Euclidean distance metric, in which case all four algorithms have similar reconstruction performance; see Figure 6(a).

For illustration, some examples of reconstructions for $\alpha=0.1$ and $\alpha=0.5$ are shown in Figure 7 and 8, respectively. The fABBA reconstruction closely follows that of ABBA and both representations capture the essential shape features of the time series.

In Figure 9 we compare the digitization runtime of fABBA and ABBA with 2-norm sorting and scaling parameter $\text{{scl}}=1$ for all problems considered in the previous section. This is done using a histogram of the ratio

$$\frac{\text{{ABBA runtime}}}{\text{{fABBA runtime}}}\vskip 2.84526pt$$

for each time series. We find that fABBA is always faster than ABBA and the speedup factor can be as large as $\approx 50$ with a mode of $\approx 30$.

In order to explain this speedup from the viewpoint of Algorithm 1, we visualize in Figure 10 the average number of distance calculations depending on the number of pieces $n$ and when different sorting methods are used, namely lexicographic sorting and 1- and 2-norm sorting. We find that in particular for norm sorting, the number of average distance calculations can be as small as $\approx 2$ and it hardly ever exceeds $\approx 10$. As a consequence, the computational cost of the aggregation is roughly linear in the number of pieces $n$. With the number of pieces scaling linearly with the length of the original time series $N$, the aggregation phase is of linear complexity in $N$, and the overall complexity of Algorithm 1 is dominated by $\mathcal{O}(n\log n)$ operations for the initial sorting.

We now perform a more exhaustive parameter study for fABBA in order to gain insights into its robustness with respect to the parameters. We consider three different initial sortings for Step 1 of Algorithm 1 and vary the $\alpha$ parameter from $0.1,\ldots,0.9$. Table 5 displays the resulting average digitization rate $\tau_{d}$, Euclidean (2-norm) and DTW distance between the reconstruction and the original time series, the total runtime for the digitization (in milliseconds), the number of distance calculations (denoted by dist), and the number of symbols used. We observe that when $\alpha$ is increased, fewer symbols are used and the digitization is performed faster. On the other hand, and as expected, the reconstruction errors measured in the Euclidean and DTW distance become larger when $\alpha$ is increased. Generally, there does not seem to be a significant difference in reconstruction accuracy for different norm sortings, but in view of distance calculations and runtime we recommend using 2-norm sorting as the default.

To demonstrate that fABBA can easily be applied to data not directly arising as time series, we discuss an application to image compression. This example also allows us to provide a different visual insight into fABBA’s reconstruction performance. Our approach is very simple: we select 12 images from the Stanford Dogs Dataset [13] and resize each image to $250\times 250$ pixels. We then reshape each $250\times 250\times 3$ RGB array into a univariate series (unfolded so that the values of the R channel come first, followed by the G and B channels) and then apply fABBA to compress this series. For the visualization of the reconstruction we simply reshape the reconstructed series into an RGB array of the original size. The results are presented in Figure 11 and Table LABEL:table:img_com, testing different compression tolerances tol while keeping the digitization tolerance $\alpha=0.001$ fixed. We find that, even for the rather crude tolerance of $\text{{tol}}=0.5$, the essential features in each image are still recognizable despite the rather high compression rate (average $\overline{\tau}_{c}=0.184$) and digitization rate (average $\overline{\tau}_{d}=0.227$). This compressed lower-dimensional representation of images might be useful, e.g., for classification tasks, but this will be left for future work.