Joseph Ndong
email: joseph.ndong@lip6.fr

A Robust Anomaly Detection Technique using Combined Statistical Methods

Keywords: Anomaly Detection, System Monitors, Kalman filter, GMM, HMM

Parametric anomaly detection is generally a three steps process where, in the first step a model of normal behavior is calibrated and thereafter, the obtained model is used in order to reduce the entropy of the observation. The second step generates an innovation process that is used in the third step to make a decision on the existence or not of an anomaly in the observed data. Under favorable conditions the innovation process is expected to be a Gaussian white noise. However, in practice, this is hardly the case as frequently the observed signals are not gaussian themselves. Moreover long range dependencies, as well as heavy tail in the observation can lead to important deviation from the normality and the independence in the innovation processes. This, results in the frequent observation that the decisions made assuming that the innovation process is a white and Gaussian results in a large false positive rate. In this paper we deal with the above issue. Our approach consists of not assuming anymore that the innovation process is Gaussian and white. In place we are assuming that the real distribution of the process is a mixture of Gaussian and that there are some time dependency in the innovation that we will capture by using a Hidden Markov Model. We therefore derive a new decision process and we show that this approach results into an important decrease of false alarm rates. We validate this approach over realistic traces.

Introduction

In the Internet, anomalous traffic behavior such as attacks, configuration changes, flash crowds, large transfer file and outages occur frequently. Large enterprise networks often have a security operations center where operators continuously monitor the network traffic hoping to detect, identify and treat anomalous events. The detection process of these events can then be used to trigger alarms to the network management system, which, in turn, trigger recovery mechanisms. Despite the recent growth in monitoring technology and in intrusion detection systems, correctly detecting anomalies in a timely fashion remains a challenging task. One of the reasons for this is that, many today's security solutions yield equipments that collect and analyze traffic from one link at a time. Similarly many research efforts consider anomaly detection on a per link basis [START_REF] Hussain | Measurement and Spectral Analysis of Denial of Service Attacks[END_REF], [START_REF] Barford | A signal analysis of network traffic anomalies[END_REF], [START_REF] Barford | Characteristics of network flow anomalies[END_REF].

To detect traffic anomalies, one typically seeks to characterize or build a model of normal behavior. In this paper we use signal processing techniques to track anomalies in situations where normal network conditions occur. The approaches used to address the anomaly detection problem are dependent on the nature of the data that is available for analysis. Network data can be obtained at multiple levels of granularity such as end-user-based (TCP/UDP data) or network-based. Traffic counts obtained from both types of data can be used to generate a time series to which statistical signal processing techniques can be applied [START_REF] Thottan | Fault detection in ip networks[END_REF], [START_REF] Wang | Detecting syn flooding attacks[END_REF]. In this work, we focus our attention on network-based (at the IP flow level) anomaly detection. After filtering out normal looking traffic, anomaly detection methods analyze the residual traffic pattern for deviations. Taking account only one link per unit time is limiting, since any flow will traverse multiple links along its path. And it is intuitive that, a flow carrying an anomalous event will appear in these links, thus increasing the evidence to detect it.

In this paper we will show that, it is more convenient to track anomalies by running the detection mechanism in all links at the same time, than in the case of per-link detection. To do this, we focus on using data from all links in an enterprise or ISP network simultaneously. Since any anomaly has to traverse several links on route to its destination, an anomaly has the potential to be visible in any of the links it traverses. Since we cannot know in advance where anomalies will originate, nor the path they will take, it is advantageous to consider the behavior of all links in the same time, when developing both a model of "normal" traffic and a method for analyzing "residuals".

Our study is built on traffic matrix scheme [START_REF] Soule | Traffic Matrix Tracking using Kalman Filters[END_REF], [START_REF] Soule | Combining Filtering and Statistical Methods for Anomaly Detection[END_REF]. Each traffic matrix entries describes the average volume of traffic, in a given time interval, that originates at a given source node (routers) and is headed towards a particular destination node. In this paper we propose to use predictions of traffic matrix behavior for the purpose of anomaly detection. Since the traffic matrix is a representation of traffic volume, the types of anomalies we might be able to detect via analysis of the traffic matrix are volume anomalies, [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF].

In [START_REF] Soule | Combining Filtering and Statistical Methods for Anomaly Detection[END_REF], we used traffic matrix at the granularity of OD flows to detect, identify and track anomalies, by comparing four schemes (a basic analysis using variance, the CUSUM and Generalized Likelihood Ration test combined, Wavelets and a fourth method based multi scale variance shift). We had shown that the basic analysis using variance performs the best above the four methods, in the case of real network data (Abilene). In this work, we show that our method perform better that that analysis using variance, with the same network data. Our study differs, to the previously cited work, in numerous ways: i) instead of using OD flows, we used granularity of link counts ii)we do not assume that the Kalman innovation process does strictly remain a zero mean gaussian process,iii) instead, we assume that the real distribution of the process is a mixture of gaussians, making sense the use of gaussian mixture modeller to show that residuals can be split in different families (where anomalies might or might not happen), iv) we use also hidden markov model in order to capture time dependencies over the different mixture components, v) we use the Viterbi algorithm to select the potential best sates path where anomalies might come from v) we use finally variance shift analysis to track the anomalies.

To start our study, first a linear dynamical system is needed to model the link counts. Thus, a model is built, where the Kalman filter algorithm is derived to estimate the state of the system, for denoising the observed data. Thereafter, we use a gaussian mixture modeller above the innovation process (obtained from the Kalman filter phase), for clustering. This makes our monitoring technique robust toward gaussianity, since we do not assume that the residual remains typically a zero mean gaussian process. These two steps constitute our learning phase. In the second step, we perform an operational phase where, first, we build a discrete sequence of finite alphabet of mixture memberships IDs (based on the residual clustered in the GMM-learning phase), using a maximum a posteriori (MAP) criteria. This discrete sequence is then used to check the sequence of hidden states (each state containing a part of our alphabet) of the HMM, using again a MAP criteria. We perform anomaly detection on this results, combining finally the use of the Viterbi algorithm (in order to find the only best sequence of our hidden states) and a variance shift test.

The organization of this paper is as follows. Section 2 describes the system monitoring design we used. Section 3 deals with the methodology we adopt in our anomaly detection scheme. In section 4, we detail our calibration method and we validate our approach by showing efficient results. Section 5 concludes our work and fixes some ideas for future works.

Description of the Monitoring System

In our paper referred in [START_REF] Soule | Traffic Matrix Tracking using Kalman Filters[END_REF], a simple basic system was designed to monitor OD flows traffic coming from network backbone. This system was produced with three functional blocks: data collection (by a NOC-Network Operations Center), a data analysis block and a decision process phase. The data analysis center is simply an operational phase built on Kalman filtering. We extend this anomaly detection system, to take account to three new components shown as part of our learning phase. In addition, our decision process is based on a combination of two underlying materials: the Viterbi algorithm and a variance test based data-point analysis. In Figure 2, the data analysis block contains a gaussian mixture modeller, a hidden markov modeller and a maximum a posteriori criteria used twice. Recall that we have made no assumptions that the residual process (on which we perform our anomaly detection scheme) is a zero mean gaussian process; instead the kalman innovation is supposed to be driven by a family of normal distributions. One can think that, this idea can be a good decision, if we know that the real distribution of the original data themselves is not purely gaussian. Based on this belief, we have introduced the use of a GMM, which is a good paradigm for modelling processes in situations where data might come from different families of exponential distributions. So, the calibration of a gaussian mixture model will help us for data clustering, then making possible to classify each residual data point into a gaussian component. After building a finite set of classes of residuals, all elements belonging to the same class might appear with the same identifier. The motivation to do so, is based on the fact that, one should easily identify a cluster (with label 1,2,3,...) and at the other hand, as we will see later, the use of a hidden markov model will be adequate to study the temporal dependencies between all the clusters. To perform this task, the gaussian mixture model is followed by a phase where we run, using a Maximum A Posteriori (MAP)criteria, an operation making us possible to form a discrete finite alphabet (sequence of symbols), using all the GMM components. At this step, all the continuous observation sequences are transformed into a discrete observation sequence, which can be easily re-used to perform a monitoring method for the purpose of anomaly detection.

