An average-compress algorithm for the sample mean problem under dynamic time warping

Computing a sample mean of time series under dynamic time warping is NP-hard. Consequently, there is an ongoing research effort to devise efficient heuristics. The majority of heuristics have been developed for the constrained sample mean problem that assumes a solution of predefined length. In contrast, research on the unconstrained sample mean problem is underdeveloped. In this article, we propose a generic average-compress (AC) algorithm to address the unconstrained problem. The algorithm alternates between averaging (A-step) and compression (C-step). The A-step takes an initial guess as input and returns an approximation of a sample mean. Then the C-step reduces the length of the approximate solution. The compressed approximation serves as initial guess of the A-step in the next iteration. The purpose of the C-step is to direct the algorithm to more promising solutions of shorter length. The proposed algorithm is generic in the sense that any averaging and any compression method can be used. Experimental results show that the AC algorithm substantially outperforms current state-of-the-art algorithms for time series averaging.

Mean time series (blue) of the two sample time series x (1) and x (2) shown in red. Both time series have a single peak but are out of phase and slightly vary in speed. We may think of x (1) and x (2) as the daily average temperature of some region during the summer at two different years. Based on this information, a typical summer of this region has a single extreme heat wave. The arithmetic mean μ = (x (1) + x (2) )/2 in Fig. 1a has two attenuated peaks suggesting that a typical summer has two moderate heat waves. In contrast, the DTW mean z in Fig. 1b captures the characteristic properties of x (1) and x (2) and shows a single peak as a representative summary of both sample peaks. (Color figure online) temporal dynamics, that is in length, speed, and shifts in phase. For example, the same word can be uttered with different speaking speeds. Similarly, monthly temperature or precipitation extremes of certain regions can differ in duration and may occur out of phase for a period of a few weeks.
To account for temporal variations in proximity-based time series mining, the dynamic time warping (DTW) distance is often the preferred choice of proximity measure [1,3,4]. An intricate problem in DTW-based time series mining is time series averaging. The problem consists in finding a typical representative that summarizes a sample of time series. Different forms of time series averaging have been applied to improve nearest neighbor classifiers [20,31,32], to accelerate similarity search [37], and to formulate k-means clustering in DTW spaces [17,30,32,36]. Figure 1 presents an example illustrating why the arithmetic mean can be inappropriate for time series averaging and motivates a concept of mean under dynamic time warping that can cope with temporal variations.
Time series averaging itself and as a subroutine of data mining tasks is inspired by the fundamental concept of mean in statistical inference. One central path in statistical inference departs from the mean, then leads via the normal distribution and the Central Limit Theorem to statistical estimation using the maximum likelihood method. The maximum likelihood method in turn is a fundamental approach that provides probabilistic interpretations to many pattern recognition methods.
This central path is well-defined in Euclidean spaces, but becomes obscure in mathematically less structured distance spaces. Since an increasing amount of non-Euclidean data is being collected and analyzed in ways that have not been realized before, statistics is undergoing an evolution [26]. Examples of this evolution are contributions to statistical analysis of shapes [5,13,24], complex objects [29], tree-structured data [14,40], graphs [15,19,23], and time series [7,17,32].
Though the volume of time series data currently collected exceeds those of the other data structures mentioned above, the concept of a mean in DTW spaces is least understood. However, a better understanding of time series averaging is the first step towards devising sound pattern recognition methods based on time series averaging such as k-means clustering.
Examples of how the lack of a clear understanding of time series averaging may lead the field astray can be found in [7,21].
As for other non-Euclidean distance spaces, the standard approach to time series averaging in DTW spaces is based on an idea by Fréchet [16]: Suppose that S = {x 1 , . . . , x N } is a sample of N time series. A sample mean of S is any time series μ that globally minimizes the Fréchet function where dtw(x, y) is the DTW-distance and U is a set of time series of finite length. The search space U typically takes two forms: 1. Unconstrained form: U is the set of all time series of finite length. 2. Constrained form: U is the set of all time series of length n.
A sample mean is guaranteed to exist in either case but may not be unique [22]. In addition, computing a sample mean is NP-hard [8].
In general, data mining algorithms such as k-means, vector quantization, and anomaly detection that rely on a sample mean typically expect a true sample mean in their objective function, that is a sample mean that minimizes the variance in a Euclidean sense or more generally the Fréchet function in any distance space. As experiments indicate [7], better sample means result in better solutions of the k-means algorithm. The k-means algorithm needs to recompute the sample mean several times. The error of suboptimal solutions propagates during evolution of the k-means algorithm resulting in possibly distorted clusters. Similarly, a suboptimal sample mean in anomaly detection could result in lower precision and recall.
Consequently, there is an ongoing research on devising heuristics for minimizing the Fréchet function. Most contributions focus on devising and applying heuristics for the constrained sample mean problem. State-of-the-art algorithms are stochastic subgradient methods [35], majorize-minimize algorithms [17,30,35], and soft-DTW [11]. In contrast, only few work has been done for solving the unconstrained sample mean problem. One algorithm is an (essentially optimal) dynamic program that finds global solutions of the unconstrained problem in exponential time [7]. A second algorithm is a heuristic, called adaptive DBA (ADBA) [28]. This algorithm uses a majorize-minimize algorithm (DBA) as a base-algorithm and iteratively refines subsequences to improve the solution quality.
The unconstrained problem is more challenging than the constrained one in the following sense: The constrained problem for length n ∈ N can be modelled as an optimization problem in a Euclidean space R n , where we have efficient local optimization techniques such as subgradient methods. 1 In contrast, the unconstrained problem can be modelled as a finite collection of constrained problems for sample means of length n = 1, . . . , n 0 where the maximum possible length n 0 of a sample mean is bounded by the Reduction Theorem [22]. The maximum length n 0 depends linearly on the lengths of the sample time series. A naive approach to address the unconstrained problem is to use a mean algorithm for the constrained problem and run it for all lengths n = 1, . . . , n 0 and return the solution with lowest Fréchet function value.
Apart from the fact that the naive approach is likely to be computationally intractable, there is the additional problem of how to choose an initial mean time series of desired length for each constrained mean algorithm. Current heuristics choose as initial point a random time series from the sample or the medoid of the sample. The initial length is therefore restricted to the lengths of the sample time series. One approach to obtain initial points of arbitrary length is generating them from a normal distribution. However, this method performed considerably worse in experiments [30].
In this work, we propose a generic average-compress (AC) algorithm for the unconstrained sample mean problem that addresses these issues. The AC algorithm repeatedly alternates between an averaging (A-step) and a compression (C-step). The A-step requires a time series as initial guess, minimizes the Fréchet function, and returns an approximate solution as output. The C-step compresses the approximation of the A-step to obtain an improved solution. The compressed solution of the C-step serves as initial guess of the A-step in the next iteration.
To ensure that the C-step returns an improved solution, a whole compression chain is efficiently computed and the best compression is selected. The AC algorithm can be initialized with a time series from the sample and iteratively seeks for promising solutions of shorter length using a compression technique. Thus, compression addresses both issues, finding a suitable length of the approximate solution for the unconstrained sample mean problem and finding suitable initial points for current state-of-the-art sample mean algorithms. In principle, any averaging algorithm and any compression method can be applied. Here, we propose a compression method that minimizes the squared DTW error between original and compressed time series. Empirical results suggest that the AC scheme substantially outperforms state-ofthe-art heuristics including ADBA. Our main contributions can be summarized as follows.
1. We propose a generic average-compress (AC) algorithm for the unconstrained sample mean problem under Dynamic Time Warping, 2. we empirically evaluate different configurations of the AC algorithm on 85 benchmark data sets from the UCR archive [10], using two compression techniques and multiple state-of-the art averaging heuristics, 3. we empirically show that on average the AC algorithm considerably outperforms current state-of-the-art averaging algorithms on the UCR benchmark data sets.
This article is organized as follows: Sect. 2 describes the AC algorithm. In Sect. 3 we present and discuss empirical results. Finally, Sect. 4 concludes with a summary of the main findings and an outlook for future research.

