Joint symbolic aggregate approximation of time series

The ABBA compression step aims to compute an adaptive piecewise linear continuous approximation of $T$, that is, to obtain time series pieces $P=[(\text{{len}}_{1},\text{{inc}}_{1}),\ldots,(\text{{len}}_{N},\text{{inc}}_{N})]\in\mathbb{R}^{N\times 2}$, followed by a reasonable digitization that results in symbolic sequence $A=[a_{1},a_{2},\ldots,a_{N}]\in\mathcal{L}^{N}$, $N\ll n$, and each $a_{j}$ is an element of a finite alphabet set $\mathcal{L}$ where $|\mathcal{L}|\ll N$. $\mathcal{L}$ can be referred to as dictionary in the procedure. The ABBA compression adaptively selects $N+1$ indices $i_{0}=0<i_{1}<\cdots<i_{N}=n$ given a tolerance tol so that the time series $T$ is well approximated by a polygonal chain going through the points $(i_{j},t_{i_{j}})$ for $j=0,1,\ldots,N$. This results in a partition of $T$ into $N$ pieces $p_{j}=(\text{{len}}_{j},\text{{inc}}_{j})$ that is determined by $T_{i_{j-1}:i_{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. Visually, each piece $p_{j}$ is represented by a straight line connecting the endpoint values $t_{i_{j-1}}$ and $t_{i_{j}}$ This partitioning criterion is the squared Euclidean distance of the values in $p_{j}$ from the straight polygonal line is upper bounded by $(\texttt{len}_{j}-1)\cdot\texttt{tol}^{2}$. For simplicity, given an index $i_{j-1}$ and starts with $i_{0}=0$, the procedure seeks the largest possible $i_{j}$ such that $i_{j-1}<i_{j}\leq n$ and

Each linear piece $p_{j}$ of the resulting polygonal chain $\widetilde{T}$ is referred to as a tuple $(\texttt{len}_{j},\texttt{inc}_{j})$, where $\texttt{inc}_{j}=t_{i_{j}}-t_{i_{j-1}}$ is the increment in value, i.e., the subtraction of ending and starting value of $T_{i_{j-1}:i_{j}}$. The whole polygonal chain can be recovered exactly from the first value $t_{0}$ and the tuple sequence $p_{1},p_{2},\ldots,p_{N}$, i.e.,

where the reconstruction error of this representation is with pinned start and end points, and can be naturally modeled as a Brownian bridge.

The next step is referred to as digitization, which we further transformed the resulting polygonal chain $\widetilde{T}$ into the symbolic representation in the form of (1).

Following (Elsworth and Güttel, 2020a), prior to digitizing, the tuple lengths and increments are separately normalized by their standard deviations $\sigma_{\texttt{len}}$ and $\sigma_{\texttt{inc}}$, respectively. After that, further scaling is employed by using a parameter scl to assign different weights to the length of each piece $p_{i}$, which denotes importance assigned to its length value in relation to its increment value. Hence, the clustering is effectively performed on the scaled tuples

In particular, if $\texttt{scl}=0$, then clustering will be only performed on the increment values of $P$, while if $\texttt{scl}=1$, the lengths and increments are clustered with equal importance.

The steps after normalization proceed with a lossy compression technique, e.g., vector quantization (VQ), which is often achieved by mean-based clustering. The concept of vector quantization can be referenced in (Gray, 1984; Dasgupta and Freund, 2009). Given an input of $N$ vectors $P=[p_{1},\ldots,p_{N}]\in\mathrm{R}^{\ell\times N}$, VQ seeks a codebook of $k$ vectors, i.e., $C=[c_{1},\ldots,c_{k}]\in\mathrm{R}^{\ell\times k}$ such that $k$ is much smaller than $N$ where each $c_{i}$ is associated with a unique cluster $S_{i}$. A quality codebook enables the sum of squared errors SSE to be small enough to an optimal level. Suppose $k$ clusters $S_{1},S_{2},\ldots,S_{k}\subseteq P$ are computed, VQ aims to minimize

where $\phi$ denotes energy function, $c_{i}$ denotes the center of cluster $S_{i}$ and $\mathrm{dist}(p_{i},p_{j})$ often denotes the Euclidean norm $\|p_{i}-p_{j}\|_{2}$. We often choose the mean center $\mu_{i}$ as $c_{i}$ for Euclidean space, i.e., $\mu_{i}=\frac{1}{|S_{i}|}\sum_{p\in S_{i}}p$, and then (5) can be written as $\texttt{SSE}=\sum_{i=1}^{k}|S_{i}|\text{Var}S_{i}$. Lyold’s algorithm (Lloyd, 1982a) (also known as k-means algorithm) is a suboptimal solution of vector quantization to minimize SSE.

The ABBA digitization can be performed by a suitable partitional clustering algorithm that finds $k$ clusters from $P\in\mathbb{R}^{2\times N}$ such that the sum of Euclidean distance SSE constructed by $C$ is minimized. The obtained codebook vectors are referred to symbolic centers here. Each symbolic center is associated with an identical symbol and each time series snippet $p_{i}$ is assigned with the closest symbolic center $c^{i}$ associated with its symbol

The symbolic centers to symbols are one-to-one mapping, denoted by $I_{d}:C\to A$ , thus the digitization $f_{d}:P\to A$ is given by

Each symbol is associated with a unique cluster. In practice, each clustering label (membership) corresponds to a unique byte-size integer value. The symbols used in ABBA can be represented by text characters, which are not limited to English alphabet letters—often more clusters will be used. Each character in most computer systems is used by the ASCII strings with a unique byte-size integer value (a unique cluster membership). Besides, it can be any combination of symbols, or ASCII representation.

Besides, it is fun to discuss compression rates in some cases. The digitization is the key to compression rate, which is the size of codebook $C$ (i.e., the number of distinct symbols $|\mathcal{L}|$) divided by the length of time series. We use $\tau_{c}$ to denote the compression rate, which is given by

The inverse symbolization refers to the process from $A$ to $\widehat{T}$, the intuition is to reconstruct time series from (1) such that the reconstructed time series $\widehat{T}$ is as close to $T$ as possible. The inverse symbolization contains three steps.

The first step is referred to as inverse-digitization, simply written as $f_{d}^{-1}$, which uses the $k$ representative elements $c_{i}$ (in terms of, e.g., mean centers or median center of the groups $S$) from codebook $C$ to replace the symbol in $A$ orderly, and thus results in a 2-by-$N$ array $\widetilde{P}$, i.e., an approximation of $P$, where each $\widetilde{p}_{i}\in\widetilde{P}$ is the closest symbolic center $c\in C$ to $p_{i}\in P$. The inverse digitization often leads to a non-integer value to the reconstructed length len, so (Elsworth and Güttel, 2020a) proposes a novel rounding method, which is referred to as quantization, to align the cumulated lengths with the closest integers. The method is as follows: start with rounding the first length into an integer value, i.e., $\widehat{\text{{len}}}_{1}:=\text{round}(\widetilde{\text{{len}}}_{1})$ and calculate the rounding error $e:=\widetilde{len}_{1}-\widehat{len}_{1}$. The the error is added to the rounding to $\widetilde{\text{{len}}}_{2}$, i.e., $\widehat{\text{{len}}}_{2}:=\text{round}(\widetilde{\text{{len}}}_{2}+e)$ and new error $e^{\prime}$ is calculated as $\widehat{\text{{len}}}_{2}+e-\widetilde{\text{{len}}}_{2}$. Then $e^{\prime}$ is involved in the next rounding similarly. After all rounding is computed, we obtain

where increments inc are unchanged, i.e., $\widehat{\text{{inc}}}=\widetilde{\text{{inc}}}$. Then, the whole polygonal chain can be recovered exactly from the initial time value $t_{0}$ and the tuple sequence (9) via the inverse-compression.

The lower reconstruction error means a higher approximation accuracy. The reconstruction error can be defined by mean squared error (MSE), which is given by

The k-means problems aim to find $k$ clusters within data in $d$-dimensional space, so as to minimize the (5). However, solving this problem is NP-hard even $k$ is restricted to $2$ (Drineas et al., 2004; Dasgupta and Freund, 2008) or in the plane (Mahajan et al., 2012). Typically, the sub-optimal k-means problem can be solved by Lloyd’s algorithm (Lloyd, 1982b). In the implementation of ABBA software³³3https://github.com/nla-group/ABBA, the k-means algorithm is performed by scikit-learn library (Pedregosa et al., 2011) which runs a few times (this is controlled by the parameter n_init⁴⁴4The default to scikit-learn is 10 before the version 1.3.2.) of Lylod’s algorithm with $D^{2}$ seeding and pick up the best result. In the setting of this paper, we found setting n_init to 1 is good enough for the ABBA performance.

As already mentioned, Lloyd’s algorithm is a widely used method to solve the k-means problem, it starts with uniformly sampling $k$ centers from data, often referred to as seeding, and then each point is allocated to a cluster with the closest center, and the mean centers of clusters are recomputed again. The procedure keeps repeating until the iteration converges.

The seeding has a great impact on the final result. The improved algorithm is combined with optimal seeding “$D^{2}$ weighting” introduced by (Arthur and Vassilvitskii, 2007), which can significantly improve Lloyd’s algorithm. Lloyd’s algorithm with $D^{2}$ weighting is called the “k-means++” algorithm. The k-means algorithm with “$D^{2}$ weighting” shows $O(\log k)$-competitive with the optimal clustering.

ABBA digitization using this clustering method has been shown incredibly slow speed, though the reconstruction error meets the needs of most applications. It is very desired to design a faster clustering alternative while retaining the original reconstruction error to an ultimate degree. For this reason, we propose a sampling-based k-means clustering algorithm that can address the above concern. The idea is to perform k-means++ on a uniform sample of data where only $r$ percent of original data is used. The algorithm is as described in Algorithm 3. Section 6 will demonstrate its performance empirically.

Figure 3 shows the simulation of k-means++ and sampling-based k-means (sampling $r=10\%$ of data) on Gaussian blobs data with 10 clusters, proceeding by increasing data sizes. The result is as illustrated in Figure 3. We can observe that sampling-based k-means runs in a fraction of the time compared to k-means++, while giving competitive performance in terms of adjusted mutual information (AMI).

The greedy aggregation is introduced in (Chen and Güttel, 2022a), which proceeds first by sorting the data and performing greedy aggregation of data into groups. The sorting order naturally avoids unnecessary computations in aggregation by triggering an early stopping. The codebook set is constructed by the mean centers of the groups resulting from the aggregation as a suboptimal solution to the k-means problem. Though its accuracy is less significant than Lylod’s algorithm, it achieves a significant speedup and the SSE is upper bounded by $\alpha^{2}(N-k)$ for $N$ data points.

Sorting is essential to the success of aggregation in our context. Since sorting can determine the starting points selection and forming of groups, even helps to discard unnecessary distance computations. A bad sorting will result in inefficiency of aggregation and bad-performed SSE. For example, (Chen and Güttel, 2022b) proposes PCA sorting which ensures the pairwise distance between $p_{i}$ and $p_{j}$ is bounded by $|\varrho_{i}-\varrho_{j}|+2\sigma_{2}^{2}$ where $\sigma_{2}$ is the second largest singular value of data matrix $P$.

As aforementioned, digitization aims to partition $P$ described in (3) into $k$ clusters $S_{1},\ldots,S_{k}$ such that (5) is minimized. The tolerance-oriented digitization enables the natural relationship between compression tol and digitization $\alpha$. In this section, we discuss a novel way to eliminate the need of choosing a parameter for fABBA digitization. The Lemma 4.2 shows the reconstruction error still ensure the pin of start and end of time series.

Proof to Lemma 4.2 is as follows:

Also, we know that $t_{i_{N}}=t_{0}+\sum_{i}^{N}\texttt{inc}_{i}$, hence, the reconstruction $\widetilde{T}$ starts and ends at the same values as $T$ so is $\widehat{T}$.

We assume variance of length and increment of pieces, denoted by $\texttt{Var}_{\texttt{len}}$ and $\texttt{Var}_{\texttt{inc}}$, are:

Here we suppose the aggregation is performed on the length and increment values (1-dimensional data) of pieces simultaneously, which is referred to as hierarchical aggregation, and we denote the digitization tolerance for length and increment $\alpha_{\texttt{len}}$ and $\alpha_{\texttt{inc}}$, respectively. Obviously, we have

We assume $\alpha_{\texttt{len}}=\alpha_{\texttt{inc}}=\alpha$, this yields

In the following, we will demonstrate that the Brownian bridge property as illustrated in (Elsworth and Güttel, 2020a) still holds in hierarchical aggregation for time series reconstruction. Though the length of each piece may not be consistent because of rounding error, we assume the length of the reconstructed time series is equal to the original length, i.e., $\sum_{i=0}^{N}\widehat{\texttt{len}}_{i}=\sum_{i=0}^{N}\texttt{len}_{i}$ (In practice, the assumption is true in most cases, but in some special cases, this does not hold true because of rounding). To simplify the modeling and facilitate the analysis, we consider only aggregating increment and assume each cluster of increment has the same mean length, i,e, $\mu_{i}^{\texttt{len}}=n/N$. The local deviation of the increment and length value of $\widehat{T}$ on piece $P_{\ell}$ from the true increment and length of $T$, which are given by

respectively.

The global incremental errors, i.e., the accumulated incremental errors, according to (15), $\ddot{e}_{i_{j},\texttt{inc}}$ are given by:

Also, we must consider the error arisen from the rounding error of length. Similarly, the global length errors, i.e., the accumulated length errors, according to (12) and (13), can be calculated as:

The global error of the reconstructed time series, denoted by $e_{i_{j}}$, is caused by errors from reconstructed length and increment. Up to this point, the global error of the reconstructed time series is still difficult to determine since the estimated error caused by the length displacement is hard to get, so we consider an approximation:

According to (14) and Lemma 4.2, the $\dot{e}_{\ell,\texttt{inc}}$ is bounded by $\alpha$, and $\mathbb{E}(\ddot{e}_{i_{j},\texttt{inc}})=0$ since they are consistent with the deviations from their respective cluster center. Also, $\ddot{e}_{i_{0},\texttt{inc}}=\ddot{e}_{i_{n},\texttt{inc}}=0$ as proved earlier. Therefore, referred to (Elsworth and Güttel, 2020a), we can model a random process of incremental errors $e_{i_{j}}$, and its associated variance:

Following (Elsworth and Güttel, 2020a), $e_{i_{j}}$ is considered to stay $\eta$ standard deviations away from its zeros mean. That is, we consider a realization

Then, with (17) and (19), we have,

The process of $e_{i_{j}}$ is modeled as a Brownian bridge following (Elsworth and Güttel, 2020a). Considering the interpolated quadratic function on the right-hand side is concave, based on the linear stitching procedure used in the reconstruction and by piecewise linear interpolation of the incremental errors from the course time grid $i_{0},i_{1},\ldots,i_{N}$ to the fine time grid $i=0,1,\ldots,n$, it is natural to deduce that

Therefore, the squared Euclidean norm of this fine-grid “worst-case” realization is upper bounded by

It has previously been established that $\texttt{euclid}(T,\widetilde{T})^{2}\leq(n-N)\cdot\texttt{tol}^{2}$ by (Elsworth and Güttel, 2020a), and based on this (4.3) is the worst-case realization of the Brownian bridge and thereby we have a probabilistic bound on the error incurred from digitization. By making $\texttt{euclid}(T,\widetilde{T})^{2}=\sum_{i=0}^{n}e_{i}^{2}$, we can smartly choose

For simplicity we can set this hyperparameter controlling the tolerance of length to be the same as that of increment, i.e., $\alpha_{\texttt{len}}=\alpha_{\texttt{inc}}=\alpha$. Therefore, the parameter of digitization is automatically determined by the compression tolerance, resulting in a non-parametric and error-bounded digitization procedure.

The procedure detailed above is referred to as auto digitization. On top of that, the method introduced in the Section 5 of (Elsworth and Güttel, 2020a) can also be used to approximate an error-bounded fABBA digitization and eliminate the need for tuning $\alpha$, however, this is not practical as it results in a linear relationship between tol and $\alpha$—the difference is as shown in Figure 4 which shows the method in (Elsworth and Güttel, 2020a) depicts a straight line (marked as green color). As a consequence, we can see our method (marked as orange color) as an improvement for choosing $\alpha$ to some degree.

The UEA Archive contains 30 multivariate time series datasets with a variety of dimensions and lengths. The datasets are very huge, therefore it is inefficient for the original ABBA and its variant fABBA to perform computations one at a time, so we select some datasets in the UEA Archive for the test, which are summarized in Table 2. Though impossible for ABBA and fABBA to symbolize time series in each dimension of multivariate time series with unified symbols, we still use them for benchmarking but without considering the symbolic consistency. In order to unify their compressed time series pieces $P$ as specified in (3), we use partitional compression for all four competing methods and reassign the subset of output corresponding to each multivariate time series dimension to ABBA and fABBA. We start with performing the partitional compression as described in Algorithm 5 with tol of 0.01, then for the two JABBA variants we perform their digitization all at once while for ABBA and fABBA we just perform their digitization on time series pieces for each dimension of the multivariate time series one at a time, and then record the runtime for their digitization, respectively.

It’s known that the more symbols are used the more accurate the reconstruction of the representation is. In order to use roughly the same number of symbols for each method as much as possible, we first perform JABBA (GA) digitization using (21) to confirm an $\alpha$ value, and a total number of symbols used for the multivariate time series, denoted by $k_{m}$, then we feed the number of symbols $k_{m}$ for JABBA (VQ) digitization (see Algorithm 3, we set $r$ to 0.5, same in the following). Then we feed the same $\alpha$ value to fABBA digitization and use $k_{m}/d$ for ABBA where $d$ is the dimension of the multivariate time series. As a consequence, the number of symbols used will be unified for ABBA, JABBA (VQ), and JABBA (GA), but not guaranteed for fABBA since its digitization is tolerance-oriented.

Table 3 showcases the average value of MSE, dynamic time warping (DTW), runtime for digitization, and the number of symbols used for each dataset. Information of compression tolerance tol used for each dataset is also given in Table 3. Accordingly, JABBA (GA) shows significant speedup over fABBA by order of magnitude though both use the same GA-based digitization. The speedup of JABBA (VQ) over ABBA is also remarkable while the reconstruction error of JABBA (VQ) is lower than ABBA in five out of the eight datasets though using sampling-based k-means is employed. Additionally, an example of reconstruction from symbolic representation for multivariate time series is presented in Figure 6, which only shows 4 out of 61 dimensions.

In this experiment, we will compare ABBA, fABBA, JABBA (GA), and JABBA (VQ) on synthetic Gaussian noise series in terms of runtime, and reconstruction accuracy with various number of time series partitions. The reconstruction accuracy is measured by MSE here.

We used Gaussian noises as the time series for benchmarking. The data generated for the test are of length 100,000 with zero mean and unit standard deviation. We first ran fABBA with $\text{{tol}}=0.01$ and $\alpha=0.05$ to compute the number of symbols it used. This simulation used 358 symbols accordingly. Second, we ran ABBA by feeding the same number of symbols fABBA used to $k$, i.e., $358$ symbols. After that, we run the JABBA (VQ) and JABBA (GA) with varying partitions by the same tol and specifying a consistent hyperparameter setting for digitization, i.e., $k=358$ and $\alpha=0.05$, respectively. The number of threads scheduled for JABBA is set the same as the partition number. The result shows the compression rate $\tau_{c}$ computed for ABBA, fABBA, JABBA (GA), and JABBA (VQ) are , respectively.

The experimental result is as exhibited in Figure 7 and Figure 8. We mainly compare the methods with the same digitization technique. We can see that there is an obvious negative correlation between reconstruction error and the number of partitions, this can be explained by the increasing partition points that will be used for reconstruction. Figure 7 also shows that JABBA (VQ) which uses sampling-based k-means achieves similar performance against ABBA regarding MSE while performing speedup by orders of magnitude, a similar result applies to JABBA (GA) and fABBA.

The parallel speedup $m$ processors is given by

where $\upsilon(m)$ is referred to as the runtime of the $m$ processors. Without loss of generality, we only evaluate the speedup of Parallelism for JABBA (GA) as shown in Figure 9. We can see the speedup $\Phi(m)$ scale almost linearly with the number of threads $m$. Since our algorithm is partially parallel in compression, which is hindered by the sequential part of the algorithm, that is, the digitization. This phenomenon can be naturally explained by Amdahl’s law which gives the theoretical speedup at a fixed workload where there are limits on the benefits one can derive from parallelizing a computation.