The GMM phase for clustering and the MAP criteria operation for discrimination, are followed by the use of an hid- den markov model, for time dependencies tracking, in the innovation process. One can think that time dependencies might involve in the innovation process, and these temporal correlation might appear between the different clusters we have already built. This means that the different families of discrete alphabet can be re-organized into different other groups, each of them containing a subset of the GMM components. Since these groups are a priori unknown, we denote them as hidden states. The use of hidden markov model will be a good paradigm to achieve this aim. Here, we perform a maximum a posteriori criteria to obtain the sequence of the hidden states of our HMM. Finally, the Figure 2 shows that, to track anomalies buried in the link counts, the use of a decision process is necessary, for management and decisions issue. In our anomaly detection philosophy, since anomalous events might be rare, we believe that it is possible to find one typically unique path, where these anomalies might be identified and detected. So, to perform the management task, we propose in this work a decision process based on two routines: the use of the Viterbi algorithm, in order to find all the variations in the residual process and, as a final step, the use of variance test for each data-point in the residual. Clearly, we will build a GMM with N components, then we use these N components to form an HMM with P states(each state containing part of the N components), and at last the Viterbi algorithm will find one unique path with a sequence of K(¡P) states and finally we use these K states to perform anomaly detection. Based on the discrete innovation sequence, the Viterbi algorithm finds the only best sequence that capture all the variations in the continuous form of the innovation process. In the case of per-link analysis, we build a Viterbi path for each data observation sequence, but in the case where we analyze all the links simultaneously, we find one single best Viterbi sequence. Our study shows thereafter that, a simple basic variance test based point-analysis suffices to identify and track the anomalous events.

Methodology

The first step in our study is to seek for a model which can monitor the system traffic dynamics. We believe that observing and estimated the real values of these metrics can help networks administrators to understand normal or potential abnormal conditions in their environment. A powerful model might identify the behaviour of normal traffic, by capturing the temporal evolution of the data, and at least can help operators to do future prediction. The full model is represented by the two linear equations, combined to form the complete specification of our system:

X t+1 =C t X t + W t Y t = A t X t + V t (1)
This block equations is the classic form for linear dynamical system (without inputs). In this model we assume that the state-noise W t and the measurement-noise V t are uncorrelated zero-mean gaussian white-noise processes with covariance matrices Q t and R t , respectively. For a full understanding of these system, we refer the reader to our work in [START_REF] Soule | Combining Filtering and Statistical Methods for Anomaly Detection[END_REF]. Since we have found a model and system equations for our system, next we need to deal properly with the different steps of our optimization algorithm for anomaly decision issue.

Equations of the Kalman filter

The first problem to solve after building a simple model to monitor link counts, is to find an optimal estimate ( Xt ) of our unobservable network states X t , given a set of measurements {Y 1 , ....., Y t }. A well-known and robust method to achieve this aim is the Kalman filtering algorithm. In our dynamical linear system, we refer to Y t as the observation vector at a specific time t. And the state of the system at time t is given by X t , let also Xt|k denotes the estimate of X t using all the information available up to time k, i.e, ∀τ < k. Xt+1 denotes the estimate of X t+1 using all the information up to time t, (this constitutes the phase predictor). The quantity Xt+1|t+1 denotes the estimate of X t+1 using all past information and the recently arrived data point at time t+1. In the other hand, P t|t denotes the covariance of the state estimate and P t+1|t indicates the covariance of the state prediction. As it is shown in its earlier elaboration, the Kalman filter addresses the problem of estimating a discrete state vector when the observations are only a linear combination of the underlying state vector. The filter runs as a predictor-corrector algorithm. As an iterative algorithm, it estimates the system state using two steps: prediction comes in the time update phase, and correction in the measurement update phase.

• Prediction step (time update equations):

In this step, the estimated state of the system at time t, Xt|t , is used to predict the state at next time t+1, Xt+1|t . And, as we know that the noise W t influences the evolution of the system at each time t, we compute only the covariance of the prediction, P t+1|t based on the updated covariance at the previous time t, P t|t , and the noise covariance at the same time, Q. The error covariance P t+1|t provides an indication of the uncertainty associated with the state estimate.

Xt+1|t = Xt|t P t+1|t =P t|t + Q (2)
• Correction step (measurement update equations):

This step updates (corrects) the state and the variance of the estimate in the previous step, using a combination of their predicted values and the new observations Y t+1 . The correctness of this update depends on the Kalman innovation Y t+1 -Xt+1|t .

Xt+1|t+1 = Xt+1|t + K t+1 (Y t+1 -Xt+1|t ), P t+1 =(I -K t+1 )P t+1|t (I -K t+1 ) T + K t+1 RK T t+1 (3)
In the measurement equations, K t+1 denotes the Kalman gain. For more details in linear dynamical system, estimation and Kalman filter techniques, we refer the reader to: [START_REF] Kalman | New results in linear filtering and predictions[END_REF], [START_REF] Kailath | Linear Estimation[END_REF], [START_REF] Wolverton | On the Linear Smoothing Problem[END_REF] and [START_REF] Raugh | Solutions to the linear smoothing problem[END_REF]. The above equations with initial conditions of the state of the system X0|0 = E[X 0 ] and the associated error covariance matrix

P 0|0 = E[( X0 - X 0 )( X0 -X 0 ) T ]
define the discrete-time sequential recursive algorithm for determining the linear minimum variance estimate known as the Kalman filter.

Gaussian Mixture Model

In Kalman filtering philosophy it is more generally assumed that the residual remains a zero mean gaussian process, however this assertion is not always true in practice. Long range dependences and heavy tails distribution can lead to non negligible deviation from the normality and the independence of the innovation process. Thus our motivation to use gaussian mixture model is based on our belief that the real distribution of the process is an ensemble (mixture) of gaussians and there is some time dependency in the innovation (which we will later study with the aid of an hidden markov model). This assertion allows us to build a method which has the ability to find anomalies in different families built on the same residual sequence. The calibration of a GMM helps us, taking as input the Kalman innovation, to build K families (clusters) and in the next step we propose the use of an hidden markov model (HMM) to classify these groups into P states, each state being formed with part of the K families. This operation have the advantage to discover the potential temporal dependencies in the innovation process.

To find the values of the GMM model parameters µ m (mean)and Σ m (variance), as well as the prior probability vector π, we are interested in maximizing the likelihood £(θ|X) = p(X; θ) of generating the known observed data (X) given the model parameters θ = {µ m , Σ m , π m }, 1≤m≤M. X denotes all the observation while θ contains all the parameters of the mixture. In other words, we hope to find θML =argmax p(x|θ). This approach is called the Maximum Likelihood (ML) framework since it finds the parameter settings that maximize the likelihood of observing the data sets. To find the best parameters of the features of θ, the Expectation Maximization (EM) iterative algorithm can be used to simplify the math considerably and numerically compute the unknown parameters. For more details about mathematical routines for EM, see [START_REF] Bilmes | A Gentle Tutorial of the EM algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models[END_REF], [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF], [START_REF] Dempster | Maximum likelihood from in-complete data via the em algorithm[END_REF].

Hidden Markov Model