Average-compress algorithm
In this section, we develop an average-compress (AC) algorithm for approximately solving the unconstrained sample mean problem. To this end, we first introduce the DTW-distance (Sect. 2.1), the concept of a sample mean under DTW (Sect. 2.2), and compressions (Sect. 2.3). Thereafter, we describe the AC algorithm in Sect. 2.4. The first condition is the boundary condition and the second condition is the step condition of the DTW-distance. We denote the set of all warping paths of order m×n by P m,n . Suppose that p = ( p 1 , . . . , p ) ∈ P m,n is a warping path with points p l = (i l , j l ) for all l ∈ [ ]. Then p defines an expansion (or warping) of the time series x = (x 1 , . . . , x m ) and y = (y 1 , . . . , y n ) to the length-time series φ p (x) = (x i 1 , . . . , x i ) and ψ p (y) = (y j 1 , . . . , y j ). By definition, the length of a warping path satisfies max(m, n) ≤ ≤ m + n.

Dynamic time warping
The cost of warping time series x and y along warping path p is defined by where · denotes the Euclidean norm and φ p and ψ p are the expansions defined by p. The DTW-distance of x and y is is called an optimal warping path of x and y. By definition, the DTW-distance minimizes the Euclidean distance between all possible expansions that can be derived from warping paths. Computing the DTW-distance and deriving an optimal warping path is usually solved via dynamic programming [34,39].

