A time series distance measure for efficient clustering of input output signals by their underlying dynamics

TL;DR

This paper introduces an extension to the Martin cepstral distance for efficient clustering of input/output time series based on their underlying dynamics, outperforming traditional methods in computational efficiency.

Contribution

The paper presents a novel extension to the Martin cepstral distance that enables fast and effective clustering of dynamical system-generated time series.

Findings

The new distance measure performs comparably to explicit system modeling.

It is significantly faster with a complexity of O(N logN).

Effective for clustering electrical circuit data.

Abstract

Starting from a dataset with input/output time series generated by multiple deterministic linear dynamical systems, this paper tackles the problem of automatically clustering these time series. We propose an extension to the so-called Martin cepstral distance, that allows to efficiently cluster these time series, and apply it to simulated electrical circuits data. Traditionally, two ways of handling the problem are used. The first class of methods employs a distance measure on time series (e.g. Euclidean, Dynamic Time Warping) and a clustering technique (e.g. k-means, k-medoids, hierarchical clustering) to find natural groups in the dataset. It is, however, often not clear whether these distance measures effectively take into account the specific temporal correlations in these time series. The second class of methods uses the input/output data to identify a dynamic system using an identification scheme, and then applies a model norm-based distance (e.g. H2, H-infinity) to find out which systems are similar. This, however, can be very time consuming for large amounts of long time series data. We show that the new distance measure presented in this paper performs as good as when every input/output pair is modelled explicitly, but remains computationally much less complex. The complexity of calculating this distance between two time series of length N is O(N logN).

Equations

$$d_E(y_m,y_n)=\sqrt{\sum _{k=0}^{N_m}{\left(y_m(k)-y_n(k)\right)}^2}.$$

$$d_{DTW}\left(y_1,y_2\right)=min\left\{\sqrt{\sum _{k=1}^K{w}_k}\right\}.$$

$$||\mathcal{H}|{|}_2=\sqrt{\frac{1}{2\pi }{\int }_{-\pi }^{\pi }Tr\left\{{\mathcal{H}}^H({\text{e}}^{i\omega })\mathcal{H}({\text{e}}^{i\omega })\right\}d\omega },$$