From the above learning phase, each GMM component (cluster) is transformed into a sequence of a finite set of alphabet (symbols as 1,2,3,...), using a maximum a posteriori criteria. This discrimination phase will help for plugging the above clusters into different a priori unknown states, using hidden markov model. Our results will show us that, clusters with anomalies appear on some states and clusters with no anomalies remain in an other states. We represent these families of states by the following collection of unknown random variables {Q 1 , Q 2 , .....Q T } (where Q t is a constant value with values in {1, 2, ..., K}). We also represent our alphabet by the known vector {O 1 , O 2 , ....., O N }. Now, the problem is resumed to find a model to produce the states and to determine the probability of each symbol being in a state. A well known model-based approach to tackle this problem, is the discrete hidden markov model (HMM) which is an exponential family of mixture model. Our choice of using HMM is based on the fact that: i) potential time dependencies in the innovation process can be modelled and captured using a finite set of a priori hidden states, each of them containing a subset of gaussian components ii) relatively efficient algorithm can be derived to solve the problems related to them, [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF], [START_REF] Bilmes | A Gentle Tutorial of the EM algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models[END_REF], [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF]. The full HMM model we used is defined by the quantity λ = (A, B, π) (where A is the transition matrix, B is the emission probabilities matrix and pi the prior probabilities). To find and estimate the best parameters of our model, we use the well-known forward-backward algorithm parameter estimation (or Baum-Welch algorithm). For more details for EM techniques related to hidden markov model, see, [START_REF] Bilmes | A Gentle Tutorial of the EM algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models[END_REF], [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF]. Thereafter, we reuse the model to find the optimal state sequence associated with the given observation sequence. We believe that this final step of our approach will allow us to capture all the variations in the innovation pro-cess. An optimal criteria we have chosen here is to find the single best state sequence (path), Q = {q 1 , q 2 , ..., q T } for the given observation sequence O = {O 1 , O 2 , ..., O T },i.e., we aim to maximize P (Q|O, λ). A formal technique for finding this unique best state sequence is the Viterbi algorithm.

Model Evaluation

Data collection : Abilene network

The Abilene backbone has 11 Points of Presence(PoP) and spans the continental US. The data from this network was collected from every PoP at the granularity of IP level flows. The Abilene backbone is composed of Juniper routers whose traffic sampling feature was enabled. Of all the packets entering a router, 1% are sampled at random. Sampled packets are aggregated at the 5-tuple IP-flow level and aggregated into intervals of 10 minute bins. This thus dictates the underlying time unit of all of our estimations and detections. The raw IP flow level data is converted into a PoP-to-PoP level matrix using the procedure described in [START_REF] Lakhina | Characterization of network-wide traffic anomalies[END_REF]. Since the Abilene backbone has 11 PoPs, this yields a traffic matrix with 121 OD flows. Note that each traffic matrix element corresponds to a single OD flow, however, for each OD flow we have a seven week long time series depicting the evolution (in 10 minute bin increments) of that flow over the measurement period. All the OD flows have traversed 41 links.

Model Validation

The anomalies injected in the Abilene data are small and high synthetic volume anomalies. We used exactly the same Abilene data as in [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF]. So for a full understanding on how the ground-truth is obtained (based on EWMA and Fourier algorithms) , we refer the reader to [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF].

System parameter identification.

Our method begins with a learning phase where we calibrate a Kalman filter for denoising, using our linear dynamical system, and the collection of observation data. In order to run the Kalman filter, we need the matrices A, Q, C et R. The matrix A is available given the routing scheme of a network. We thus only need to obtain Q, C and R. Recall that one of our aims is to prove that, learning all the links simultaneously is better than using one link at a time. Calibrating the system for one link at a time is not very straightforward. Here to estimate the parameters θ={C, Q, R}, we use an autoregressive (AR) model of order p. The method is based on a stepwise least squares algorithm which uses a QR factorization of a data matrix to evaluate, for a sequence of successive orders, a criterion (here Schwarz's Bayesian Criterion and Akaike's Final Prediction Error) for the selection of the model order p, and to compute the parameters of the AR model of the optimum order. It is sufficient to learn the model's parameters using a sample of one week measurements from the data links studied. We run this method based on the notes described in [START_REF] Neumaier | Estimation of parameters and eigenmodes of multivariate autoregressive models[END_REF], and developed in [START_REF] Schneider | Algorithm 808: ARfit -A Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models[END_REF]. We form the matrix C using the p parameters of the AR(p). The noise variance of the model is used to fix Q. The matrix R is built as a function of Q. In practice, R is obtained by dividing the measurement noise by some constant. In our experiments we have found this constant in the interval [2:4]. In the case where we take account of all links at the same time, we find the same parameters θ, using N consecutive samples of link counts, and we do our calibration using the EM algorithm. In practice, in order to apply the EM procedure, initial values are required for the parameters and these were simply guessed by examining portions of the completely observed series. It is a convenient idea to examine several different sets of starting values, since the EM algorithm may reach different kinds of stationary values corresponding to local rather than global maxima. For more details of the EM techniques, see [START_REF] Shumway | An Approach to Time Series Smoothing And Forecasting Using the EM Algorithm[END_REF], [START_REF] Shumway | Dynamic Linear Model With Switching[END_REF].

Summary of the Results.