Sample means under DTW
Note that S is a multiset that allows multiple instances of the same elements. A sample mean of S is any time series that minimizes the Fréchet function [16] F : serves as a measure of variability of S. Here, the search space U takes one of the following two forms: (i) U = T and (ii) U = T m . We refer to (i) as the unconstrained and to (ii) as the constrained sample mean problem. Note that the constrained formulation only restricts the length of the candidate solutions, whereas there is no length restriction on the sample time series to be averaged. A sample mean exists in either case but is not unique in general [22]. This result implies that F attains its infimum (has a unique minimum). However, computing a sample mean is NP-hard [8]. The implication is that we often need to resort to heuristics that return useful solutions within acceptable time.
We briefly describe two state-of-the-art algorithms for the constrained sample mean problem: a stochastic subgradient method (SSG) [35] and a majorize-minimize algorithm (DBA) [17,30]. For a detailed description of both algorithms, we refer to [35].
To present the update rule of both algorithms in a compact form, we introduce the notions of warping and valence matrix as proposed by [35]. Suppose that p ∈ P m,n is a warping initialize solution z ∈ T m 3: initialize best solution z * = z 4: repeat 5: reshuffle order of sample time series 6: for i ← 1 to N do 7: compute optimal warping path p i of z and x i 8: compute valence matrix V i of p i 9: compute warping matrix W i of p i 10: update solution z according to the rule 11: Suppose that z and x are time series of length |z| = m and |x| = n. Then W warps x onto the time axis of z. Each diagonal element v ii of V counts how many elements of x are warped to element z i . Stochastic Subgradient Algorithm. Subgradient methods for time series averaging have been proposed by [35]. Algorithm 1 outlines a vanilla version of the SSG algorithm with constant learning rate η. In practice, more sophisticated stochastic subgradient variants such as Adam [27] are preferred. The input of Algorithm 1 are a learning rate η, a length-parameter m of the constrained search space, and a sample x 1 , . . . , x N of time series to be averaged. The output is a time series with lowest Fréchet variation that has been encountered during optimization.
Majorize-Minimize Algorithm. Majorize-minimize algorithms for time series averaging have been proposed in the 1970s mainly by Rabiner and his co-workers with speech recognition as the primary application [32,42]. The early approaches fell largely into oblivion and where successively rediscovered, consolidated, and improved in a first step by Abdulla et al. [2] in 2003 and then finalized in 2008 by Hautamaki et al. [17]. In 2011, Petitjean et al. [30] reformulated, explored, and popularized the majorize-minimize algorithm by Hautamaki et al. [17] under the name DTW Barycenter Averaging (DBA). Algorithm 2 describes the DBA algorithm. It takes a length-parameter m and a sample of time series as input and returns the candidate solution of the last iteration as output. The DBA algorithm terminates after a finite number of iterations in a local minimum of the Fréchet function [35].

