Change-Point Detection with Feature Selection in High-Dimensional Time-Series Data

1 Introduction

Change-point detection, which is the problem of detecting abrupt changes in time-series data, is attracting a lot of attention in the artificial intelligence and data mining communities [Basseville et al., 1993; Brodsky and Darkhovsky, 1993; Kifer et al., 2004], and there are various types of real-world applications such as fraud detection in cellular systems [Murad and Pinkas, 1999], intrusion detection in computer networks [Yamanishi et al., 2000], irregular-motion detection in vision systems [Ke et al., 2007], music segmentation [Desobry et al., 2005], and sentiment analysis from Twitter data [Liu et al., 2013]. Recently, the problem of change-point detection from high-dimensional time-series data such as music and multi-sensor data has been attracting increasing attention [Desobry et al., 2005; Liu et al., 2013]. Compared with traditional change-point detection problems, the number of features d tends to be larger than the number of data points n, and it includes a large number of noise features. This makes the change-point detection problem from high-dimensional timeseries data challenging.

An effective change-point detection approach would be to use the divergence between probability distributions of data in the past and data in the future at time t. Specifically, we regard the time t point as a change-point if the divergence between two distribution is significantly large. Various change-point detection methods have been proposed based on this concept including generalized likelihood ratio (GLR) and cumulative sum approaches [Gustafsson, 1996; Basseville et al., 1993]. In these approaches, the logarithm of the likelihood ratio between two probability distributions is used as a measure of change-point detection, where each probability density is estimated independently by density estimation. However, since density estimation is known to be a difficult problem [H¨ardle et al., 2004; Huang et al., 2007], density estimation based approaches tend to perform poorly. Moreover, since high-dimensional timeseries data includes a large number of noise features, the density estimation accuracy tends to be degraded by noise.

To avoid using density estimation, direct density-ratio based change-point detection approaches have been proposed [Kawahara and Sugiyama, 2009; Liu et al., 2013]. These approaches estimate the ratio of probability distributions directly without using density estimation. Direct density-ratio estimation has been actively studied by the machine learning community and techniques include kernel mean matching (KMM) [Huang et al., 2007], the Kullback-Leibler importance estimation procedure (KLIEP) [Sugiyama et al., 2008], WKV [Nguyen et al., 2010], and unconstrained leastsquares importance fitting (uLSIF) and its robust extension called relative uLSIF (RuLSIF) [Kanamori et al., 2009; Yamada et al., 2011]. These direct density-ratio estimation methods have exhibited the optimal convergence rate for nonparametric density-ratio estimation. However, as with density estimation based methods, the accuracy of the density-ratio estimation is likely to be degraded by noise features.

Change-point detection with stationary subspace analysis (SSA), which is a dimensionality reduction method, is a promising change-point detection method for multivariate time-series data [Blythe et al., 2012]. SSA factorizes a multivariate time-series data into stationary and non-stationary sources, and the change-points can be detected in a nonstationary subspace. Since SSA can reduce the dimensionality of data without losing the abrupt change characteristic, it can significantly improve change-point detection performance. However, since SSA needs to compute the log of a covariance matrix which is singular when (d < n), it needs a large number of training samples (n ≫ d) to accurately factorize stationary and non-stationary sources. Therefore, an SSA based change-point detection algorithm is not applicable to high-dimensional change-point detection problems.

In this paper to deal with high-dimensional time-series data, we integrate feature selection into change-point detection measure estimation. More specifically, we present a supervised learning based change-point detection approach in which we use the separability of past and future data at time t (they are labeled as +1 and -1) as a plausibility of a changepoints [Desobry et al., 2005; Hido et al., 2008]. If we can separate the past and future data sets easily, we regard them as having different probability distributions. On the other hand, if two data sets are not separable, we regard them as having the same probability distribution. Using this framework as a basis, we propose a detection measure called the additive Hilbert-Schmidt independence criterion (aHSIC), which is given as the weighted sum of the HSIC scores where HSIC is a kernel-based independence measure [Gretton et al., 2005]. Here, each HSIC score is computed from samples of a feature and its corresponding pseudo binary labels. An advantage of the aHSIC score over existing detection measures is that it can incorporate feature selection during detection measure estimation. That is, we first select features that are responsible for an abrupt change in a supervised manner, and then compute an aHSIC score by using those selected features. Thus, compared with traditional detection measures, aHSIC is more robust as regards noise features than existing measures and it is suited for high-dimensional time-series data. Experiments on synthetic and real-world human activity data showed that the proposed methods are promising.