The first result we want to show is the ability of our method to track the behavior of link counts (total byte per unit time) over time. In Figure 2, we show the real and inferred link counts for our model. The evolution of the traffic and estimates are shown for a seven weeks duration for each observation vector. The calibration is applied only once, for the case where we analyze all the links at the same time (we have shown only the figures for this case). For the case of per-ink analysis, we have performed the calibration for each sequence separately. Now we are looking at the performance of our two methods, in the Abilene network. To validate our GMM model, we calibrate one GMM with a set of r components (r=2,3,4,5...) using the EM algorithm, and the decision to select the best model is done by analyzing the variance performed for each component in the mixture. The model with the lowest variance is chosen. We find in our computations, that the data residual can be organized into three (r = 3) distinct clusters. Thereafter A maximum a posteriori criteria is used to build, over the clusters, a finite alphabet of symbols, where we run the hidden markov model. To train the hidden markov model, one must ensure that that the different hidden states are well separated. This means that one should have a transition matrix with higher probabilities in its main diagonal. In our study we obtain two (2) states with the transition and observation matrices shown as follow (case where the The HMM results show clearly that, the state 1 is composed with almost entirely the symbol#3 and the state 2 is a mixture of two symbols (1 and 2) with 49% of probability of presence of symbol#1 and 51% for symbol#2. In the philosophy of anomaly detection theory, generally it is assumed, that anomalies might be rare; base on this assumption we can ask this question: if anomalous events occur, do they might come from state #1 or state #2 or both? To answer this question, we believe that all the changes in the mean of the residual (abrupt change, slow or high variations) can be tracked by a combination of symbols yielding in one state. In other words, one can think that one state must be classify as normal state and the other as abnormal state. Clearly, we believe that the state 1 containing the symbol #3 corresponding to the component with mean closely equal to zero might be labelled as the "normal state" and the state 2 as the "abnormal state". To confirm our intuition, we run the Viterbi algorithm for all the sequences of discrete alphabet and we discover one unique sequence composed only by the symbols in the state #2. In Figure 3, we show how the lower and higher variations in the residual are well tracked by the Viterbi path (composed only by symbol#1 and symbol#2). For the purpose of plotting, we multiply each symbol by a same constant). At this time, we can argue that, if anomalies exist they might be caught either by the two symbols simultaneously, either by one symbol only. In addition, in our study we discover that anomalous events never happen in the cluster (#3) with mean closely equal to zero. Combining this theory with a basic variance test performed for each data point, allows us to detect, identify and track anomalies in the different links. Among the 41 links, only 13 have drawn anomalies. In the Figure 5 we plot the tracking operation for some links. For all the figures, the red points corresponds to the symbol#1 and the blue ones to the symbol#2. The solid blue line is 1 In the case of per-link analysis we obtain an transition matrix and an observation matrix for each link the threshold for the part of residual corresponding to the symbol#1, and the dashed red line is the threshold for the part of residual corresponding to the symbol#2. We run the HMM more than one time to show that in some cases, the Viterbi path can be a sequence with only the presence of one symbol. And all anomalies yielding in the link in question are tracked inside the residual corresponding to that symbol. We can see some examples in Figure 4. In these cases, we can plot the results in continuous time, what we can not do in Figure 5. In our running scheme, we found that the coefficient that we multiply the standard deviation (to fix the threshold) of the residual variance lies in the interval [2 : 6]. Another result is about the performance obtained in analyzing the trade-off between false positive and false negative rates. We examine the entire traffic for each method at two levels: in level 1 we run the methods using one link at a time and in level 2 we consider all the links at the same time. We then can compute one false positive percentage and one false negative percentage for each threshold configuration scheme. The performance of the two methods on the Abilene data is depicted in the ROC (Receiver Operating Characteristic) curve of Figure 6. In Figure 6(a), we can see clearly that: i) when the analysis is done on a perlink basis, the hmm-variance-based anomaly detection performs better than a basic variance analysis. We can see that with a false positive rate (FPR) of 1.3% the HMM misses no anomalies (100% detection rate), while the variance test catches about 80% of the anomalies for the same FPR. ii) in the Figure 6(b), we show that analyzing all the links count simultaneously is more advantageous than taking account for one link at a time. One can observe that the HMM catches all the anomalies (0% of false negative) with a false positive rate of 0.006, while in the same time the variance test catches 85% of anomalies with the same FPR.

Conclusion

Our work has shown the robustness of combining mathematical tools for Anomaly detection issue. We begin by performing a Kalman filter on data selected in granularity of link counts. Then a GMM is used to organize the residual into different gaussian components, since we have as- sumed that the real distribution of the innovation process is a mixture of normal densities. Thereafter a MAP criteria is used to build an alphabet in which we perform a hidden markov model to find the temporal dependencies above the different gaussian components. Finally, we use the Viterbi algorithm to show that anomalous events can be tracked using a single sequence of our HMM states where we perform a variance test for each data point. In our study, we show that our anomaly detection methodology can provide appropriate management tools and can help network administrators to surveying their environment, since it can achieve higher detection rate with few false positive and the results are more convenient in the case of all-links analysis than in per-link based analysis. We have found that anomalies are capture by a part of the alphabet and one can ask this question: what about the remaining symbol (here symbol #3) ? Is it uncertain that anomalies may never happen in the cluster with mean closely equal to zero ? At the other hand, we must notice that the determination of the constant used to fix the variance threshold is heuristic and not straightforward. One would prefer the fixation of this value in an obvious and automatic manner. We will point out these ideas and questions in future works.

Figure 1 .

 1 Figure 1. Architecture of the monitoring system.

Figure 2 .

 2 Figure 2. Real(red) and estimated (blue) links obtained using Kalman filter.

Figure 3 .

 3 Figure 3. Track of the variations in the residual with the Viterbi path.

Figure 4 .Figure 5 .

 45 Figure 4. Anomalies detection in 2 links using only symbol 1.

Figure 6 .

 6 Figure 6. ROC curve using Abilene data.


Daniela Brauckhoff
email: brauckhoff@tik.ee.ethz.ch

Kave Salamatian

Martin May
email: may@tik.ee.ethz.ch

Applying PCA for Traffic Anomaly Detection: Problems and Solutions

Spatial Principal Component Analysis (PCA) has been proposed for network-wide anomaly detection. A recent work has shown that PCA is very sensitive to calibration settings, unfortunately, the authors did not provide further explanations for this observation. In this paper, we fill this gap and provide the reasoning behind the found discrepancies.

First, we revisit PCA for anomaly detection and evaluate its performance on our data. We develop a slightly modified version of PCA that uses only data from a single router. Instead of correlating data across different spatial measurement points, we correlate the data across different metrics. With the help of the analyzed data, we explain the pitfalls of PCA and underline our argumentation with measurement results. We show that the main problems that make PCA difficult to apply are (i) the temporal correlation in the data; (ii) the non-stationarity of the data; and (iii) the difficulty about choosing the right number of components. Moreover, we propose a solution to deal with the most dominant problem, the temporal correlation in data. We find that when we consider temporal correlation, PCA detection results are significantly improved.

I. INTRODUCTION

Principal Component Analysis (PCA) has been first proposed as a method for traffic anomaly detection in [?]. While being known in other domains before, Lakhina et al made its application very popular in the networking community. Subsequent publication confirmed the excellent performance of PCA and proposed extensions to it [?]. Only recently, it has been shown by Ringberg et al [?] that PCA is very sensitive to its parameter settings. The authors have reported about instability problems encountered when using PCA, however, they failed in providing precise reasons for their observation.

In this paper, we will provide the missing explanations for the encountered problems. In particular, we revisit PCA-based approaches for anomaly detection from a signal processing point of view. During the application of the PCA method on our dataset, we found similar inconsistencies as those reported in [?]. Further investigating the results, we found that the main problem of PCA, as used today, is that it does not consider the temporal correlation of the data. The fact that the data is temporally correlated indeed breaks the underlying hypothesis of the PCA-based anomaly detection method. The main contributions of this paper are that we (i) show what kind of problems arise when PCA is not carefully applied to anomaly detection; and (ii) provide a profound theoretical explanation for the encountered problem; and (iii) provide correction to the previously published method and alter it to a efficient anomaly detection mechanism. We validate our improved method by applying it to real network traffic with well known and identified anomalies.

The paper is structured as follows. In section II, we revisit the basics of anomaly detection methods and introduce the dataset used throughout this paper. Then, we report on the application of PCA to this data. The first results show very poor performance of the detector and by investigation the reason for it, we find the some basic features necessary for anomaly detection are not provided with PCA. For example, we find that the decision variable does not have a mean of zero -unfortunately PCA should be applied to zero mean random variables only. We propose two mechanisms that help to get the mean closer to zero and we'll show how these mechanisms are applied to the dataset.

In section III, we revisit the PCA theory and develop its basic properties. We thereafter extend the analysis to stochastic processes and explain why the simple PCA is not applicable and should be replaced by Karhunen-Loeve (KL) expansion. We describe this expansion an develop a Galerkin-based approach to calculate the KL expansion from a finite number of samples of the process. We then use the KL expansion to develop the predictive model that could be used for anomaly detection. We present the characteristics of the predictive models and propose a new methodology for choosing the optimal number of components in the KL expansion.

Section IV applies the KL extension method developed in section III to our data and examines therewith the reasons for the bad performance of classic PCA. We are able to validate that the source of the poor performance really is the temporal correlation. We further show that the non-stationarity is a critical issue and recalibration of the model is mandatory for good anomaly detection performance. Finally we describe the predictive models used for anomaly detection and open way for a signal processing approach to anomaly detection where the main challenge becomes to design filters that will be adapted to anomaly features.

II. PRACTICAL EXERCISE: APPLYING PCA FOR ANOMALY DETECTION ON OUR NETWORK

The main goal of our research is to develop efficient anomaly detection methods that are convenient for medium-sized ISPs (as for example, the SWITCH network, AS559). One of the most frequently used method for anomaly detection is the Principal Component Analysis based subspace method developed in [?], [?]. It was therefore, an obvious choice for a starting point to apply this method to our data. In this section, we give an introduction to anomaly detection in general, describe our dataset, and present the results of the PCA-based anomaly detection method on this data.

A. About anomaly detection

An anomaly detector consists essentially of two components: (i) an entropy reduction component and (ii) a decision component applying statistical tests to a decision variable issued from the first step. The entropy reduction step is here to simplify the second step.

A predictive model that forecasts the value of the parameter to monitor as X[k] is frequently used to build the entropy reduction. This results in an error signal that has a smaller entropy than the initial signal, but still retains the information required for anomaly detection. A decision variable D[k] used in the second step, is derived as a function of the prediction error e

[k] = X[k] -X[k].
In PCA based method, we build the predictive model assuming the projection in an orthogonal subspace obtained through application of PCA is a good predictor of the signal (see section III for a detailed description).

In the second step, we apply a statistical test to a decision variable that depends on the prediction error. Generally, the statistical test step checks the null hypothesis {H 0 : the decision variable is compatible with a prediction in conformity with observation}, i.e.,we reject the null hypothesis if the likelihood of the decision variable is below a threshold. The statistical test needs to have the distribution of the decision variable to obtain its likelihood. Two types of errors can occur during the statistical test: false negatives when one assumes that the null hypothesis is valid when there is an anomaly (the likelihood of observing the given variable being smaller than the threshold even if they have been an anomaly), and false alarm when deciding to refute the null hypothesis when there is no anomalies. The Neyman-Pearson theorem about statistical tests [?] defines a fundamental trade-off between false alarm probability and true negative probability; larger thresholds lead to lower false negative rates but larger false alarms rate and smaller thresholds result in higher false negative rates and smaller false alarms rate. The Receiver Operating Characteristics (ROC) Curve combining the two parameters in one value [?] captures this essential trade-off. The ROC curve is indeed a good performance metric for the two steps of an anomaly detector. A good entropy reduction technique is one that generates a decision variable with a good detection vs. false alarm rate tradeoff, i.e.,achieves a high detection rate with few false alarms.

For PCA based anomaly detectors, [?] proposed to use a non-linear function of the squared error

Q[k] = e[k] T e[k] as decision variable D[k] D(Q[k]) = Q[k] θ h (1)
[scale=0.5]plots/ROCcurve1.eps Fig. 1. ROC curve resulting from application of PCA method to three variation of traffic metrics. This non-linear function has been tailored to make its distribution converging as closely as possible to a gaussian distribution. In [?], the authors refer to the work of Jensen [?] showing that, under the hypothesis that elements of e[k] are statistically independent and follow a gaussian distribution, the distribution of D(Q[k]) converges to a normal distribution with known mean and variance. [?] gives formulas for deriving from PCA characteristics the mean and variance of D[k], as well as the parameters θ and h. It is therefore, possible to normalize D[k] to a gaussian random variable with zero mean and a variance of one. In the forthcoming, we will use the normalized D[k] as decision variable. Because of space restrictions, we are not giving all the details of this formula in this paper, but refer the reader to [?], [?] and [?].

B. Data Set

We use for our experiments 3 weeks of Netflow data coming from one of the peering links of a medium-sized ISP (SWITCH, AS559). These data were recorded in August 2007 and comprised a variety of traffic anomalies happening in the daily operation such as network scans, denial of service attacks, alpha flows, etc. In this dataset we distinguish between incoming and outgoing traffic, as well as UDP and TCP flows. For each of these four categories, we computed seven commonly used traffic features: byte, packet, and flow counts, source and destination IP address entropy, as well as unique source and destination IP address counts. All metrics were obtained by aggregating the traffic at 15 minute intervals resulting in a 28 × 192 data matrix per measurement day.

We applied the anomaly detection methods to a matrix x of size t × f with t observations from f features. We used the first two days (394 time samples) of this dataset for model calibration. Thereafter, we applied the anomaly detection method to all time samples. anomalies in the data were using available manual labeling methods: visual inspection of time series and top-n queries on the flow data. This resulted in 28 detected anomalies event in UDP and 73 detected in TCP traffics.

It is noteworthy that our dataset is different in nature from the one used in [?], [?]. There, the data were collected from different network border routers. Also, [?] only analyzes the traffic volume. In the follow-up publication [?], the authors extended their dataset with the 4 entropy values of the source and destination IP addresses as well as the source and destination port numbers. In this work, we collected the data at a single link. As in the work of Lakhina et al., the observed metrics are correlated in time and in space. Our observations and statements are valid whenever both spatial and temporal correlation exist; they can be extended to the case studied by [?], [?].

C. PCA application results

We then applied the methodology described in [?] to our data. With this method, we first apply the classical PCA as described in section III-A to the vector of metrics x[1 : 194] containing the first two days of metrics. For this purpose, we derive a spatial correlation matrix as Γ = 1 193 x[1 : 194]x[1 : 194] T . Then, we apply the SVD decomposition to the data, resulting in a basis change matrix. We construct a model using only the 8 top principal components (out of 28 possible). We choose component numbers following the methodology proposed in [?], [?] to englobe more than 95% of the variance in the initial metric (for detailed description see section III). We compute the term

Q[k] = e[k] T e[k]
from prediction errors, and we derived the normalized decision variable D(Q[k]) as described above.

Because of lack of space, we show only the results for the UDP traffic. Results for TCP traffic are analogous and are available in [?]. The ROC curve obtained by strict application of the PCA method as described in [?] leads to the figure labeled "Before mean removal" in Fig. 1.

Unexpectedly, the ROC curve showed very bad anomaly detection performance for the classic PCA method ("Before mean removal"). At best, we achieve a detection rate of 90% and this comes with the burden of 49% of false alarms. We detect only 15% of anomalies without any false alarm. The results obtained for TCP are even worse. Clearly, with such a performance this anomaly detection method is not useful in practice. After checking all the steps and ensuring that we closely followed the methodology proposed in [?], we had to back trace all causes that could have led to the bad anomaly detection performance.

On general terms, the results are in accordance with what was reported in [?]. It shows that tuning PCA to operate effectively in practice is difficult. Specifically, it shows high sensitivity of the false positive rate to the number of principal components chosen and relates this to the possibility that large anomalies pollute the normal subspace. Even [?] acknowledged this last point when it used another decision variable in place of the one defined in Eq. 1 to deal with the normal space pollution. However, [?], [?] failed to identify the causes of the problems they described. Particularly, a misleading term generated some confusion. In linear filter theory, any linear filter separates the space of signals into two orthogonal and dual subspaces: a subspace of the signals filtered by the filter and a subspace of the signal that pass through the filter. Calling the PCA method the subspace method is a kind of tautology. Any entropy reduction technique based on a predictive linear filter is a subspace method; PCA is just one subspace method with a very particular subspace structure. This means that in all methods, a component of the anomaly that maps to the normal space can pollute it. This is why the misdetection/false alarm trade)-off occurs! Stating that the PCA method has stability problems because of normal subspace pollution is a trivial statement as any failing prediction-based anomaly detector is failing because of normal subspace pollution! Analyzing if PCA is a suitable method or not requires less trivial reasoning. One need first to evaluate the assumption of the PCA method itself to ensure that the problem is not coming from bad application of the method. This is indeed the core contribution of our research.

D. Fault diagnostics of classic PCA

We began to back trace the inconsistencies by verifying if the normalized anomaly decision variable (Fig. 2) was following our theoretical expectations. We were expecting to have a decision variable with mean zero and variance 1. Actually, we observed that the mean of the decision variable over the 192 values used for calibration was 10.58 in place of zero! This bias could elucidate why the anomaly detection performance was bad as we apply the thresholds with an assumption of zero mean. Moreover, we found (see Fig. 2) that the convergence to a gaussian distribution is dubious, especially for the tails of the distribution where anomaly detection methods are looking at.

The first possible explanation was a well known argument, that was not stated in previous papers on anomaly detection: PCA should be applied to zero mean random variables. Meaning that one needs to first get the mean as close as possible to zero before applying PCA. We describe later in this paper how a mean close to zero can be obtained even for real-time operation. We verified if ignoring the non-zero-mean could have lead to our observations. We show in Fig. 1 the ROC curve obtained by applying PCA to the metrics after reducing the mean value. This second ROC curve shows a significant improvement of the false alarm rate. But still, the detection rate performance is below expectation. We checked the bias in the decision variable and found that it is reduced but is still important as the mean remains equal to 4.8 in place of zero.

A closer look at the proof of convergence in [?] shows that even if the convergence is robust toward non-gaussianity of the underlying variables, it is heavily dependent on the independence condition between the terms in e [k]. We verified if we can validate this condition. We estimated the correlation between the error terms and found that there is a high correlation between some of these terms. So the hypothesis of independence should be fully rejected! We conclude that the observed bias comes therefore from the lack of independence between error terms.

However, PCA theory predicted that these errors terms should be independent. By pursuing further the investigation, it appeared that the problem is coming from the method proposed in [?], [?]. They applied PCA to the spatial correlation alone (see above the formula for obtaining Γ) and did not remove the temporal correlation that is intrinsic in the used data. This temporal correlation leads to the prediction error terms and hence to the observed correlation.

E. Removing daily and weekly pattern

To check if we were going in the right direction, we tried to reduce the temporal correlation by removing the daily and weekly pattern. The observed metrics exhibited, as expected, a diurnal and a weekly pattern so that they were clearly nonstationary. Fortunately, this pattern is highly predictable. We used the first week of traffic as a reference to derive a weekly profile. We computed it by applying a smoothing filter that consisted of a centered moving average windows of size 20, to the observed metrics during a single week. This resulted in a smoothed version of the traffic metrics that we used as weekly profile. We subtracted this profile from the observed metrics to remove the daily and weekly pattern. However, the resulting signal has still a non-zero mean. We explained earlier that it is essential to apply PCA on zero mean signals. Whenever removing the mean for a given dataset is a trivial task, it becomes more difficult to remove for an online signal in real time. For this purpose, we used a zero-mean filter that estimates the mean from a window of 20 past samples and subtracts it from the signal. This yields signals that have a mean very close to zero.

We applied the same PCA method as above to the metrics after removing the weekly pattern and the reduction of the mean. This resulted in the third ROC curve also plotted in Fig. 1. This ROC curve shows a major improvement in anomaly detection quality. We achieve a detection rate of above 95% with 22% of false alarm rate. Moreover, the bias in the decision variable decreased from 10.58 to 2.2. This validated the hypothesis that the problems we were seeing are coming from the temporal correlation. A verification of the independence condition on error terms showed that the correlation decreases significanlt, but it remains too high to accept the independence assumption. Similar results for TCP traffic can be found in [?].

As explained before, the main goal of the entropy reduction step is to generate a decision variable that could thereafter lead to suitable ROC curves. It seems that applying PCA as developed in [?], [?] fails to achieve this goal. The main reason is the inability of PCA to generate a set of independent error terms. In the following, we revisit the theoretical foundation for PCA and examine PCA from a signal processing point of view.

III. WHAT IS PCA: A SIGNAL PROCESSING VIEW

A browse in the literature shows two closely related but different interpretations of PCA:

• As an efficient representation that transforms the data to a new coordinate system such that the projection on the first coordinate contains the greatest variance, the projection on second coordinate has the second greatest variance, and so on. • As a modeling technique using a finite number of terms of an orthogonal serie expansion of the signal with uncorrelated coefficients. Interestingly, the literature mainly motivates the application of PCA to networking anomaly detection by the first interpretation. However, this application is indeed following the two last interpretations. This has resulted in some erroneous interpretation and practices that have widely spread among the community. We are devoting this paper to describing these erroneous practices and to present a way of correcting them. Let's describe the above two interpretations with more details.

A. PCA : a suitable data representation

Let's suppose that we have a column vector of correlated random variables X = (X 1 , . . . , X K ) T ∈ R k . One observes these random variables through N independent realization vectors x i = (x i 1 , . . . , x i K ) T , i = 1, . . . , N and arranges them in a N × K observation matrix x with each row containing an observation vector x i . We are searching for the most "suitable" non-canonical basis (e 1 , . . . , e K ) for the vector space R K to represent the random variables X.

For the class of random variables that are linear (they can be decomposed to a linear combination of independent linear random variables) and have as sufficient statistics their means and covariances (i.e.,means and covariances entirely describe their joint probability distributions.), the most suitable basis is the one that maximizes the variance of each projected component. One very popular case where these two assumptions hold is when (X 1 , . . . , X K ) are jointly gaussian. Nonetheless, the literature is full of examples where using such an orthonormal basis results in erroneous interpretation because the linearity or the sufficiency of mean and covariance is not valid.

Under assumption of linearity and sufficiency of mean and variance, the most suitable basis is (φ 1 , . . . , φ K ), where φ i is an eigenvector of the covariance matrix of X defined as

Σ = ❊ (X -µ)(X -µ) T ) (µ is a column vector containing
the means of X i ). We derive these eigenvectors by solving the following linear equation:

Σφ i = λ i φ i (2)
where λ i are the eigenvalues of the covariance matrix. As the covariance matrix is positive definite, this equation has at most K positive eigenvalues and K different orthonormal eigenvectors. The basis change matrix U = [φ 1 , . . . , φ K ] contains in its columns, the eigenvectors φ i . Solving the above problem is called in matrix theory the Singular Value Decomposition (SVD) of the covariance matrix.

It is noteworthy that U is a basis change matrix only when X is zero mean, and in general one has to work with X = Xµ in place of X, i.e.,the coordinate change is ỹ = U x. This last point is frequently overlooked in the literature, and not taking care of it could lead to large errors when using PCA1 . In the forthcoming we will assume that we have taken care of this obvious precaution so we can drop the ∼. For real time operation, removing the mean can be approximated by using a zero-mean filter as described in section II-E.

After applying PCA one can rewrite the initial vector of random variables X in the new coordinate system as:

X = K i=1 Y i φ i (3)
where Y i are jointly independent random variables with ❊ Y i = 0 and ❱ar Y i = λ i . PCA replaces the correlated random variables X by a vector of independent random variables Y that are linearly equivalent. The independence of Y i is therefore an essential property as this is the main reason that PCA representation is "suitable". The above discussion remains theoretical, and in practice one has a set of observations and has to pick the suitable basis. Whenever the dataset under study is not flagrantly in contradiction with the conditions of mean and variance sufficiency, and linearity we can apply PCA and find a convenient representation of the data.

First one has to estimate the covariance matrix using the popular sum of product formula Σ = 1 N -1 x T x. Because of independence between observations this gives a reliable estimation. Thereafter applying the SVD factorization is just a straightforward and mechanical step that could provide the needed basis as well as the basis transform matrix.

B. Extension of PCA to a vector of stochastic processes

Before describing PCA as a modeling technique, let's extend it to stochastic processes. This extension is mandatory as the signals used for anomaly detection are samples of stochastic processes that have temporal as well as spatial correlations. Let's assume we have a K-vector of zero mean stationary stochastic processes X(t) = (X 1 (t), . . . , X K (t)) T with a covariance functions σ i,j (τ ) = ❊ X i (t)X j (tτ ) defined over an interval [a, b]. The extension 2 to multi-dimension of the Karhunen-Loeve (KLT) theorem [?] states that one can rewrite these processes as a serie expansion (named the KL expansion):

X l (t) = K i=1 ∞ j=1 Y l i,j Φ i,j (t) (4)
where Y l i,j are pairwise independent random variables and Φ i,j (t) are pairwise orthogonal deterministic (non-random) functions defined on [a, b] , i.e., b a Φ i,j (t)Φ * m,n (t)dt = 0 for i = m or j = n. Generally, the basis function Φ i,j (t) are rescaled such that b a |Φ i,j | 2 (s)ds = 1. This theorem extends the PCA to a vector of stochastic process. Eq. 4 is the equivalent of Eq. 3. The family of deterministic functions Φ i,j (t) is an orthonormal basis for the space of linear stochastic processes and the random variables Y l i,j are coordinates of the stochastic process X l (t) in this new space. We can formally derive the basis functions Φ i,j (t) by solving the following set of linear integral equations:

K i=1 b a σ i,l (s)Φ i,j (s -t)ds = λ l,j Φ l,j (t), j > 0.
These complex equations are simply the equivalent of equations 2. The random variables Y l i,j are obtained by projecting 2 The KL theorem was initially defined for one dimensional stochastic processes each stochastic process over an eigenfunction:

Y l i,j = b a X l (s)Φ i,j (s)ds
The KL expansion considers the temporal correlation (between time t and t + τ ) as well as the spatial correlation (between process X i (.) and X j (.)). This results in a more complex analysis than the simple PCA described earlier. However, this higher complexity is unavoidable because of the temporal correlation. Not taking it into account could lead to the errors described in section II-C.

In practice, we have only access to a finite set of samples observed each T time unit from the vector of stochastic processes and we have to implement the KL expansion only using them. In the forthcoming, we will use the notation [k] to represent the discrete version of a time continuous process sampled at time kT . Let's assume that we have n samples of the multidimensional stochastic process and the covariance values σ i,j (τ ) can be assumed as negligible for τ > N T . We can therefore truncate the KL expansion to N terms in place of the infinite number of terms needed normally. The Galerkin method [?] transforms the above integral equations to a matrix problem that could be solved by applying the SVD technique. This makes possible the derivation of KL expansion using only a finite number of samples. The Galerkin method generates a set of eigenvectors in a KN dimensional vector space, that are time-sampled version Φ i,j [k] = Φ i,j (kT ) of the originally continuous function Φ i,j (t). Finally, we obtain a discrete version of the KL expansion as :

X l [k] = K i=1 N j=1 Y l i,j Φ i,j j[k]. (5)
We first have to estimate the spatio-temporal correlation matrix. Let's construct a KN × (n -N ) observation matrix:

x =               x 1 (1) . . . x 1 (n -N ) x 1 (2) . . . x 1 (n -N + 1) . . . . . . . . . x 1 (N ) . . . x 1 (n) x 2 (1) . . . x 2 (n -N ) . . . . . . . . . x 2 (N ) . . . x 2 (n) . . . . . . . . . x K (1) . . . x K (n -N ) . . . . . . . . . x K (N ) . . . x k (n)               The matrix Σ = 1 n-N -1 x T
x is a KN × KN matrix that contains all the needed spatio-temporal covariance estimates. It is noteworthy that because of temporal correlation one needs more data to estimate correctly the covariance here than for the indepedent case we had in section III-A.

The Galerkin method consists of applying PCA as described in section III-A to this large matrix. This results in KN eigenvectors of length KN Φ i,j [.] that are used to construct a basis transform matrix U . The coefficients Y l i,j are obtained by applying the basis change transform y = U x. As can be seen applying KL expansion to K stochastic process entails diagonalizing a KN × KN matrix (in place of a K × K matrix in section III-A). However, this added complexity is unavoidable when one has to deal with correlated observations.

C. PCA as a modeling method

Up to now, we described the KL expansion as a tool for creating an equivalent (in probability) and suitable representation of vector of stochastic processes. If we neglect some of the smaller terms of the expansion (terms with small values of ❱ar Y l i,j ), we obtain a linear approximation of the initial process in a smaller dimension vector space. The discrete expansion in Eq. 5 is therefore approximated as:

Xl (kT ) = L i=1 M j=1 Y l i,j Φ k i,j . (6)
where M < N and L < K. This approximation has a noteworthy optimality property. Among all approximations defined over a linear space of dimension LM , this is the linear approximation with the smallest approximation error variance (❱ar X(t) -X(t) . The basis change transform becomes a KN × LM matrix U LM that contains the LM eigenfunctions Φ i,j [.] in its columns. This is the theoretical basis to use the KL expansion as a non-parametric and generic technique for modeling a large class of processes where we cannot reject the linearity and sufficiency of mean and variance (see section III-A).

The non-parametric nature of the above modeling technique is simultaneously its strength and Achilles heel; a nonparametric method is not based on any precise form of the distribution (out of the linearity and the mean and variance sufficiency) meaning it is more robust. At the same time being non-parametric means that no prior knowledge can be incorporated into the model.

Before going further, let's first give some details about the obtained model. The expansion in Eq. 6 provides a synthesis method for generating an approximated process Xl [k] by a bank of M L filters with Finite Impulse Response equal to Φ i,j [k], k = 0; . . . , KN ; each filter being excited by the random variable input Y l i,j . By predicting the values of the realization of the KN random variables Y l i,j by applying the basis change matrix to observation X[.], we can use this synthesis filter as a predictive filter. This is the approach followed in PCA and KL expansion based anomaly detectors. The anomaly detection literature has foreseen this approach.

Let's assume that the K linear stochastic processes in vector X[k] are linear processes, i.e.,one can represent them using a dynamic state space representation as

X[k +1] = AX[k]+ǫ[k],
where X[k] is a KN dimension vector constructed by concatenating N vectors (X[k], . . . , X[k -N ]) and ǫ[k] is a vector of KN independent and identically distributed (iid) random variables. Now let's assume that the process vector X[k] is approximated by a finite KL expansion with LM terms. There is therefore a U LM basis transform matrix that maps

X[k] into the new coordinate : Y R [k] = U LM X[k] (Y R [k] being the reduced coordinate vector of dimension LM ). The inverse projection can be found through X[k] = U T M L Y R [k]
(as U T is a Hermitian matrix). By replacing this term in the dynamic state space representation and by algebraic manipulation we obtain the following state space model for the approximated processes:

Y R [k + 1] = A M Y R [k] + B M ǫ[k] X[k] = U M L Y R [k] + D M ǫ[k] (7)
where

A M = U T LM AU LM ; B M = U T LM [(AS + I) | -S] with S = -A -1 + U LM A -1 M U T LM ; D M = S.
The resulting model has LM uncorrelated state variables (in Y) in place of KN state variables in the initial process (X).

The model in Eq. 7 shows precisely the effect of removing some terms from the KL expansion: this replaces the initial system dynamic by a new system with the above parameters. We can easily compare the initial system dynamic to the obtained model through three equivalent comparison approaches: impulse response, zeros and poles location, and frequency response. We can compute the poles and zeros locations by deriving the Z-transform of the transfer function of the processes. The initial process transfer function is -A) -1 ; the N K poles of the initial system are the roots of the equation det(zI -A) = 0 and all the zeros are at z = 0. The approximated system transfer function is