Compressions
Let x ∈ T be a time series of length n. A compression of x is a time series x of length m ≤ n that maintains some desirable problem-specific properties of x. By definition, x is also a compression of itself.
There are numerous compression methods such as principal component analysis, discrete Fourier transform, discrete wavelet transform, and many more. Here, we consider two simple methods: adaptive scaling (ADA) and minimum squared DTW error (MSE). Minimum Squared DTW Error Compression. The second compression method computes a time series of a given length such that the squared DTW error is minimized. Let x ∈ T be a time series of length n and let m < n. We call each x ∈ argmin dtw(x, z) 2 : z ∈ T m an MSE compression of x of length m. Observe that the MSE compression problem for x is the constrained sample mean problem of the sample S = {x}. Algorithm 4 outlines MSE compression.
It is not hard to see that for a compression x , an optimal warping path p between x and x warps every element of x to exactly one element of x , that is, φ p (x) = x. Thus, we can write p = (1, 1), . . . , (i 1 , 1), (i 1 + 1, 2), . . . , MSE compression is also known as adaptive piecewise constant approximation [9] and as segmentation problem [38]. It can be solved exactly via dynamic programming in O(n 2 m) time [6]. Moreover, the dynamic program allows to find all n compressions (for each length m = 1, . . . , n) in O(n 3 ) time by running it once for m = n (as it is done in Algorithm 4).
Interestingly, MSE compression is also related to one-dimensional k-means clustering. To see this relationship, consider an optimal warping path between the compression x and the original time series x as in Eq. (1). Then, the squared DTW error is minimal for x l = (x d l−1 +1 + · · · + x d l )/i l . Thus, finding an MSE compression x of length k can also be seen as a one-dimensional k-means clustering problem, where every cluster consists of a consecutive subsequence of elements in x. Indeed, the dynamic program in Algorithm 4 is the same as for one-dimensional k-means [41] (without previously sorting the elements in x).
To reduce the computational complexity, several heuristics and approximations for MSE compression have been proposed [9,25,38,43]. Also ADA compression can be regarded as a heuristic for MSE compression since it greedily averages two consecutive elements.
To conclude, with MSE compression, we consider an exact solution method to a sound compression problem and with ADA compression, we consider a fast heuristic. Among the various heuristics, we have chosen ADA compression because it has been successfully tested for improving approximate solutions of the constrained sample mean problem [30].

The average-compress algorithm
In this section, we assemble the pieces of the previous sections and propose a generic averagecompress (AC) algorithm for approximately solving the unconstrained sample mean problem. AC Algorithm. The AC algorithm alternates between averaging (A-step) and compression (C-step). For this purpose, any averaging algorithm and any compression method can be used. Algorithm 5 depicts the generic procedure. The input of the algorithm is a sample S of time series and an initial guess z ∈ T . It then repeatedly applies the following steps until termination: 1. A-step: approximate sample mean z ← average(S, z).