$$||\mathcal{H}|{|}_{\infty }=\underset{\omega \in [0,\pi [}{\max }|\mathcal{H}({\text{e}}^{i\omega })|.$$

$$c_y={\mathcal{F}}^{-1}(\log ({\Phi }_y)),$$

$$||\mathcal{H}|{|}_C=\sum _{k=0}^{N}k{\left(c_y(k)\right)}^2.$$

$$d_C(y_i,y_j)=\sum _{k=0}^{\max \{N_i,N_j\}}k{\left(c_{y_i}(k)-c_{y_j}(k)\right)}^2,$$

$$c_y(k)=c_u(k)+c_h(k).$$

$$d_{C_e}$$

$$\sum _{k=0}^{\min \{N_i,N_j\}}k{\left(c_{h_i}(k)-c_{h_j}(k)\right)}^2.$$

Full text

Traditionally, two ways of handling the problem are used. The first class of methods employs a distance measure on time series (e.g. Euclidean, Dynamic Time Warping) and a clustering technique (e.g. k-means, k-medoids, hierarchical clustering) to find natural groups in the dataset. It is, however, often not clear whether these distance measures effectively take into account the specific temporal correlations in these time series. The second class of methods uses the input/output data to identify a dynamic system using an identification scheme, and then applies a model norm-based distance (e.g. $H_{2}$, $H_{\infty}$) to find out which systems are similar. This, however, can be very time consuming for large amounts of long time series data.

We show that the new distance measure presented in this paper performs as good as when every input/output pair is modelled explicitly, but remains computationally much less complex. The complexity of calculating this distance between two time series of length $N$ is $\mathcal{O}(N\log{N})$.

I INTRODUCTION

Time series clustering is an important topic in modern research. State-of-the-art clustering methods of other data types are often not suited for this high-dimensional, temporally correlated data structure. Clustering is the task of finding groups with similar elements in a dataset and consists of three components: a similarity measure based on relevant data features, a clustering algorithm and an evaluation criterion. While the latter two components might carry over, defining a good distance measure is a difficult problem, especially if one is interested in the dynamics of the generating dynamical system of the time series.

Representing the time series as single-input single-output (SISO) linear time invariant (LTI) deterministic dynamical systems further generates problems of its own, as the contributions of the input signal and the impulse response of the system are convolved in the time domain. It is thus not intuitively clear how these two contributions can be separated, for example when one is interested only in the dynamics of the system and not in the specific input signal.

This problem grows ever more relevant as large scale big data time series problems grow more prevalent in areas like finance, medicine, or the industrial internet of things, where clustering is important in tasks like anomaly detection [7, 12]. A typical industrial problem contains several hundred sensors per machine, tens of machines per plant, and several plants per industrial player, collecting data every few seconds, for months or even years of operation time. This results in datasets of several million time points for thousands of series. Clustering techniques should thus scale well.

In Section II we look at state-of-the-art clustering methods for time series from two perspectives, starting from a dataset containing input/output time series pairs, generated by different SISO LTI dynamical systems. From a machine learning point of view, we use an automated clustering method with an off-the-shelf time series distance such as the Euclidean distance or Dynamic Time Warping (DTW). From a system identification point of view, we apply norms such as the $H_{2}$ or $H_{\infty}$ norm to compare systems estimated from the data. We find that these techniques either are very fast, but give poor results, or perform well, but are computationally expensive.

Next, in Section III, we look at the Martin cepstral distance [3, 8], which combines insights from systems theory into a distance measure that can be computed on the raw data. This metric was defined for SISO ARMA models (i.e. LTI models that use white noise as an input signal).

The main contribution of this paper is an extension of the cepstral distance measure, that incorporates deterministic input signals, and allows to calculate distances between a broader class of SISO LTI dynamical systems. It thus allows to cluster time series by dynamics, but remains computationally much simpler than explicitly estimating models.

Subsequently, we apply this new distance measure in Section IV to an application on electrical circuits, where we generate a dataset consisting of input/output signal pairs, and the problem is to identify which data belong to which generating system. Finally, we conclude the paper and provide some paths for future research in Section V.

II EXISTING METHODS

Existing methods to cluster time series employ a clustering technique, together with some distance measure. The author of [6] discerns three types of distance measures: measures based on raw data, measures based on features of the time series and measures based on models. For the scope of this paper, we will focus on the first and the latter (as the distance measure we propose combines elements of these two broad classes). We present two raw data distance measures, the Euclidean metric and the Dynamic Time Warping metric [5], and two model-based distance measures, connected to the $H_{2}$-norm and the $H_{\infty}$-norm. In the next section, we will introduce and extend the cepstral distance [3, 8], which combines the efficiency of the raw data distance measures with the insight in generative dynamics of the model norms, and thus has representations both as a raw data distance and as a model-based one.

II-A Raw Data Distance Measures

In what follows we will define $u_{m}$ to be the input signal of the $m$-th element of a dataset, $y_{m}$ is the corresponding output signal and $u_{m}(k)$ or $y_{m}(k)$ is the value at timepoint $k$ of respectively the input and output of the $m$-th element of the input/output dataset. Time series from element $m$ start at $k=0$ and end at $k=N_{m}$. The system that generated an output from a given input will be called the generating (dynamical) system.

II-A1 Euclidean Distance

Definition 1**.**

The Euclidean distance, $d_{E}(\cdot,\cdot)$ treats the time series as a vector, and applies the element-wise Euclidean vector distance between two time series of same length $N_{m}$, defined as

[TABLE]

Advantages

• •

The Euclidean distance is easy to calculate, allowing for very efficient computation and clustering.

• •

No system identification step is needed.

Disadvantages

• •

There is no clear link between this distance measure and the generating system.

• •

This measure treats the time series as a vector, and ignores the temporal correlations in the data.

• •

This measure does not allow to compute distances between time series of different length.

• •

This measure does not take the input into account.

II-A2 Dynamic Time Warping

Dynamic Time Warping (DTW) [5, 11] is an algorithm that tries to locally align time series, by warping them such that the Euclidean distance between the warped time series is minimal. Mathematically, this warping, and the measure that is found in this way, can be described as follows.

Given two output signals, $y_{1}$ and $y_{2}$, of length $N_{1}$ and $N_{2}$ respectively, a matrix $M$ is constructed, where the $(l,m)$-th element of $M$ is defined as $M_{(l,m)}=(y_{1}(l)-y_{2}(m))^{2}$. A warping path, $W=w_{1},w_{2},\ldots,w_{k},\ldots,w_{K}$ is then defined, with each $w_{k}=\left(M_{(l,m)}\right)_{k}$ an element of matrix $M$ and $max(N_{1},N_{2})\leq K<N_{1}+N_{2}-1$.

The path is subject to the boundary conditions $w_{1}=M_{1,1}$ and $w_{K}=M_{N_{1},N_{2}}$ (i.e. the path starts in one corner of the matrix and ends in the opposite one), has to be continuous, in such a way that two consecutive elements $w_{k}$ and $w_{k+1}$ are maximally one column and one row apart, and has to be monotonously increasing in its indices, i.e., that in going from $w_{k}$ to $w_{k+1}$, column nor row number can decrease.

Definition 2**.**

We are now interested in the warping path $W_{DTW}$ that minimizes the cost function

[TABLE]

The sum over this path is then the DTW distance between the time series.

Though this algorithm is computationally expensive due to the combinatorial nature of the problem, several lower bounds have been devised that can be implemented efficiently. In what follows, we use the Keogh Lower Bound [5] as an efficient approximation to the DTW distance.

Advantages

• •

The DTW distance takes into account (part of) the local temporal correlations.

• •

No system identification step is needed.

• •

Lower bounds on the distance are reasonably efficient.

• •

This measure allows to calculate distances between time series of different length.

Disadvantages

• •

There is no clear link between this distance measure and the generating system.

• •

The DTW distance as such is expensive to calculate.

• •

This measure does not take the input into account.

II-B Model-based Distance Measures

We use the same notation as in subsection II-A. The generating system of the input/output pair $(u_{m},y_{m})$ will be denoted by $M_{m}$, and its corresponding transfer function will be written $\mathcal{H}_{m}$. Based on a model norm $||\cdot||$, the distance between two models $M_{i}$ and $M_{j}$ is defined as $||\mathcal{H}_{i}-\mathcal{H}_{j}||$.

II-B1 $H_{2}$-norm

Definition 3**.**

The $H_{2}$-norm, $||\mathcal{H}||_{2}$, of a discrete-time system $M$ with transfer function $\mathcal{H}$ is defined as

[TABLE]

where $Tr\{\}$ denotes the trace, the superscript $\cdot^{H}$ denotes the Hermitian conjugate and $i$ denotes the imaginary unit.

The $H_{2}$-norm can be seen as the root-mean-square of the system response to a normalized white noise input. It is thus a measure of the power, or steady-state variance of this response. The $H_{2}$-norm will be infinite for unstable systems.

Advantages

• •

The $H_{2}$-norm provides a physically interpretable way to characterize underlying dynamics of time series.

• •

This norm allows to calculate distances between time series of different length.

• •

This norm takes the input data into account.

Disadvantages

• •

A system identification procedure is needed, which is both difficult to automate and often computationally expensive (at least more expensive than the raw data measures).

II-B2 $H_{\infty}$-norm

Definition 4**.**

The $H_{\infty}$-norm, $||\mathcal{H}||_{\infty}$, of a discrete-time system $M$ with transfer function $\mathcal{H}$ is calculated as

[TABLE]

This norm thus measures the maximal gain of the frequency response and is called the gain of the system. It becomes infinite for systems with poles on the unit circle.

Advantages

• •

The $H_{\infty}$-norm provides a physically interpretable way to characterize underlying dynamics of time series.

• •

This norm allows to calculate distances between time series of different length.

• •

This norm takes the input data into account.

Disadvantages

• •

A system identification procedure is needed, which is both difficult to automate and often computationally expensive (at least more expensive than the raw data measures).

III CEPSTRAL DISTANCE

In this section we take a closer look at an insightful distance measure on ARMA models, which can be interpreted both as a raw data distance measure and as a model norm: the Martin cepstral norm [3, 8]. We first give a very concise review of the cepstral norm in the stochastic case, then proceed with an extension that allows us to incorporate information about the deterministic input signal.

III-A Original Cepstral Norm

Based on the power spectral density, $\Phi_{y}$, of a signal $y$, we can define its power cepstrum, $c_{y}$ as

[TABLE]

where $\mathcal{F}^{-1}$ denotes the inverse Fourier transform. This produces a series of coefficients, $c_{y}(k)$, with integer $k\in[0,N]$, where $N$ denotes the length of time series $y$.

Definition 5**.**

The cepstral norm, $||\mathcal{H}||_{C}$, of model $M$ with transfer function $\mathcal{H}$, and output $y$ is defined as

[TABLE]

For ARMA models it was proven in [3] that there are multiple methods to calculate this norm: it can be derived from the subspace angles of the output Hankel matrices of the generating system, from the mutual information of the output space of a system, and from a combination of poles and zeros of the transfer function of the model. Moreover, equation (6) allows us to calculate the norm straight from raw data, without the need to identify the underlying systems. We can thus connect the cepstral norm to a raw data distance measure in the following sense:

Definition 6**.**

The cepstral distance, $d_{C}(y_{i},y_{j})$, between two time series, $y_{i}$ and $y_{j}$, is defined as

[TABLE]

where $\max\{N_{i},N_{j}\}-\min\{N_{i},N_{j}\}$ zeros are added at the end of the cepstrum of length $\min\{N_{i},N_{j}\}$.

Advantages

• •

The cepstral distance has an interpretation in terms of the generating model of the time series.

• •

The cepstral distance is easy to calculate, allowing for very efficient computation and clustering.

• •

No system identification step is needed.

• •

This measure allows to calculate distances between time series of different length.

Disadvantages

• •

This distance measure can only take information coming from a stochastic input into account.

III-B Extended Cepstral Distance

The cepstrum, defined in the previous section, finds its roots in homomorphic signal processing [9, Chapter 10]. In this type of processing, the original time series data, which often involves complex multiplicative operators like convolutions, is mapped, through a non-linear mapping, to a different domain, that allows for linear filtering. The cepstrum, as in equation (5), is a good example. The convolution in the time domain changes into a multiplication by calculating the power spectral density. Applying a logarithmic transformation then turns the multiplication in frequency domain into an addition. Finally, the inverse Fourier transform takes the problem back to (a transformed version of) the time domain. Equation (5) is thus effectively a method to transform the convolution into an addition.

This allows us to take the output, and separate the contributions from the input signal (which was the main disadvantage left in the cepstral distance, see subsection III-A) and the impulse responses of the system. Indeed, defining the cepstrum coefficients of the input signal $u$ as $c_{u}(k)$, and the contribution to the cepstrum coefficients of the transfer function $\mathcal{H}$ as $c_{h}(k)$, we can write

[TABLE]

Based on input/output signal pairs, we now have a measure of the underlying generating system dynamics by looking at $c_{h}(k)=c_{y}(k)-c_{u}(k)$.

Definition 7**.**

The extended cepstral distance, $d_{C_{e}}((y_{i},u_{i}),(y_{j},u_{j}))$, between two input/output pairs of time series, $(y_{i},u_{i})$ and $(y_{j},u_{j})$, with respective transfer functions $\mathcal{H}_{i}$ and $\mathcal{H}_{j}$, is defined as

[TABLE]

Note that, for now, this distance measure does not have the whole theoretic framework with connections to subspace angles, mutual information and generating system parameters.111These theoretical equivalences will be researched and most of them proven to carry over in a forthcoming paper, where we will also try to connect the extended cepstral distance to an extended cepstral model norm. However, it is clear that the $c_{h}(k)$ can only come from the generating system dynamics, and thus the distance measure tells us something about these systems, even if it is still unclear what exactly is measured.

We propose this extended cepstral distance as a way to efficiently cluster input/output data by their generating dynamics.

Advantages

• •

The extended cepstral distance is linked to the generating model of the time series.

• •

The extended cepstral distance is easy to calculate, allowing for very efficient computation and clustering.

• •

No system identification step is needed.

• •

This measure allows to calculate distances between time series of different length.

• •

This measure takes the input into account.

Disadvantages

• •

The interpretation of the measure in terms of system parameters and properties is not immediately clear, thus the theoretical framework of the original cepstral distance does not carry over trivially.

IV APPLICATION ON ELECTRICAL CIRCUITS

IV-A Simulation Set-Up

To test the proposed techniques, we simulate data coming from electrical circuits. We start out by modelling two circuits with the same topology, but different values for the R, L, and C components. The topology was taken from a course on linear physical systems analysis [2]. The network topology and the values of the components are shown in Figure 1. The input of the system is the current $i_{u}$, the output is the voltage over $L_{2}$, $e_{y}$. State-space models of order 3 are then written down for these networks.

We provide both systems with 200 different input signals (100 outputs of LTI models of order 15, 50 multisine waves corrupted by Gaussian white noise with standard deviation of 0.1 and 50 white noise signals), and measure the output signals. This generates a dataset of 400 input/output signal pairs (200 inputs times 2 models). The question at hand is whether we can use this input/output data, and only this data, to determine which pairs were generated by the same system, i.e. cluster the dataset in two groups, defined by the generative dynamics.

We will do this using the distance measures defined in section II and subsection III-A, keeping in mind that we use the Keogh Lower Bound [5] as an efficient approximation to DTW. We then compare to the technique developed in subsection III-B. There, the power spectral density is estimated by using Welch’s method [13], which provides a stable approximation222Note that, for longer time series (i.e. $2^{10}$ and beyond), the Fast Fourier Transform [1] provides a clean enough output to work on. We could thus speed up the algorithm even further for longer series. of the Fourier transform for short time series. In the Appendix, we give a pseudo-code overview of how the distance measure is calculated, as well as a link to a minimal working example of the simulations discussed and a complexity analysis of the algorithm.

The performance of these simulations will be measured by the Adjusted Rand Index (ARI) [4, 10], which is a similarity measure between partitions. The ARI compares two partitions, $S_{1}$ and $S_{2}$, by calculating the ratio of pairs that have the same partitioning status (i.e. belonging to the same partition or not) in both $S_{1}$ and $S_{2}$ to the total amount of data pairs, then adjusting the resulting ratio by subtracting the expected value, to account for guessing (i.e. a partitioning that is the result of random guessing is assigned an ARI of 0). An ARI of 1 corresponds to perfectly similar partitions.

We compare the partitions generated by a hierarchical clustering method, cut-off at two clusters, using distance matrices generated by the different distance measures of section II and section III versus the ground truth (i.e. the time series was generated by the system with parameters $S_{1}$ or with parameters $S_{2}$, as in Figure 1).

IV-B Results

The results for the set-up in the previous subsection are shown in Figure 2, which shows the average and standard deviation for the ARI of the simulation results, and Figure 3, which shows the average and standard deviation for the execution time of the simulations.

It is clear that the extended cepstral distance gives the best results. In fact, it manages to cluster the simulated input/output pairs perfectly every time. This is, of course, to be expected, as this distance measure was tailored specifically to take into account the dynamics of the underlying model333We redid the experiments for generating systems of higher order, and the extended cepstral distance still performed best. Results were omitted., and nothing but those dynamics. The reasons why it performs better than the other measures will be explained in what follows, and we will again use the distinction between raw data and model-based distances measures from Section II.

IV-B1 Raw Data Distance Measures

The reason why the other raw data distance measures do not perform well on the problem at hand, is because they do not take into account the information from the input signal. Indeed, the dynamics of the output are dominated by the input, due to the way the inputs were designed (i.e. the models generating the inputs are of higher order than the models describing the electrical circuits). The other distance measures are thus dominated by contributions coming from the input to cluster the time series, as they cannot separate the different contributions.

If we only use white noise inputs, we see, on the left hand side in Figure 4, that the original cepstral distance performs better.444In fact, the original and extended cepstral distance are equivalent in this case. Indeed, the cepstrum of white noise is only non-zero in its zeroth component, which is not taken into account in the sum in equations (7) and (9), which coincide in that case. The Euclidean and DTW distances still do not deliver good results when detecting the difference in dynamics.

There is thus no hope to achieve better results by taking the input signal into account in the case of the Euclidean distance or the DTW distance. Indeed, the distances look at the shape of the signal, rather than its generative dynamics. DTW is better at this job [5], but, as we can see from Figure 3, also has a big disadvantage: it takes a lot of time to compute, especially for long time series, where it even surpasses the model-based distance measures in computation time.

Based on these results, the extended cepstral distance is thus preferred to cluster input/output signals based on the dynamics of their generating models.

IV-B2 Model-based Distance Measures

The model-based distance measures show better results than the raw data distance measures, and this again is to be expected. Indeed, the model-based measures take the input information into account and thus manage to peel out the information on the system that generated the input/output pair. However, since a priori we have no information on the order of the underlying system, we arbitrarily have to set a model order. In this case, we estimated transfer functions of order 5. If we share the information on the correct model order (3) with the system identification algorithm, the performance of the model norms increases, as on the right hand side of Figure 4.

There exist, of course, schemes to determine appropriate model orders, and more effort can be put in correctly identifying the underlying model. However, as can be seen from Figure 3, the model norm techniques are already several orders of magnitude slower than the extended cepstrum distance measure. For problems concerning large amounts of long input/output-pairs, as can be found in realistic problems in process industry (see, for example, [7], where more than 250 sensors make a measurement every 5 minutes for 6 months), this becomes highly impractical.

The extended cepstral distance is thus preferred over explicitly identifying systems, because of both being easier to automate, and taking less time to compute.

V CONCLUSION AND FURTHER RESEARCH

We have devised a distance measure that is as insightful as a model norm-based distance, yet remains computationally much simpler than explicitly estimating models. It allows to meaningfully cluster large input/output signal pair datasets based exclusively on the dynamics of the generating systems. We have tested it on a simulation of data coming from electrical circuits, where we started from two electrical circuits with a current as input and a voltage difference over an inductor as output. We provided both circuits with 200 different inputs, resulting in 400 input/output pairs.

We then showed that the proposed measure performs as well as model-based distances on estimates of the generative systems, but is much easier to calculate and that other distance measures (Euclidean, DTW) perform much worse.

We furthermore show that, in the stochastic input case, the extended distance proposed in this paper reduces to the original cepstrum distance, which was proven ([3, 8]) to be equivalent to a model norm. This gives hope that the extended distance could also be linked to a model norm. Research that looks into this link is currently under way and will be discussed in a forthcoming paper.

The results indicate the extended cepstral distance measure does a good job of capturing the dynamics of input/output pairs. An application to a real-life dataset is needed to validate the effectiveness in practice, but for the simulated problem at hand, the distance measure succeeded in perfectly distinguishing different dynamics based on raw data alone.

APPENDIX

A pseudo-code overview of the algorithm is shown in Algorithm 1. A minimal working example of the simulations performed in Section IV is available on GitHub.555https://github.com/Olauwers/Extended-Cepstral-Distance

Calculating the extended cepstral distance amounts to estimating the power spectral density of both input and output by Welch’s method [13] (employing the FFT, which is of $\mathcal{O}(n\log{n})$, with $n$ the length of the windows considered in Welch’s method), taking the logarithm of the resulting vector, and then applying an inverse Fourier transform (employing the IFFT, running in $\mathcal{O}(N\log{N})$ time, with $N$ the length of the time series) on them. In the end, we then apply a weighted Euclidean distance on the results.

The complexity of calculating the extended cepstral distance between two time series is thus $\mathcal{O}(N\log{N})$, with $N$ the length of the time series.