X (z) E(z) = (zI
X (z) E(z) = U M L (zI -A M ) -1 B M + D M ;the poles being the M L roots of det(zI -A M ) = 0 and the zeros the roots of the equation det(B M + (zI -A M )U T M L D M ) = 0.
Having the poles and zero locations we have a full characterization of the initial and modeled system.

In practice, we frequently do not have the state-space model of the observed process to be able to apply equation 7. One can infer the state space model following the Maximum Likelihood method used in [?] and thereafter use the above described approach. Another simpler approach consists of using spectral estimation techniques [?] to derive the spectra from the initial process. This results in an estimated spectra for the initial process. We derive the approximated process X

[k] easily as X[k] = U T M L U M L X[k]
. By applying the same spectral estimation to the approximated signal one can compute the spectra for the approximation and compare it with the spectra of the initial process. We will illustrate this approach later in this paper.

This results in a criterion to select the number of terms to keep in the KL expansion. [?], [?] it suggested to choose enough terms to capture 95% of the variance of the initial process. However, this approach does not take care of the features that will remain in the approximate model. A complete criterion retains a large enough number of terms so that the approximated process exhibits essential features that are of interest for anomaly detection. We will illustrate this approach later in section IV-D. Another outcome of the model description is that it enables the model designer to use its a priori knowledge to introduce or emphasize features that he believes would be useful for anomaly detection. This enables a refinement approach where in a first attempt we derive a

[width=]plots/spe u dp.eps model through the non-parametric KL expansion method and refine it through inspection of its features. We can refine the model could by classical control theory techniques by adding/removing/moving some of the poles or zeros such that the frequency response or the impulse response has suitable properties.

IV. VALIDATION

In this section, we will present the results of applying the KL expansion to the dataset described in section II-B. We first confirm that the temporal correlation is indeed the source of the bad performance of PCA based anomaly detection on our dataset. For this purpose, we have derived the KL expansion of our data set after preprocessing it only by removing the mean of the metrics.

A. Validation of KL expansion

Fig. 3 shows the ROC curves obtained for different values of the temporal correlation range N . The figure is plotted in semilog to present a better comparison between different values of N . All ROC curves are obtained with enough terms to capture 95% of the total variance in the model. In all cases, no more than two expansion terms were neede for this setting. The comparison of ROC curves (see Fig. 3) shows a considerable improvement in performance of the anomaly detection with use of KL expansion with N = 2, 3 and thereafter a decrease for N = 4. The same is observed for N > 4. Particularly, for N = 3 on can detect 97% of anomalies are with no more than 16% of false alarm rate.

To verify if the problems described in section II-C are really coming from the bias in the decision variable, we checked the value of this bias. We observed that augmenting N decreases the bias from 4.8 when N = 1 to 0.4 when N = 3. We plotted in Fig. IV-A, the distribution and a normal plot of the normalized decision variable for the model with N = 3. With increasing N the bias decreases (even if the marginal gain saturates). Moreover, the correlation between the terms in e[k] decreases significantly with augmenting N . That validates our hypothesis that relates the bad anomaly detection performance of PCA to temporal correlation.

B. Effect of non-stationarity

It remains another issue to solve: for N ≥ 4 the performance decreases even if the bias in the decision variable decreases. We investigated this observation and found that most of the mis-detections for N ≥ 4 happen at the end of the second week and during the third week of our dataset while they were spread in the three weeks for N < 3. A possible explanation is the stationarity issue: when N increases, the model contains more parameters and becomes more sensitive to the stationarity of the traffic metrics. This means that [scale=0.48]plots/ROCcurve2.eps Fig. 3. ROC curves resulting from KL expansion to the dataset with only mean removal for different value of the temporal correlation range.

[width=0.9]plots/spe u dp2.eps we can expect that anomaly detection performance decreases with time and this decrease is more pronounced for larger N . To further check this explanation, we plotted in Fig. 5 the ROC curves obtained only on the first week of anomaly data for increasing values of N . These ROC curves are now following the theoretical intuition: with augmenting N the anomaly detection performance ameliorates. But, it still begins to decrease for N = 10. This is compatible with our previous hypothesis about the effect of non-stationarity on models with increasing N .

Non-stationarity and model recalibration are indeed important issues that need careful evaluation and analysis. It is out of the scope of this paper and will be the subject of a forthcoming paper.

C. Interest of weekly trend removal

We introduced in section II-E a method to remove the daily and weekly trend. We showed also in Fig. 1 that this approach has promising performance. By further applying the KL expansion to the metrics pre-processed with mean and weekly trend removal we obtain the ROC curves plotted in Fig. 6. For the sake of comparison, we added in the plot the best performance ROC curve derived before trend removal with N = 3. The ROC curves in Fig. 6 show that weekly trend removal alone attains a performance close to the best model obtained without weekly trend removal. One can achieve better performance by adding a KL expansion: one can detect all anomalies with only 15% of false alarms. We can detect 93% of them with less than 7% of false alarms. Moreover, one can detect 55% of anomalies without any false alarms. Analogously to what was described in section IV-B, by augmenting N to 3 the performance degrades. We investigated this and found results analogous to what was described earlier.

D. Analysis of the predictive models

We developed in section III-C an approach for analyzing the model resulting from KL expansion. Particularly, we described an approach to evaluate the model resulting from PCA or KL expansion based on inferring the signal spectra at the output of the predictive model. We will illustrate the approach on our data in this section. We are showing in Fig. 7 the spectra of our metrics estimated using the Yule-Walker approach [?]. We also provide the spectra of the same metric predicted by the model resulting from the KL expansion with different numbers of kept components. The figure shows that applying the popular heuristic that consists in maintaining enough components to englobe 95% of variance is not satisfactory as it results in [scale=0.48]plots/ROCcurve4.eps [scale=0.48]plots/ROCcurve3.eps Fig. 6. ROC curve resulting from KL expansion on the metrics after weekly trend removal. a spectra that is far from the initial spectra. Specifically the predicted signal contains a higher level of high frequencies. Moreover, the low frequency part of the spectra is not well modeled. This means that using the 5 component model will result in overestimating high frequencies and missing the anomalies that are happening in lower frequencies. However, using a 12 component filter gives a better approximation of the initial signal. This leads to a finer methodology for choosing how many components of the KL expansion to use. We should choose enough components to have a good approximation of the initial signal spectra particularly in the region of the spectrum where expected anomalies will occurs.

The pole-zero diagram of the initial sytem along with the a predictive model with 8 components are plotted in Fig. IV-D. The plot shows how the predictive model approximates the initial system by positioning its 8 poles and zeros in the zplane. If an anomaly detector designer wishes to detect a particular type of anomaly signals he might use the above Pole-zeros diagram and move slightly some of the poles or zeros to emphasize specific parts of the transfer function spectra to achieve the needed goal.

The approach here becomes similar to the design of an equalizer in sound processing: the predictive filter acts as an equalizer that should filter out essential features of the anomaly signal such that these features appears in the error signal. One can then use the predictive model resulting from the KL expansion and ameliorate it to be more sensitive to specific anomalies that will happen in specific regions of the spectra. This novel anomaly detector design methodology opens new perspectives for introducing a priori knowledge in the design process. We will follow up this approach in our forthcoming papers.

V. CONCLUSION

This paper began with a very practical problem: how to apply the popular PCA method in real world anomaly detection. We found that direct application of the PCA method results in poor performance in terms of ROC curves; we investigated the problem and found that the main source of the problem is the bias coming from correlation in prediction error terms. After a detailed theoretical analysis, it appears that the correct framework is not the classical PCA but rather the Karhunen-Loeve expansion. We have presented the KL expansion and have provided a Galerkin method for developing a predictive model. This method has thereafter been applied to data traces from the Switch network and we have shown that an important improvement is attained when temporal correlation is considered. We also have developed a methodology for designing anomaly detection predictive filters and have shown its application to real data.

It is noteworthy that this paper is not claiming that PCA is the only way to do anomaly detection. By showing that PCA is based on a predictive filter -just as the Kalman filter is -we open way for a signal processing approach to anomaly detection where the main challenge becomes to design filters that will be adapted to anomaly features. To the best of our knowledge, this is the first paper to introduce such an approach and to remove completely the curtain of black magic that draped the application of PCA to anomaly detection.

Fig. 2 .

 2 Fig. 2. Normal plot of normalized decision variable D[k] for UDP traffic obtaining over the data used for calibrating the PCA model

Fig. 4 .

 4 Fig. 4. Normal plot of normalized decision variable D[k] for N = 3

Fig. 5 .

 5 Fig. 5. ROC curves obtained over the first week of data for different values of N .

It is noteworthy that even if [?], [?] did not state clearly the necessity of removing the mean, they have used zero mean signals in their implementation code

ACKNOWLEDGMENTS

The authors would like to thank Mark Crovella for his helpful discussions that have led to some aspects of this paper. We also thank Anukool Lakhina for giving access to his code and the Abilene data he used in his papers.