Evaluate solution:
(a) Select the shortest compression z * ∈ C(z) such that F(z * ) ≤ F(z ) for all z ∈ C(z). (b) If F(z * ) < F(z), set z ← z * and go to Step 1, otherwise terminate. Line 7 computes the complete compression chain C(z) that consists of all |z| compressions of z of lengths 1 to |z|. To accelerate the algorithm at a possible expense of solution quality, sparse compression chains can be considered.
In the following, we explain why and under which conditions compression is useful. To simplify our argument, we assume that AC uses an averaging algorithm for the constrained sample mean problem (such as SSG or DBA). In this case, the length m of the initial guess restricts the search space of AC to the set T ≤m of all time series of maximum length m. The choice of the length-parameter m via the initial guess is critical. If m is too small, the search space T ≤m may not contain an unconstrained sample mean. For a given sample S, the Reduction Theorem [22] guarantees the existence of an unconstrained sample mean of a length at most m S = x∈S |x| − 2(|S| − 1). Consequently, we can safely constrain the search space to T ≤m S for solving the unconstrained sample mean problem. Then, a naive approach to minimize the Fréchet function on T ≤m S is to solve m S constrained sample mean problems on T m S , . . . , T 1 and then to pick the solution with lowest Fréchet variation. When using state-of-the-art heuristics for the m S constrained problems, the naive approach is computationally infeasible.
The purpose of compression is to substantially accelerate the intractable naive approach at the expense of solution quality. Instead of solving all m S constrained problems, the AC A-Step: z ← average(X , z) 7: C-Step: C(z) ← compress(z) 8: //*** Evaluate solution 9: z ← argmin F(z ) : z ∈ C(z) ∪ {z} 10: if F(z) < f * or (F(z) = f * and |z| < * ) then 11: f * ← F(z) 12: l * ← |z| 13: z * ← z 14: until convergence 15: return z * algorithm uses compressions to select a few promising search spaces is compressed in order to determine the next search space T m i . The length-parameter m i of the next search space T m i corresponds to the length of the compression z i of z i−1 with lowest Fréchet variation. Obviously, this idea only accelerates the naive approach if the length m i of the best compression is substantially smaller than the length m i−1 of the previous solution.
The theoretical upper bound m S provided by the Reduction Theorem [22] is usually very large such that existing state-of-the-art heuristics for solving the constrained problem on T ≤m S are computationally intractable. In this case, also the AC algorithm using such a heuristic would be infeasible. However, empirical results on samples of two time series of equal length n suggest that the length of an unconstrained sample mean is more likely to be less than n [7]. Similar results for larger sample sizes are unavailable due to forbidding running times required for exact sample means. For solving constrained sample mean problems, it is common practice to choose m within the range of the lengths of the sample time series [30]. Within this range, experimental results showed that an approximate solution of a constrained sample mean can be improved by reducing its length using adaptive scaling [30]. These findings suggest to choose the length-parameter m within or slightly above the range of lengths of the sample time series.

Experiments
Our goal is to assess the performance of the proposed AC algorithm. For our experiments, we use the 85 data sets from the UCR archive [10]. Appendix 1 summarizes the parameter settings of the mean algorithms used in these experiments.

Comparison of ADA and MSE
We compared the performance of ADA and MSE as compression subroutines of the AC algorithm. We applied the following configurations: Column I refers to the number of iterations of the repeat-until loop of the AC Algorithm. Compression schemes with * iterations run until convergence. We applied DBA and the four AC algorithms to approximate the class means of every UCR training set. 2 To assess the performance of the mean algorithms, we recorded the percentage deviations, ranking distribution, and space-saving ratios. Here, we used the solutions of the DBA algorithm as reference. The percentage deviation of a mean algorithm A is defined by where z DBA is the solution of the DBA algorithm and z A is the solution of algorithm A. Negative (positive) percentage deviations mean that algorithm A has better (worse) Fréchet variation than DBA. The ranking distribution summarizes the rankings of every mean algorithm over all samples. The best (worst) algorithm is ranked first (last). The space-saving ratio of algorithm A is A positive (negative) space-saving ratio means that the solution z A is shorter (longer) than z DBA . Table 1 summarizes the results. The top table shows the average, standard deviation, minimum, and maximum percentage deviations from the Fréchet variation of the DBA algorithm (lower is better). The table in the middle shows the distribution of rankings and their corresponding averages and standard deviations. The best (worst) algorithm is ranked first (fifth). Finally, the bottom table shows the average, standard deviation, minimum, and maximum space-saving ratios (higher is better).
All AC variants improved the solutions of the DBA baseline by 4.6% to 7.0% on average and 45%(±2%) in the best case. By construction, an AC solution is never worse than a DBA solution. The best method is DBA-MSE with average rank 1.0 followed by DBA-MSE 1 and DBA-ADA with average ranks 2.4 and 2.5, respectively. These three methods clearly outperformed DBA-ADA 1 proposed by Petitjean et al. [30].
On the majority of data sets the AC algorithm converged after two or three iterations to a local optimum. In rare cases, four iterations were required. The number of iterations does neither correlate with the length of the time series nor with the sample size. The main improvement of DBA-MSE and DBA-ADA occurs at the first iteration. The first compression appears to be most aggressive, because it collapses flat regions. Subsequent averaging steps slightly shift the compressed candidate solutions. These shifts may create new flat regions, which can be further compressed to obtain better solutions.
We also considered the lengths of the approximated means. Recall that all mean algorithms started with the same initial guess (medoid). The length of a DBA solution corresponds to