2 ProblemFormulation

In this section, we formulate our change-point detection problem based on supervised learning framework [Desobry et al., 2005; Hido et al., 2008].

Let x(t) 2 Rd be a d-dimensional sample at time t and

are samples with length n extracted from time series data in a sliding-windowmanner at time t.

Let us consider two non-overlapped sequences X(t) and X(t + n), and we annotate the samples in X(t) as y = 1 and the samples in X (t+n) as y = -1. We denote the augmented

where yi 2 ¶+1↪ -1♦ is a pseudo binary label. See Figure 1 for an illustrative example of two intervals X(t) and X(t +

n). Note, in high-dimensional time-series data, the number of dimensions d is usually bigger than that of samples 2n (i.e., d > 2n).

Then, our change-point detection strategy is to compute a dependency score between input x and output y from data Z(t), and use it as the change-point plausibility:

where D(Z(t)) is a dependency measure that takes a large value if x and y are dependent and ⌧ is a threshold that controls the sensibility/robustness tradeoff. Note, since y can be regarded as a step function like sequence that exhibits a change at time t, a large dependency value means that X(t) and X(t + 1) are separable. That is, we can regard X(t) and X(t +n) as samples from different distributions. In contrast, if the dependency score is small, then we can regard X(t) and X(t + n) as samples from the same distribution.

3 Proposed Change-Point DetectionMethod

The detection performance depends strongly on the dependency measure. In particular, since high-dimensional data tends to include a large number of noise features, a key issue is to compute a detection measure for an abrupt change using only important features. Thus, we incorporate feature selection into our change-point detection measure. To the best of our knowledge, this is the first work to perform feature selection in change-point detection for high-dimensional data.

In this section, we first propose our dependency measure, and then show a way to estimate it.

3.1 Additive Hilbert-Schmidt Independence Criterion

We propose an additive HSIC (aHSIC) score as a dependency measure D(Z(t)):

where uk(t) = [xk(t - n + 1)↪ xk(t - n + 2)↪ ▷ ▷ ▷ ↪ xk(t + n - 1)]> 2 R2n is the k-th feature for all samples, > is the matrix transpose, y = [1↪ ▷ ▷ ▷ ↪ 1↪-1↪ ▷ ▷ ▷ ↪-1]> 2 R2n is the pseudo binary label vector, HSIC(uk↪ y) = tr( ¯K(k) ¯L) is a kernel-based independence measure called the (empirical) Hilbert-Schmidt independence criterion (HSIC) [Gretton et al., 2005], tr(·) denotes the trace, ↵1↪ ▷ ▷ ▷ ↪ ↵d ≥ 0, Pk=1 ↵k = 1, ¯K(k) = ΓK(k)Γ and ¯L = ΓLΓ are centered and normalized Gram matrices, K(k) i↪j = K(xk↪i↪ xk↪j) and Li↪j = L(yi↪ yj) are Gram matrices, K(x↪ x0 ) and L(y↪ y0 ) are kernel functions, Γ = In - 1n 1n1>n is the centering matrix, In is the n-dimensional identity matrix, and 1n is the n-dimensional vector with all ones.

Note, HSIC always takes a non-negative value, and is zero if and only if two random variables are statistically independent when a universal reproducing kernel [Steinwart, 2001] such as a Gaussian kernel is used. That is, if the k-th feature uk is independent of y (i.e., the k-th feature is not important for an abrupt change), HSIC(uk↪ y) takes a small value.

A novelty of the aHSIC score is that it is possible to measure dependency based solely on features that are related to an output y if we set ↵ appropriately. That is, if we can select features that are responsible for abrupt changes, the aHSIC score is independent of noise features.

3.2 HSIC Lasso

In aHSIC, the choice of ↵ parameter is a key issue. A simple heuristic is to set 1n for all ↵. However, for high-dimensional time series data, a few features are important and rest are noise. That is, equal weighting of the HSIC scores is not a suitable choice for high-dimensional time-series.

In this paper, we propose using HSIC Lasso for estimating the ↵ parameter [Yamada et al., 2012]:

where k · kFrob is the Frobenius norm and k · k1 is the `1 norm. The first term in Eq.(1) means that we are regressing the output kernel matrix ¯L by a linear combination of featurewise input kernel matrices ¶ ¯K(k) ♦dk=1 . After estimating ↵, we normalize each element of ↵ as ↵k ↵k◁Pdk=1 ↵k. The first term in Eq.(1) can be rewritten as

Thus, if all features uk are mutually independent (i.e., HSIC(uk↪ ul) = 0↪ 8 k↪ l), we can rewrite Eq.(2) as

This means that, minimizing the objective function of HSIC Lasso corresponds to maximizing the aHSIC(Z) score. Thus, using HSIC Lasso to estimate ↵ is a natural choice.

Statistical Interpretation of HSIC Lasso: From Eq.(2), if the k-th feature uk has high dependence on output y, HSIC(uk↪ y) takes a large value and thus ↵k should also be large. On the other hand, if uk and y are independent, HSIC(uk↪ y) is close to zero and thus such ↵k tends to be removed by the `1-regularizer. That is, relevant features that have strong dependence on output y tends to be selected by the HSIC Lasso. That is, features that are important for an abrupt change are selected.

Furthermore, if uk and ul are strongly dependent (i.e., redundant features), HSIC(uk↪ul) is large and thus either ↵k or ↵l tends to be zero. Thus, redundant features tend to be removed by HSIC Lasso.

Overall, HSIC Lasso tends to find non-redundant features with strong dependence on output y, which is a preferable property for a change-point detection measure.

Computational Property: An important computational property of HSIC Lasso is that the first term in Eq.(1) can be rewritten as

where vec(·) is the vectorization operator. This is the same form as plain Lasso with n2 and d are the numbers of samples and optimization parameters, respectively.

To solve this Lasso problem, a technique called the dual augmented Lagrangian (DAL) was shown to be computationally highly efficient [Tomioka et al., 2011]. Because DAL can also incorporate the non-negativity constraint without losing its computational advantages, we can directly use DAL to solve our HSIC Lasso problem.

A MATLAB R© implementation of the HSIC Lasso is available from http://www.kecl.ntt.co.jp/icl/ls/members/ myamada/hsiclasso.html.

3.3 Group HSIC Lasso

If we have a prior knowledge of the features, we can utilize this knowledge to select features in HSIC Lasso. For example when detecting a change in human activity from sensors attached to the hands and legs, it is reasonable to select a group of features (i.e., sensors) that are important in relation to an abrupt change.

To this end, we propose using the group-type lasso regularizer [Meier et al., 2008; Zou and Hastie, 2005] for estimating ↵ as

where ↵ = [↵>1 ↪ ▷ ▷ ▷ ↪ ↵>G]>, ↵g is the gth group of variables, and G is the number of groups. This group-type lasso problem can also be efficiently solved by the DAL package with the non-negativity constraint.

3.4 Kernel Selection

HSIC is a kernel based independence measure; the independence criterion changes with respect to a kernel parameter and/or kernel types. That is, the input kernel parameters should be fixed for all features. Thus, we first normalize the input features with a standard deviation at 1 and then use the same kernel parameters for all kernels.

For input x, we use the Gaussian kernel,

where we set σx = 1. For output y, we use the delta kernel,

4 RelatedMethods

Here, we review related change-point detection measures.

4.1 KLIEP

Let us assume that samples in X(t) are drawn i.i.d. from a distribution with a density p(x) and samples in X (t + n) are drawn i.i.d. from a distribution with density p0 (x). Then, we can use divergence as the plausibility of the change-points. More specifically, if the divergence D(p(x)♣♣p0 (x)) is larger than the threshold ⌧, we can regard the point as the changepoint.

A popular choice for the divergence function is KullbackLeibler (KL) divergence:

A naive approach for estimating the KL-divergence is to estimate the probability densities p(x) and p0 (x) separately using a kernel density estimation [H¨ardle et al., 2004] and then take the ratio. However, density estimation is known to be a hard problem, and the KL-divergence estimation can be poor. To mitigate this, an alternative KL divergence estimation method is proposed, where the density-ratio p(x) p0 (x) is directly estimated without going through the density estimations by using the Kullback-Leibler Importance Estimation Procedure (KLIEP), which has been proved to achieve the optimal non-parametric convergence rate in a mini-max sense [Sugiyama et al., 2008].

In [Kawahara and Sugiyama, 2009], a KLIEP based KLdivergence estimator was used for change-point detection, and it outperformed existing density estimation based methods. However, since KLIEP uses all the features to compute the divergence, the KL-divergence estimator tends to perform poorly when there are many noise features.

4.2 uLSIF/RuLSIF

Recently, computationally efficient density-ratio estimation method called relative unconstrained Least-Squares Importance Fitting (RuLSIF) has been proposed for estimating relative Pearson-divergence[Yamada et al., 2011]:

where q↵(x) = ↵p(x) + (1 - ↵)p0 (x) is called ↵-mixed density and r↵ is the relative density ratio:

The novelty of the relative density-ratio is that it is always bounded above by 1↵ , and it has been shown that the convergence rate for estimating relative density ratio is faster than that of the standard density-ratio p(x) p0 (x) [Yamada et al., 2011]. This is a favorable property in practice. Note, when ↵ = 0, it is reduced to p(x) p0 (x) (uLSIF) [Kanamori et al., 2009].

It has been shown that the relative PE divergence estimator is suited to change-point detection [Liu et al., 2013]. However, as with the KL-divergence estimator, estimation of relative PE divergence tends to be poor when there are many noise features and a small number of training samples.

4.3 Kernel Change Detection

Kernel change detection (KCD) is also a detection method that does not require density estimations [Desobry et al., 2005].

Let βt↪1 and βt↪2 be coefficients of a one-class support vector machine (OSVM) computed from X(t) and X(t + 1) respectively, Kt↪12 is the Gram matrix computed from support vectors obtained from X (t) and X (t + 1), Kt↪11 is the Gram matrix computed from support vectors obtained from X(t), and Kt↪22 is the Gram matrix computed from support vectors obtained from X (t+1). Then, the dissimilarity measure used in KCD is given as

This dissimilarity measure is shown to be asymptotically equivalent to the Fisher ratio in the Gaussian case.

KCD is robust to outliers, since outliers are removed by OSVM. However, since KCD uses all the features for detecting change-points, it can perform poorly when there are many noisy features.

4.4 Change-Point Detection using SSA

Change-point detection using stationary subspace analysis (SSA) is a promising method for multivariate time-series data [Blythe et al., 2012].

Let us divide the entire time series X into X1↪ ▷ ▷ ▷ ↪ XN sets and define the hypothesis of testing non-stationarity as

where N(θ↪ ⌃) is the Gaussian distribution with mean θ and variance ⌃.

Then, the goal of the SSA based change-point detection is to test whether time-series data are stationary. More specifically, if we can reject the null hypothesis H0, we can consider that there is a change-point. On the other hand, if we cannot reject the null hypothesis H0, we consider there to be no change-point.

The test statistics based on SSA can be written as

where bθ bBi = bBbθi and b⌃ bBi = bBb⌃i bB>, bθi = 1 ♣Xi♣ Pt2Xi x(t), b⌃i = Pt2Xi (x(t)- bθi)(x(t)- bθi)>, and bB 2 Rm⇥d (m < d) is given as the solution of the following optimization problem:

It has been shown that SSA-based change-point detection performs well on multivariate time-series data. However, the SSA based change-point detection algorithm needs a large number of training samples (i.e., n ≫ d), since it includes the log of the covariance matrix which is singular when (d < n). In addition, SSA needs several non-overlapped time-series. Thus, the SSA-based change-point detection algorithm is not applicable to high-dimensional change-point detection problems.

5 Experiments

In this section, we investigate experimentally the performance of the proposed and existing feature selection methods using synthetic and real-world human activity datasets.

5.1 Setup

We compare the performance of the proposed methods with that of RuLSIF, KLIEP, and KCD. With RuLSIF and KLIEP, we use publicly available codes. For KCD, we use LIBSVM [Chang and Lin, 2011] to compute OSVM, where we use a Gaussian kernel with σ = pd◁2 and the regularization parameter of OSVM ⌫ = 0▷5. We experimentally fix λ at 0▷01 for HSIC Lasso and λ at 10▷0 for Group HSIC Lasso. For all the methods, we fix the window size at 2n = 40.

In this paper, we compare the performance of change-point detection methods objectively in terms of the receiver operating characteristic (ROC) curves and the area under the ROC curve (AUC) values. Note that, detection at t is regarded as correct if there exists a true alarm at step t⇤ such that t 2 [t⇤ - 10↪ t⇤ + 10].

5.2 Synthetic Datasets

In this section, we illustrate the behavior of the proposed additive HSIC based change-point detection method using synthetic datasets. We generate time-series data so that one feature includes abrupt changes and the remaining features are noise.

We use the following two synthetic multivariate time-series datasets, which contain manually inserted change-points:

• (a) Data1 (Jumping mean): The following 1-dimensional auto-regressive model borrowed from [Yamanishi and Takeuchi, 2002] is used to generate 1000 samples

where ✏t ⇠ N(0↪ 1). The initial values are set as x1(1) = x1(2) = 0. A change-point is inserted at every 100 times steps by setting the noise mean θ at time t as

where M is a change-point index such that 100(M - 1) + 1  t  100M.

Then, we generate a noise vector (x2(t)↪ ▷ ▷ ▷ ↪ x50(t))> ⇠ N(049↪ I49) and concatenate it to x1(t) as x(t) = [x1(t)↪ x2(t)↪ ▷▷▷↪ x50(t)]>.

• (b) Data2 (Scaling variance): The same auto-regressive model is used as with Data1, but a change-point is inserted at every 100 time steps by setting the noise standard deviation σ at time t as

Then, we generate a noise vector (x2(t)↪ ▷ ▷ ▷ ↪ x50(t))> ⇠ N(049↪ I49) and concatenate it to x1(t) as x(t) = [x1(t)↪ x2(t)↪ ▷▷▷↪ x50(t)]>.

Figure 2-(a),(c),(e),(g) show part of x1(t) and the corresponding aHSIC score for each data set. The vertical dotted red lines in those figures denote the true change-points. As can be seen, the proposed approach can correctly detect change-points in both cases. Figure 2-(b),(d) show the estimated ↵1 values. It is clear that the true feature is successfully selected by HSIC Lasso. Figure 2-(f),(h) show the ROC curves. The experimental results show that proposed method compares favorably with existing methods.

Table 1 shows the mean computational time of the proposed method over the Data1 data set. As can be observed, the computational time of aHSIC score is reasonable. Note that OSVM is implemented with C, while the proposed method is implemented with Matlab. Thus, by implementing the proposed method with C/C++, we can boost the computational speed of the proposed method.

5.3 Real-World Human Activity Dataset

The proposed feature selection based approach is very useful for human activity change detection problems. For example when detecting the change from "Standing" to "Brushing teeth", it is reasonable to only use right/left hand information for detection.

In this section, we report the change-point detection of human activity data set [Maekawa and Watanabe, 2011]. The data set contains 14 actions from 61 subjects, where each subject wore three-axis acceleration sensors on the both hands, the waist, and the right thigh. Then, we computed the mean, energy, entropy, and main frequency component (F0) for each axis (12 features for a sensor and 48 features in total). For Group HSIC Lasso, we use four groups namely the right/left hands, waist, and right thigh.

In this experiment, we first randomly concatenate all the motion sequences and check whether the change-point is correctly detected. Figure 3 shows the ROC curves (averaged over 10 runs). The AUC values of the aHSIC, aHSIC (Group), RuLSIF, KLIEP, and KCD were 0.896, 0.907, 0.853, 0.820, and 0.869, respectively. Moreover, paired t-tests were conducted, and we observe that aHSIC and its group version outperform existing methods at p = 0▷01(1%). The experimental results show that the proposed method compares favorably with existing methods.

6 Conclusion

In this paper, we proposed a change-point detection method with feature selection for high-dimensional time-series data. We adopted the supervised change-point detection approach [Hido et al., 2008] in which we use the separability of the past and current data sets (they are labeled +1 and -1, respectively) as the change-point detection measure. Based on this framework, we proposed a new change-point detection measure called the additive Hilbert-Schmidt Independence Criterion (aHSIC), which is defined as the weighted sum of HSIC values between each feature and its corresponding pseudo binary labels. An advantage of the proposed method over existing methods is that it is estimated by using features that are important for an abrupt change. That is, the proposed approach is more robust to noisy features than existing methods. Through extensive experiments on synthetic and real-world human-activity dataset, we demonstrated the promise of the proposed change-point detection method.

Following the current line of research, there are several issues to be pursued if we are to further improve the change point detection performance. For example, in this paper, we assume that separable features are useful for detecting an abrupt change. Thus, if there is a feature that is independent of an abrupt change and has high separability, false positives of the proposed method can be increased. We have already verified experimentally that the proposed method performs well, but the issue definitely constitutes interesting future work. Moreover, investigating the performance of the proposed method over other real datasets such as music change detection is also a challenging future work.

Acknowledgments

The authors thank Dr. Katsuhiko Ishiguro and Mr. Koh Takeuchi for their valuable comments.