Comparison of AC algorithms
The goal of the second experiment is to compare the performance of the following mean algorithms:  Table 2 summarizes the results using the same legend as in Table 1. The percentage deviations and rankings suggest that the three AC variants DBA-MSE, SSG-MSE, and ADBA-MSE performed substantially better than the corresponding base algorithms DBA, SSG, and ADBA, respectively. The SSG-MSE algorithm performed best with an average rank of 1.4, followed by ADBA-MSE (2.9), SSG (3.1), and DBA-MSE (3.3). Interestingly, ADBA performed worse than DBA-MSE. Both, ADBA and DBA-MSE, are based on the DBA algorithm. The difference between both algorithms is that ADBA compresses and expands selected subsequences of the current DBA solution, whereas DBA-MSE only compresses the current DBA solution. The results indicate that simple MSE compression on the entire sequence appears to be a better strategy than ADBA's compression and expansion schemes on selected subsequences. Notably, SSG performed best among the three base averaging algorithms DBA, SSG, ADBA, and performed even better than DBA-MSE. These results are in contrast to those presented in [28], where ADBA outperformed SSG (and also DBA). Our findings confirm that the performance of SSG substantially depends on a careful selection of an optimizer (such as Adam) and a proper choice of the initial learning rate.
Next, we examine the length of the solutions. Note that SSG also does not alter the length of its initial guess such that ρ ss (DBA) = ρ ss (SSG) = 0. The bottom table shows that MSE compression schemes reduce the length of the solutions obtained by their corresponding base algorithm (DBA, SSG, and ADBA). The space-saving ratios of the AC variants are roughly independent of the particular base algorithm for mean computation (0.43-0.45). Notably, the base algorithm ADBA is more likely to compress rather than to expand the DBA solutions. This finding is in line with the hypothesis that an exact mean is typically shorter than the length of the sample time series [7].

Qualitative analysis of mean algorithms
The goal of this section is to qualitatively analyze the behavior and phenomena behind the different types of mean algorithms. For this, we considered DBA, ADBA, and the AC variant SSG-MSE relative to an exact dynamic program (EDP) proposed by [7]. Since the sample mean problem is NP-hard, we only considered a sample with two sample time series of length 24 from the Chinatown data set [12]. The two sample time series slightly differ in their amplitudes and on a high abstraction level, they have the following features in common: Both start with a wide valley followed by a peak, a flat plateau-like valley at a high altitude until they finally end with a descent. Figure 2 shows the sample time series and the sample means returned by the four algorithms. The four algorithms differ in the level of feature abstraction and in susceptibility of spurious features, whereby lower level feature representations are more prone to spurious features than higher level ones: The EDP exhibits the common shape of both sample time series. In addition, it filters out variations of speed by condensing the mean to length 16. The SSG-MSE algorithm more aggressively condenses a solution than EDP resulting in a more compact and higher level description of a mean of length 13. In contrast to EDP, the AC variant SSG-MSE has smoothed out the flat plateau-like valley. We note that this "over-compression" is shape-dependent but not length-dependent. As indicated by Fig. 2, over-compression is more likely for flat and less likely for steep regions. The solution of ADBA more moderately condenses a solution than SSG-MSE and EDP resulting in a lower level representation of length 18. In addition, ADBA includes a spurious plateau at the beginning that occurs only in the upper sample time series. Finally, DBA aims at capturing the common features of both sample time series with respect to a predefined length (here 24). The resulting solution contains more low level features than the other approaches and includes spurious features which occur in only one of both time series. Finally, we hypothesize that spurious features may also occur in exact solutions when two sample time series do not share many common features. Figure 2 shows the error of each mean algorithm as the percentage deviation of their Fréchet variations from the minimum Fréchet variation. The Fréchet variation measures the amount of dispersion of a sample of time series. Such a measure is, for example, important in evaluating k-means clustering using validation indices based on the Fréchet variation for each cluster. The errors of the heuristics differ substantially with EDP ranked first followed by SSG-MSE (1.0%), ADBA (10.1%), and DBA (19.3%). We hypothesize that low level and spurious features could result in erroneous measures of dispersion. These errors then propagate to pattern recognition methods based on time series averaging such as k-means clustering.

Application: k-means clustering
In this experiment, we investigated how the quality of a mean algorithm affects the quality of a k-means clustering. Let S = {x 1 , . . . , x n } ⊆ T be a set of n finite time series. The goal of k-means is to find a set Z = {z 1 , . . . , z k } of k centroids z j ∈ T such that the k-means error is minimized. We used DBA, SSG, ADBA, DBA-MSE, SSG-MSE, and ADBA-MSE for computing the set Z of centroids. We applied the six variants of k-means to 70 UCR data sets and excluded 15 UCR data sets due to overly long running times (see Appendix 1). We merged the prespecified training and test sets. The number k of clusters was set to the number of classes and the centroids were initialized by the class medoids. Table 3 summarizes the results. The top table presents the average, standard deviation, minimum, and maximum percentage deviations from the respective minimum k-means error (lower is better). The percentage deviation of k-means algorithm A for data set D is defined by where Z A is the set of centroids returned by algorithm A and Z D is the best solution obtained by one of the six k-means algorithms. The bottom table shows the distribution of rankings and their corresponding averages and standard deviations. The best (worst) algorithm is ranked first (sixth).
The results show that the AC approach substantially improved all k-means variants using one of the base averaging methods (DBA, SSG, ADBA). Notably, SSG-MSE performed best with an average percentage deviation of 1.5% and an average rank of 1.4, followed by ADBA-MSE (3.8% and 2.9). The average percentage deviations of DBA and DBA-MSE are substantially impacted by the results on a single data set (DiatomSizeReduction). Removing the DiatomSizeReduction data set yields an average percentage deviation of 13.2 for DBA and 4.6 for DBA-MSE, whereas the other average percentage deviations remain unchanged up to ±0.1%. These findings confirm the hypothesis raised by Brill et al. [7] that better mean algorithms more likely result in lower k-means errors.

Conclusion
We formulated a generic average-compress algorithm for the unconstrained sample mean problem in DTW spaces. Starting with an initial guess of sufficient length, the AC algorithm alternates between averaging and compression. In principle, any averaging and any compression algorithm can be plugged into the AC scheme. The compression guides the algorithm to promising search spaces of shorter time series. This approach is theoretically justified by the Reduction Theorem [22] that guarantees the existence of an unconstrained sample mean in a search space of bounded length. Experimental results show that the AC algorithm substantially outperforms state-of-the-art heuristics for time series averaging. In addition, we observed that better averaging algorithms yield lower k-means errors on average. Open research questions comprise application of the AC scheme to the empirical analysis of alternative compression methods for the AC algorithm and reducing its computational effort.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 Hyperparameter settings
In all experiments, we selected the sample medoid as initial guess of a mean algorithm. The DBA algorithm terminated after convergence and latest after 50 epochs (cycles through a sample). The ADBA algorithm terminates subsequence optimization when the sum of the scaling coefficients changes its sign and latest after 50 iterations. The SSG algorithm terminated after 50 iterations without observing an improvement and latest after max (50, 5000/n) epochs. As optimization scheme, SSG applied Adam [27] with β 1 = 0.9 as first and β 2 = 0.999 as second momentum. To cope with the problem of selecting an initial learning rate, we used the procedure described in Algorithm 6. The input is a sample S of size n. The output is the best solution found. The algorithm terminates if the solution did not improve for two consecutive learning rates and latest if √ n/2 i ≤ 10 −6 .

A.2 Data sets excluded from k-means experiments
The following list contains all UCR data sets excluded from k-means clustering due to computational reasons: