Christian Callegari
email: christian.callegari@iet.unipi.it

Sandrine Vaton
email: sandrine.vaton@enst-bretagne.fr

Michele Pagano
email: m.pagano@iet.unipi.it

A New Statistical Approach to Network Anomaly Detection

Keywords: Intrusion Detection System, High Order Markov Chain, Mixture Transition Model, Non-Homogeneous Markov Chain

In the last few years, the number and impact of security attacks over the Internet have been continuously increasing. To face this issue, the use of Intrusion Detection Systems (IDSs) has emerged as a key element in network security. In this paper we address the problem considering a novel statistical technique for detecting network anomalies. Our approach is based on the use of different families of Markovian models (namely high order and non homogeneous Markov chains) for modeling network traffic running over TCP. The performance results shown in the paper, justify the proposed method and highlight the improvements over commonly used statistical techniques.

I. INTRODUCTION

In the last few years Internet has experienced an explosive growth. Along with the wide proliferation of new services, the quantity and impact of attacks have been continuously increasing. The number of computer systems and their vulnerabilities have been rising, while the level of sophistication and knowledge required to carry out an attack have been decreasing, as much technical attack know-how is readily available on Web sites all over the world.

Recent advances in encryption, public key exchange, digital signature, and the development of related standards have set a foundation for network security. However, security on a network goes beyond these issues. Indeed it must include security of computer systems and networks, at all levels, top to bottom.

Since it seems impossible to guarantee complete protection to a system by means of prevention mechanisms (e.g. authentication techniques), the use of an Intrusion Detection System (IDS) is of primary importance to reveal intrusions in a network or in a system. IDSs are usually classified on the basis of several criteria [START_REF] Kemmerer | Intrusion Detection: a Brief History and Overview[END_REF].

State of the art in the field of intrusion detection is mostly represented by misuse based IDSs. Considering that most attacks are realized with known tools, available on the Internet, a signature based IDS could seem a good solution. Nevertheless hackers continuously come up with new ideas for the attacks, that a misuse based IDS is not able to block. This is the main reason why our work has focused on the development of an anomaly based IDS. In particular our goal is to reveal intrusions carried out exploiting TCP bugs, by using Markovian models (high order and non homogeneous Markov chains) to describe the behavior of network traffic.

The use of first order homogeneous Markov chain is a wellknown approach to detect two distinct kinds of "anomalies": masqueraders (analyzing the command stream of a host) and intruders (analyzing the evolution of TCP flows in the network traffic) [START_REF] Ye | Robustness of the Markov-Chain for Cyber-Attack Detection[END_REF].

Vardi and Ju in [START_REF] Ju | A Hybrid High-order Markov Chain Model for Computer Intrusion Detection[END_REF] describe the use of high order Markov chains to detect masqueraders at the host level, and in [START_REF] Schonlau | Computer Intrusion: Detecting Masquerades, NISS[END_REF] [START_REF] Ye | First-order Versus High-order Stochastic Models for Computer Intrusion Detection[END_REF] the authors compare performance of first order models and "generic" high order models.

After an extensive survey, to the best of our knowledge, there is no work directly related neither to the use of high order Markov chains to detect anomalies in the TCP traffic nor to the application of non-homogeneous Markov chains to anomaly detection in general.

Moreover no study at all compares the performance achievable with Markov chains of different orders and with a simple "independent" model..

The paper is structured as follows: next section provides a detailed description of the implemented system, while the subsequent section presents the experimental results. Finally section 4 concludes the paper with some final remarks..

II. SYSTEM DESIGN

In this section we provide a detailed description of the proposed anomaly based NIDS.

The aim of our work is to perform a comparison between several statistical models, which can be used to describe the behavior of TCP connections. More in detail we take into account the use of:

• first order homogeneous Markov chains

A. Training Phase

To build the model which represents the "normal" behavior of the network, the system needs a training phase during which it analyzes some network traffic, supposed to be attack free. The system analyzes raw traffic traces in libpcap format, the standard used by publicly available packet sniffer software, as Tcpdump or Ethereal. First of all the IDS performs a filtering phase so that only TCP packets are passed as input to the detection blocks.

The IDS only considers some fields of the packet headers, more precisely the IP source address, the IP destination address, the source port number, the destination port number, and the TCP flags. The IP addresses and the port numbers are used to identify a connection, while the value of the flags is used to build the profile. Experimental results have shown that the stochastic models associated to the different applications strongly differ one from the other. Thus, before constructing the model, the system isolates the different services, on the basis of the server port number, and the following procedure is realized once per each application. After that the IDS reconstructs the single connections on the basis of the 5-tuple (source and destination addresses, source and destination ports, and protocol).

A value s i is associated to each packet, according to the configuration of the TCP flags:

s i = syn + 2 • ack + 4 • psh + 8 • rst + 16 • urg + 32 • f in (1)
Thus each "mono-directional" connection is represented by a sequence of symbols s i , which are integers in {0, 1, • • • , 63}.

The training phase, as well as the detection phase, varies according to the stochastic model we are taking into account.

1) ECDF: In the case of the stationary ECDF the training phase simply consists of evaluating the probabilities P(s i ) that the TCP flags assume the value s i , independently of the position of the packet in the TCP connection.

For the non-stationary ECDF the system has to compute the probabilities P j (s i ) that the TCP flags of the j th packet of the connection assume the value s i . Taking into account the nature of the security attacks, for reducing the complexity of the system, we have decided to evaluate such probabilities only for the first 10 packets of a connection, i.e. j = 1, 2, • • • , 10.

2) Markov Chains: In the case of Markovian models the symbols s i are considered as the states of a hidden discrete time finite state Markov chain.

Since not all the TCP flags configurations are observable in real traffic, the system only considers the states observed in the training phase. Moreover, to take into account the possibility that some new flags configurations could be observed during the detection phase, a rare state is added. This procedure allows us to reduce the cardinality of the state space from 64 (all the possible configurations of the six TCP flags bits) to a number K, usually smaller than ten.

Then the system estimates the transition probabilities of the Markov chain.

Since the computation of such probabilities is quite straightforward in the case of first order Markov chains (homogeneous and non-homogeneous), in the following we consider a Markov chain of order l. The main problem related to this kind of models is the "explosion" of the number of parameters, which grows exponentially with the order, according to the rule K l (K -1). This entails the need of a parsimonious representation of the transition probabilities. The approach used in this paper is the Mixture Transition Distribution (MTD) model, first proposed in [START_REF] Raftery | A model for high-order Markov chains[END_REF]. Under the MTD model, the transition probabilities of an l th order Markov chain can be expressed as follows:

P(C t = s i 0 |C t-1 = s i 1 ,C t-2 = s i 2 , • • • ,C t-l = s i l ) = ∑ l j=1 λ j r(s i 0 |s i j ) (2)
where C t represents the state of the chain at step t and the quantities

R = {r(s i |s j ); i, j = 1, 2, • • • , K} Λ = {λ j ; j = 1, 2, • • • , l} (3)
satisfy to the following constraints:

r(s i |s j ) ≥ 0; i, j = 1, 2, • • • , K ∑ K s i =1 r(s i |s j ) = 1 ∀ j = 1, 2, • • • , K (4)
λ j ≥ 0; j = 1, 2, • • • , l and ∑ l j=1 λ j = 1 (5)
A consequence of the use of the MTD model is the reduction of the number of parameters from K l (K -1) to K(K -1) + l -1. To take into account the presence of the "rare" state (labelled K), we have to fix the following quantities:

r(rare|s i ) = ε, ∀i = 1, 2, • • • , K and ε small (ε = 10 -6 ) r(s i |rare) = (1 -ε)/(K -1), ∀i = 1, 2, • • • , K -1 (6)
According to the MTD model the log-likelihood of a sequence

(c 1 , c 2 , • • • , c T ) of length T is LL(c 1 , c 2 , • • • , c T ) = ∑ K i 0 =1 • • • ∑ K i l =1 N(s i 0 , s i 1 , • • • , s i l )log ∑ l j=1 λ j r(s i 0 |s i j ) (7)
where

N(s i 0 , s i 1 , • • • , s i l ) represents the number of times the transition s i l → s i l-1 → • • • → s i 0 is observed. Maximum like- lihood estimation (MLE)
of the chain parameters requires to maximize the right hand side of eq. ( 7), with respect to R and Λ, taking into account the constraints (4) and [START_REF] Ye | First-order Versus High-order Stochastic Models for Computer Intrusion Detection[END_REF].

Since the original solution [START_REF] Raftery | Estimation and modelling repeated patterns in high-order Markov chains with the mixture transition distribution (MTD) model[END_REF] seems to be too much computationally demanding, we have applied the procedure proposed in [START_REF] Ju | A Hybrid High-order Markov Chain Model for Computer Intrusion Detection[END_REF], which consists in an alternate maximization with respect to R and to Λ. This process leads to a global maximum, since LL is concave in R and Λ. For the part when R is fixed, we maximize LL with respect to Λ, and vice-versa. In the first step (estimation of Λ) we have used the sequential quadratic programming, while the second maximization step (estimation of R) is a linear inverse problem with positivity constraints (LININPOS) that we have solved applying the expectation maximization (EM) algorithm [START_REF] Vardi | From Image deblurring to Optimal investments: Maximum Likelihood Solutions for Positive Linear Inverse Problem[END_REF]. Since the first maximization step is quite trivial, in what follows we discuss the second step, i.e. the estimation of the matrix R, with the vector Λ fixed. First of all we have re-indexed the loglikelihood in the following way:

φ(s i 0 , s i 1 , • • • , s i l ) = 1 + l ∑ j=0 (s i j -1)K l-j → k (8)
which takes to

N(s i 0 , s i 1 , • • • , s i l ) → a k and l ∑ j=1 λ j r(s i 0 |s i j ) → b k (9)
Thus, at first we estimate the quantities b k (MLE) and then we solve the linear system

b k = l ∑ j=1 λ j r(s i 0 |s i j ) (10)
which is a LININPOS problem. At this point, the log-likelihood can be expressed as:

k l+1 ∑ k=1 a k logb k ( 11
)
where [START_REF] Ju | A Hybrid High-order Markov Chain Model for Computer Intrusion Detection[END_REF]

k l+1 ∑ k=1 a k = T -l and k l+1 ∑ k=1 b k = K l (12)
Thus a simple Lagrange method argument shows that the loglikelihood is maximized when bk

= a k ∑ k a k ∑ k b k = a k T -l K l , ∀k (13)
or, equivalently, when

l ∑ j=1 λ j r(s i 0 |s i j ) = K l T -l N(s i 0 , s i 1 , • • • , s i l ), ∀(i 0 , • • • , i l ) (14)
Thus, if we consider these equations as a linear system subject to the constraints (4), we obtain a LININPOS problem, which can be solved, in the sense of the minimum Kullback-Leibler distance, using the EM algorithm. More in detail we have [START_REF] Iusem | A New Smoothing-Regularization Approach for a Maximum-Likelihood Estimation Problem[END_REF]

A B R = K l T -l • N 1 ( 15
)
where

R T = (r(s 1 |s 1 ), r(s 2 |s 1 ), • • • , r(s K |s 1 ), • • • , r(s K |s K )) = (r 1 , r 2 , • • • , r K 2 ) (16)
are the unknowns, r(s i |s j ) = r i+K( j-1) , and

N =                N(s 1 , s 1 , • • • , s 1 ) N(s 1 , s 1 , • • • , s 2 ) . . . N(s 1 , s 1 , • • • , s K ) . . . N(s K , s K , • • • , s 1 ) . . . N(s K , s K , • • • , s K )                =                N 1 N 2 . . . N K . . . N 1-K+K l+1 . . . N K l+1                (17)
where

N(s 0 , s 1 , • • • , s l ) = N i , i = φ(s i 0 , s i 2 , • • • , s i l ) A = {a i j } K l+1 ×K 2 where a i j = ∑ l k=0 λ k I[ j = i 0 + K(i k -1)], (i 0 , • • • , i l ) = φ -1 (i) (18)
The matrix B looks like:

B =      1, • • • , 1 0, • • • 0 0, • • • 0 0, • • • , 0 0, • • • 0 1, • • • , 1 0, • • • 0 0, • • • , 0 0, • • • 0 0, • • • 0 . . . 0, • • • , 0 0, • • • 0 0, • • • 0 0, • • • 0 1, • • • , 1      K×K 2 (19)
At this point the EM iteration step is the following:

r j ← a • j a • j + b • j r j (A, K l T -l • N, R) + b • j a • j + b • j r j (B, 1, R) (20)
where

r j (W , u, v) ≡ v j w • j ∑ i w i j u i ∑ k w ik v k , j = 1, 2, • • • , K 2 (21)
for matrix W = {w i j } and vectors u = {u i }, v = {v i }, and

a • j = ∑ i a i j = ∑ (i 0 ,••• ,i l )=φ -1 (i) l ∑ k=0 λ k I[ j = i 0 + K(i k -1)] = l ∑ k=0 λ k ∑ i 0 • • • ∑ i l I[ j = i 0 + K(i k -1)] = l ∑ k=0 λ k K l-1 = K l-1 (22) and b • j = ∑ i b i j = 1.
The choice of the initial values for R and Λ is a key point. Experimental tests have shown that good results are obtained choosing λ i = 1/l, i = 1, 2, • • • , l and setting R to the first order transition probabilities.

B. Detection Phase

Once the training phase has been performed, the IDS has a model of the "normal" behavior of the network, represented by the computed profile.

As for the training phase, the input is given by raw traffic traces in libpcap format, which are processed so as to extract sequences of TCP flags configurations.

Thus, given an observed sequence (c 1 , c 2 , • • • , c T ), the system has to decide between the two hypotheses:

H 0 : {(c 1 , c 2 , • • • , c T ) ∼ computed model} H 1 : {anomaly} (23)
The problem is to choose between a single hypothesis H 0 , which is associated to the estimated stochastic model, and the composite hypothesis H 1 , which represents all the other possibilities. No optimal result is presented in the literature about this decision theory problem, thus the best solution is represented by the use of the Generalized Likelihood Ratio (GLR) test [START_REF] Mood | Introduction to the Theory of Statistics 3rd ed[END_REF].

Since the problem is quite straightforward for ECDF, in the following we only consider the case of Markovian models, for which the GLR test is defined as follows:

H(X) = H 0 i f X < ξ H 1 i f X > ξ (24)
where the threshold ξ is set by means of MonteCarlo simulations and the quantity X is given by:

X = Max v =u L(c 1 , c 2 , • • • , c T |Λ v , R v ) L(c 1 , c 2 , • • • , c T |Λ u , R u ) 1 T (25)
where the vector (Λ u , R u ) represents the parameters corresponding to the model computed during the training phase (hypothesis H 0 ) and the component 1 T is introduced to take into account that each observed sequence is characterized by a different length T . It is worth noticing that this test is equivalent to decide on the basis of the Kullback-Leibler divergence between the model associated to H 0 and the one computed for the observed sequence.

III. EXPERIMENTAL RESULTS

In this section we compare the performance of the different statistical models over the 1999 DARPA evaluation project [START_REF] Lippmann | The 1999 DARPA Off-Line Intrusion Detection Evaluation[END_REF] data set.

For sake of brevity, in the following we only present the results related to the Telnet traffic, since they appear to be representative of the overall performance.

To test the correctness of the computed models we have calculated the log-likelihood function of some sequences. Figure 1 corresponds to a "normal" connection. As expected from the theory, the function decreases almost linearly with the number of packets; its slope is equal to the entropy of the model, which, for first order Markov chain, is defined as:

H (MC) = ∑ i ∑ j π(i)P(s j |s i ) log P(s j |s i )
where π(i) is the stationary distribution of the Markov chain. The given definition can be easily extended for higher order Markov chains.

On the other hand the effect of an anomaly is an abrupt jump in the log-likelihood function, as highlighted by figure 2. Both these figures refer to a first order model, but the behavior of the log-likelihood function does not significantly vary with the order of the Markov chain.

To evaluate the performance we have used a Receiver Operating Characteristic (ROC) curve, which plots detection rate vs. false positive rate, obtained varying the value of the threshold ξ.

Figure 3 shows the ROC curves for Markov chains of different orders. We have considered Markov chains of order up to 4, since higher orders imply a heavy processing time, not suitable for on-line detection. Since the results obtained using a model based on a Markov chain of order 1 are already very good for these traffic traces, it is not easy to realize that we achieve some improvements with high order models. To be noted that the ROC curves are almost ideal, since we have a detection rate close to 100% with a negligeable false alarm rate. Nevertheless the zoomed area inside the figure shows that with the model of order 4 we are able to achieve the best results, obtaining a detection rate of 53% with a false alarm rate which is about one half of that related to the Markov chain of order 1.

The following figure shows the performance of the ECDF, while figure 5 presents a comparison between the first order Markov chains and the time dependent models described in the paper.

Since a detection performed analyzing only the first 10 packets of each connection is obviously worse than the one based on the entire connections, also the time independent model has been computed only considering the first ten packets of each connection. It is easy to conclude that the homogeneous Markov chain achieves a detection rate almost 10% bigger than the other two models.

This apparent paradox can be justified by the fact that the non-homogeneous models have been computed with a relatively short, and so "incomplete", training phase. Indeed, on one side the whole training data set has been used to compute only one homogeneous model, while on the other side, the same quantity of data is partitioned into ten subsets corresponding to the first ten steps in the time evolution of each connection. In particular this can lead to almost deterministic probabilities for the first steps of the non homogeneous models, thus a single flag configuration at step i, present in the training data set only at steps j = i (and hence captured by the time independent model), may generate a false alarm.

Finally, we have taken into account that an intrusion should be detected as soon as the anomaly appears. Thus, in figure 6, we show the performance of the homogeneous Markov chain model as a function of the number of analyzed packets for each connection (both for building the model and for the detection phase). The results highlight that good performance are achieved with a small number of packets, demonstrating that such statistical models are suitable for on line anomaly detection.

IV. CONCLUSIONS

In this paper we have presented an anomaly based network intrusion detection system, which detects anomalies using statistical characterizations of the TCP traffic. We have compared several stochastic models, such as first order homogeneous and non-homogeneous Markov chains, high order homogeneous Markov chains, and stationary and non-stationary ECDF. We have detailed the estimation of the parameters of the models and we have shown the results obtained with the DARPA 1999 data set. The performance analysis has highlighted that the best results are obtained with the use of homogeneous Markov chains and that some improvements can be achieved using high order Markovian models: for instance, 4th order Markov chains lead to the same detection rate of first order models, with almost one half of false alarms.

Moreover, we have shown that, since only a small quantity of packets is sufficient to reveal intrusions in the TCP traffic, this kind of approach is suitable for on line detection.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Log-likelihood function of a "normal" connection

Fig. 3 .Fig. 5 .Fig. 6 .

 356 Fig. 3. Performance of Markovian models of different orders

V. ACKNOWLEDGMENTS

This work was partially supported by the Euro-NGI Network of Excellence funded by the European Commission and partly by the RECIPE project funded by MIUR.


Lionel Fillatre

Igor Nikiforov

Sandrine Vaton

Pedro Casas

Sequential Non-Bayesian Network Traffic Flows Anomaly Detection and Isolation

Keywords: Statistical change detection/isolation, nuisance parameters, network traffic flows, parsimonious model, invariant decision algorithm

Sequential detection (and isolation) of unusual and significant changes in network Origin-Destination (OD) traffic volumes from simple link load measurements is considered in the paper. The ambient traffic, i.e. the OD traffic matrix corresponding to the non-anomalous network state, is unknown and it is considered here as a nuisance parameter because it can mask the anomalies. Since the OD traffic matrix is not recoverable from the simple link load measurements, the anomaly detection is an ill-posed decision-making problem. The method discussed in this paper consists of finding a linear parsimonious model of ambient traffic (nuisance parameter) and detecting/isolating anomalies by using an invariant decision algorithm. An optimal sequential algorithm has been discussed in our previous publication, the main goal of the present paper is to discuss a simple "snapshot" algorithm based on the last vector of observations.

Introduction and motivation

The traffic demand over a network is typically described by a traffic matrix that captures the amount of traffic transmitted between every pair of ingress and egress nodes in a network, also called the Origin-Destination (OD) flows. A volume anomaly is a sudden change in an OD flow's traffic (for example, due to denial-of-service, viruses/worms, external routing reconfigurations, etc.) that spans multiple physical links of the network. The reliable detection/isolation of these unusual and significant changes in the OD traffic matrix is an important issue for network operation. The direct measurements (flow-level data) need high hardware requirements to be collected and processed network-wide [START_REF] Coates | Internet tomography[END_REF]. Consequently, the Simple Network Management Protocol (SNMP), which is a widely deployed standardized protocol, is preferred in practice to measure link loads and obtain some information on the traffic matrix. The challenge lies in the ill-posed nature of the problem: the number of unknown OD flows is much larger than the number of SNMP measurements [START_REF] Ahmed | Multivariate online anomaly detection using kernel recursive least squares[END_REF].

Several approaches are proposed to remedy this problem. The first group of methods consists in detecting anomalies in SNMP measurements without taking in consideration the traffic matrix. Such methods typically use times series, ARIMA models among others, [START_REF] Krishnamurthy | Sketch-based change detection: methods, evaluation, and applications[END_REF]; [START_REF] Thottan | Anomaly detection in IP networks[END_REF]; Tartakovsky et al. (2006aTartakovsky et al. ( , 2006b)); [START_REF] Ahmed | Multivariate online anomaly detection using kernel recursive least squares[END_REF]) to model the SNMP measurement evolution in time and to detect deviations. In [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF], the authors propose to decompose the SNMP measurements on a Principal Component Analysis (PCA) basis. These methods can detect anomalies by monitoring each link but they do not exploit the routing matrix, i.e. the linear mathematical relationship between the OD traffic matrix and the SNMP measurements.

The second group of methods [START_REF] Cao | Time-varying network tomography: router link data[END_REF]; [START_REF] Coates | Internet tomography[END_REF]; [START_REF] Tebaldi | Bayesian inference on network traffic using link count data[END_REF]) exploits this linear mathematical relationship. These approaches typically assume that the traffic matrix is well approximated by a known statistical model. Such a method requires a well known prior to be efficient, which is not always feasible in practice. In [START_REF] Zhang | Network anomography[END_REF], the authors study a large number of methods based on several kinds of model for SNMP measurements (wavelets, PCA among others) and OD flows (ARIMA time series) to detect anomalies. A major drawback of these methods is the absence of theoretical results on the optimality of the studied methods with respect to the statistical hypotheses testing theory: no results are proposed about the design of an optimal test which maximises the probability to detect anomalies.

Finally the last group of methods consists in using the Kalman filtering technique [START_REF] Soule | Traffic matrices: balancing measurements, inference and modeling[END_REF]) to track the evolution of the traffic matrix in time and to detect changes in the OD flows. Strictly speaking, the ill-posed nature of the measurement model makes the Kalman filter not observable and the Kalman filtering efficiency strongly depends on the initialization, which is a serious limitation in practice.

The contributions of this study are the following. First, a parsimonious linear model of non-anomalous OD flow volumes ("ambient" traffic) is proposed. This model can be used in two ways, either to estimate the OD flow volumes or to eliminate the non-anomalous "ambient" traffic from the SNMP measurements in order to provide residuals sensitive to anomalies. Secondly, since a few anomaly-free SNMP measurements (at most one hour of measurements) is sufficient to obtain a reliable model of the OD flows, the proposed method is well adapted to the case of non stationary measurements and to dynamic routing which is a typical situation in practice. Finally, an optimal invariant detection algorithm is proposed to detect anomalies directly from SNMP measurements. x(1, 2)

SNMP Measurements

x(1, 2)

x(1, 3)

x(1, 3) Let us consider a network composed of r nodes and n monodirectional links [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF]; [START_REF] Zhang | Network anomography[END_REF]) where y(ℓ) denotes the volume of traffic (typically in bytes) on the link ℓ at time t. To simplify the notations, the subscript t is omitted. The link volumes are provided by SNMP measurements. Let x(i, j) be the OD traffic demand from node i to node j at time t. This situation is shown in Fig. 1. The traffic matrix X = {x(i, j)} is reordered in the lexicographical order as a column vector X = (x(1), . . . , x(m))

T where X T denotes the transpose of the matrix X and m = r 2 is the number of OD flows. Let us define an n × m routing matrix A = (a(ℓ, k)) where 0 ≤ a(ℓ, k) ≤ 1 represents the fraction of OD flow k volume that is routed through link ℓ. This leads to the linear model Y = A X where Y = (y(1), . . . , y(n))

T . Without loss of generality, the known matrix A is assumed to be full row rank, rank (A) = n.

Problem statement : volume anomaly detection

The problem consists in detecting (/isolating) a significant volume anomaly in an OD flow x(i, j) by using only SNMP measurements y(1), . . . , y(n). For example, in Fig. 1, it is necessary to detect a sudden increase of the traffic volume x(1, 3) by using y(1), y(2), y(3) (typical figures of anomaly to be detected in SNMP measurements are presented in Lakhina et al. ( 2004)).

As it has been mentioned above, the main problem with the SNMP measurements is that n ≪ m. To overcome this difficulty a parsimonious linear model of non-anomalous traffic has been used. The derivation of this model includes two steps: i) description of the ambient traffic by using a spatial stationary model and ii) linear approximation of the model by using piecewise polynomial splines. The idea of the spline model is that the non-anomalous (ambient) traffic at each time t can be represented by using a known family of basis functions superimposed with unknown coefficients : i.e. X t ≈ Bµ t , where the m × q matrix B is assumed to be known and µ t ∈ R q is a vector of unknown coefficients such that q < n. Finally, it is assumed that the model residuals together with the natural variability of the OD flows follow a spatial Gaussian distribution, which leads to the following equation:

X t = Bµ t + ξ t (1)
where ξ t ∼ N (0, Σ) is a Gaussian noise, with the m × m spatial diagonal covariance matrix Σ = diag(σ 2 1 , . . . , σ 2 m ). The advantages of a parametric model are the following ones. First, a non-parametric basis, typically the PCA basis, can be used to generate the matrix B but this solution needs direct OD flow measurements (infeasible in practice) and the PCA basis depends on the period when the measurements are made (see [START_REF] Ringberg | Sensitivy of PCA for traffic anomaly detection[END_REF]). Secondly, the parametric detection method performances are better than the non-parametric ones provided that the adopted model is accurate enough. Hence, the link load measurement model is given by the following linear equation :

Y t = A Bµ t + Aξ t = Gµ t + ζ t + [η = possible anomaly], (2)
where Y t = (y(1), . . . , y(n)) T and ζ t ∼ N (0, AΣA T ). Without any loss of generality, the resulting matrix G = A B is assumed to be full column rank. Typically, when an anomaly occurs on OD flow j at time t 0 (change-point), the vector η has the form η = ǫ a(j) where a(j) is the j-th normalized column of A and ǫ is the intensity of the anomaly. The goal is to detect the presence of an anomalous vector η not explicable by the ambient traffic model X t ≈ Bµ t and to isolate the contaminated OD flow j.

The problem of the anomaly detection/isolation based on a sequential decision rule [START_REF] Nikiforov | A generalized change detection problem[END_REF][START_REF] Nikiforov | A simple recursive algorithm for diagnosis of abrupt changes in random signals[END_REF][START_REF] Nikiforov | A lower bound for the detection/isolation delay in a class of sequential tests[END_REF]) has been discussed in our previous publication (see Fillatre et al. (2007b)) for a time correlated noise sequence (ζ t ) t≥1 . In case of a strong autocorrelation in the residual process (ζ t ) t≥1 and a large signal-to-noise ratio, the statistical characteristic of the snapshot algorithm is very close to the characteristic of the optimal sequential test because the mean detection delay is almost equal to one observation. For this reason let us consider now a simple and efficient detection scheme based on the step by step hypotheses testing sometimes called "snapshot".

Snapshot detection : hypotheses testing

Since the matrix Φ = AΣA T is known, the testing problem consists of choosing between the two alternatives at time t:

H 0 = {Z t ∼ N (θ + Hµ t , I n ); θ = 0, µ t ∈ R q }
(3)

H 1 = {Z t ∼ N (θ + Hµ t , I n ); θ = 0, µ t ∈ R q }, (4)
with

Z t = Φ -1 2 Y t , H = Φ -1 2 G, θ = Φ -1 2 η, Φ -1 2
is the square-root matrix of Φ -1 , Φ -1 is the inverse of Φ and I n is the identity matrix of size n. Here µ t is considered as a nuisance vector parameter since i) it is completely unknown, ii) it is of no interest for the system in charge to detect the anomaly θ and iii) it can mask the anomalies.

Let K α = {φ : sup µ t ∈R q Pr θ=0,µ t (φ(Z t ) = H 1 ) ≤ α} be the class of tests φ : R n → {H 0 , H 1 } with upper-bounded maximum false alarm probability, where the probability Pr θ,µ t stands for the vector of observations Z t being generated by the distribution N (θ + Hµ t , I n ) and α is the prescribed probability of false alarm. The power function β is defined as the probability of detection: β(θ; µ t ) = Pr θ =0,µ t (φ(Z t ) = H 1 ). The subtlety of the above mentioned hypotheses testing problem consists of choosing between H 0 and H 1 with the best possible performance indexes (α, β) while considering µ t as a nuisance parameter.

Anomaly detection methodology

It is easy to see that the problem remains invariant under the group of translations G = {g : g(Z) = Z + Hc , c ∈ R q }. Let S be the family of surfaces S = {S c : c > 0} with S c = {θ : P ⊥ H θ 2 2 = c 2 }. Then, it is shown (Fillatre and Nikiforov (2007a)) that the test

φ * (Z t ) = H 0 if Λ(Z t ) = P ⊥ H Z t 2 2 < λ α H 1 else , (5)
where P ⊥ H = I n -H(H T H) -1 H T and the threshold λ α is chosen to satisfy the false alarm bound α, Pr θ=0,µ t (Λ(Z t ) ≥ λ α ) = α, is Uniformly Best Constantly Powerful (UBCP) 1 in the class K α over the family of surfaces S. The statistics Λ is distributed according to the χ 2 law with nq degrees of freedom. This law is central under H 0 and non-central under H 1 with the non-centrality parameter θ T P ⊥ H θ.

  Fig. 1. Detection of unusual changes in OD traffic volumes. A step-wise anomaly appears at time t0 in the OD flow x(1, 3) and it is routed on links y(1) and y(2).

A test φ * ∈ Kα is UBCP on S if 1) β φ * (θ ′ ) = β φ * (θ ′′ ), ∀θ ′ , θ ′′ ∈ Sc;

2) β φ * (θ) ≥ β φ (θ), ∀θ ∈ Sc, ∀c > 0 for any test φ ∈ Kα which satisfies 1).

Conclusion

The results of numerical experiments with real data will be shown during the presentation of the paper. They confirm the quality of the proposed methods.


Pedro Casas
email: pedro.casas@telecom-bretagne.eu

Lionel Fillatre
email: lionel.fillatre@utt.fr

Sandrine Vaton
email: sandrine.vaton@telecom-bretagne.eu

Multi Hour Robust Routing and Fast Load Change Detection for Traffic Engineering

Traffic Engineering (TE) has become a challenging mechanism for network management and resources optimization due to the uncertainty and the difficulty to predict current traffic patterns. Recent works have proposed robust optimization techniques to cope with uncertain traffic, computing a stable routing configuration that is immune to demand variations within certain uncertainty set. However, using a single routing configuration for long-time periods can be highly inefficient. Even more, the presence of abnormal and malicious traffic has magnified the network operation problem, claiming for solutions which not only deal with traffic uncertainty but also allow to identify faulty traffic. In this paper, we propose two complementary methods to tackle both problems. Based on expected traffic patterns, we adapt the uncertainty set and build a multi-hour yet robust routing scheme that outperforms the stable approach. For the case of anomalous and unexpected traffic, we propose a fast anomaly detection/isolation algorithm which relies on a novel linear spline-based model of traffic demands to identify traffic problems and decide routing changes. This algorithm is optimal in the sense that it minimizes the decision delay for a given mean false alarm rate and false isolation probabilities. Both proposals are validated using real traffic data from two Internet backbone networks.

I. INTRODUCTION

Traffic engineering (TE) represents a major issue for network operators in today's scenario. TE allows the optimization of network resources usage through multiple mechanisms.

In this work, we focus on routing optimization over an Autonomous System (AS). This optimization is becoming increasingly difficult due to the dynamic nature of current traffic. Traffic demands present two different components or behaviors: on one hand, a stable and predictable component due to normal traffic usage patterns (e.g. daily demand fluctuation); on the other hand, an abrupt and unpredictable behavior due to unexpected events, such as network equipment failures, flash crowd occurrences, security threats (e.g. denial of service attacks, virus propagation), external routing changes (e.g. inter-AS routing through BGP) and new spontaneous overlay services (e.g. P2P applications). We use the term volume anomaly [START_REF] Lakhina | Diagnosing Network-Wide Traffic Anomalies[END_REF] to describe these unexpected network events (large and sudden link load changes).

Recent works [START_REF] Ben-Ameur | Routing of Uncertain Traffic Demands[END_REF]- [START_REF] Wang | COPE: Traffic Engineering in Dynamic Networks[END_REF] have proposed a new perspective to the routing optimization under traffic uncertainty: the Robust Routing (RR) approach. In a robust approach of TE, demand uncertainty is taken into account directly within the routing optimization, computing a single routing configuration for all demands within an uncertainty set. While this routing configuration is not optimal for any single traffic matrix (TM) within the set, it minimizes the worst case performance over the whole set. RR provides performance guarantees (i.e. worstcase bounds) for all possible traffic variations within the uncertainty set. However, applying a single robust configuration in the presence of highly variable traffic raises a difficult question: how should this uncertainty set be defined? Larger sets cover a broader group of possible demands, but at the cost of routing inefficiency. On the other hand, tighter sets produce more efficient routing schemes, but subject to poor performance guarantees. RR presents another drawback: it does not solve the problem of faulty traffic identification. The early detection and isolation of unusual and significant changes in traffic demands allows not only to perform an accurate routing reconfiguration but also provides additional information for improving network operation.

A. Related Work

There is a large literature on traffic engineering with uncertain traffic demands. Traditional algorithms rely on a small group of expected TMs (representative traffic demands from past observations) or estimated TMs to compute optimal and reliable routing configurations. An extreme case is presented in [START_REF] Roughan | Traffic Engineering with Estimated Traffic Matrices[END_REF], where routing is optimized for a single estimated TM and it is then applied for long-time periods (24hs periods). Traffic uncertainty is characterized by multiple TMs in [START_REF] Zhang | Optimal Routing with Multiple Traffic Matrices: Tradeoff between Average case and Worst case Performance[END_REF], [START_REF] Zhang | On Optimal Routing with Multiple Traffic Matrices[END_REF] (e.g. set of TMs from previous day, same day of previous week, etc.), and different ways to find optimal routes for the set are presented. Given the dynamic nature of present demands, this perspective is no longer suitable for current scenario [START_REF] Casas | An Adaptive Multi Temporal Approach for Robust Routing[END_REF]. A different approach is provided by online reactive algorithms: TeXCP [START_REF] Kandula | Walking the Tightrope: Responsive yet Stable Traffic Engineering[END_REF] and MATE [START_REF] Elwalid | MATE: MPLS Adaptive Traffic Engineering[END_REF] both balance load in realtime, responding to instantaneous traffic demands. Their main goal is to avoid network congestion by adaptively balancing the load among paths, based on measurement. Reactive routing presents a desirable property, that of keeping routing adapted to current traffic. However, these adaptive algorithms present poor performance under significant and abrupt traffic changes [START_REF] Wang | COPE: Traffic Engineering in Dynamic Networks[END_REF]. A third category of algorithms consists in Stable Robust Routing techniques [START_REF] Ben-Ameur | Routing of Uncertain Traffic Demands[END_REF]- [START_REF] Juva | Robust Load Balancing[END_REF]. In [START_REF] Ben-Ameur | Routing of Uncertain Traffic Demands[END_REF], the authors capture traffic variations by introducing a polyhedral set of demands, ap-plying linear programming techniques to compute an optimal stable routing for all demands within this set. [START_REF] Johansson | Data-driven Traffic Engineering: techniques, experiences and challenges[END_REF] applies this robust technique to compute a robust MPLS routing configuration without depending on TM estimation. Oblivious Routing [START_REF] Applegate | Making Intra-Domain Routing Robust to Changing and Uncertain Traffic Demands: Understanding Fundamental Tradeoffs[END_REF] also defines linear algorithms to optimize worst-case performance for different sizes of traffic uncertainty sets, aiming to handle dynamic changes. [START_REF] Juva | Robust Load Balancing[END_REF] analyses the use of robust routing through a combination of traffic matrix estimation and its corresponding estimation error bounds, in order to shrink the uncertainty set. The drawback of stable robust routing is its inherent dependence on the definition of the uncertainty set: larger sets allow to handle a broader group of traffic demands, but at the cost of routing inefficiency; conversely, tighter sets produce more efficient routing schemes, but subject to poor performance guarantees.

As regards anomaly detection in data networks, the problem has been extensively studied. In this section, we will just overview those works that have motivated our signalprocessing based detection algorithm. Signal processing techniques have been applied to the anomaly detection field [START_REF] Barford | A Signal Analysis of Network Traffic Anomalies[END_REF]- [START_REF] Cheng | Use of Spectral Analysis in Defense Against DoS Attacks[END_REF]. The usual behavior of data flows is modeled by several approaches: spectral analysis, time series analysis, wavelets decomposition, etc. Anomalies correspond to deviations from the usual behavior of the data flows. The general flaw of these algorithms is the lack of stability over time of the proposed traffic models, as well as the absence of optimality conditions for the detection in most cases. A second class of methods related to our model concerns statistical hypotheses testing [START_REF] Lakhina | Diagnosing Network-Wide Traffic Anomalies[END_REF]- [START_REF] Thottan | Anomaly Detection in IP Networks[END_REF]. When data flows are parametrically modeled, the design of optimal algorithms is possible. Nevertheless, nonparametric approaches are particularly studied because of the lack of parametric models, and these approaches are often suboptimal. The detection/isolation of traffic anomalies problem was previously treated in [START_REF] Lakhina | Diagnosing Network-Wide Traffic Anomalies[END_REF], using a TM decomposition on the Principal Component Analysis (PCA) basis. However, this approach presents a major stability problem: the PCA basis depends on the measurement period, rendering it unstable over time.

B. Contributions of the Paper

We propose two novel and complementary approaches to deal with current dynamic traffic demands, separately treating both traffic uncertainty sources. For expected traffic fluctuations, we present a time varying approach of RR that outperforms the current stable approach: the Multi-Hour Robust Routing (MHRR). We preserve the virtues of RR, but change the routing configuration during time. The uncertainty set is optimally divided into several uncertainty sub-sets that better adapt to real traffic loads and a stable robust routing scheme is computed for each sub-set. The partitioning algorithm allows to calculate the exact moments when routing changes must be performed. For the case of unpredictable traffic behavior, we propose a novel volume anomaly detection/isolation algorithm to identify traffic problems and decide routing changes. To overcome the limitations of the PCA approach [START_REF] Lakhina | Diagnosing Network-Wide Traffic Anomalies[END_REF], we propose a non data-driven traffic model which remains sta-ble over time. This model allows to separate normal from anomalous traffic, based on simple link load measurements. Both proposals are validated using real traffic data from two backbone networks, the Internet2 Abilene backbone network and a private international Tier-2 network.

The remainder of this paper is organized as follows. In Section II, we recall the basic aspects of the robust routing approach. Section III presents the theoretical background and empirical evaluation of the MHRR. The proposed algorithm and traffic model for anomaly detection/isolation are introduced and validated in section IV. Finally, Section V concludes this work.

II. ROBUST ROUTING

Let us consider a network topology defined by a set of n nodes and L = {1, . . . , r} links with capacities in C = (c 1 , c 2 , . . . , c r ). The TM demand d = {d i,j } denotes the traffic flow between every node i and node j (i = j) of the network. We re-arrange d as a column vector, d = {d k, k=1..m }, where d k represents the traffic flow transmitted by OD pair k (OD-flow k) and m = n × (n -1) is the number of OD pairs. Let N = {OD 1 , . . . , OD m } be the set of OD pairs. Link's information y l represents the total traffic (i.e. aggregated OD flows) through link l in a certain period of time. This information is available from router's MIB variables and it is usually collected every 5' periods via SNMP [START_REF] Feldmann | Deriving Traffic Demands for Operational IP Networks[END_REF]. Traffic demands and links' traffic are related through the routing matrix R, a r × m matrix R = {r l,k } where 0 r l,k 1 represents the fraction of OD demand k routed through link l:

y = R × d. (1)
with y = {y l, l=1..r }. Routing optimization depends on the underlying data transport mechanism; we will focus on path-based routing such as MPLS. This optimization consists of minimizing certain performance metric associated with traffic demand. Throughout this work we consider maximum link utilization (MLU) as the routing performance criterion. Overloaded links tend to cause QoS degradation (e.g. larger delays and packet losses, throughput reduction, etc.), so MLU represents a reasonable measure of network performance. The are many other performance metrics that could be used instead of MLU, like path's end-to-end delay or mean link utilization; setting the focus too strictly on the MLU can often lead to longer average traffic paths and thus adversely affect the mean. However, we use the MLU as it is the most commonly applied criterion and it represents an easy to understand performance metric. For a given routing matrix R = {r l,k } and a traffic demand d, the MLU (u max ) is defined as the maximum of the ratio between link load and link capacity:

u max (C, d, R) = max l∈{1...r} k r l,k • d k c l = max l∈{1...r} y l c l (2)
Let P (k) be the set of possible paths for OD demand k. Let x k p be the proportion of traffic demand d k that flows through path p ∈ P (k), 0

x k p 1. Finally, let x k l be the proportion of traffic demand d k that flows through link l ∈ L, 0

x k l 1. We define D as the uncertainty set where traffic demand may vary. This set can be defined in different ways, depending on the available information: link load measurements and historical routing, a set of previously observed TMs {d (3)

p∈P (k) x k p 1 ∀ k ∈ N p∈P (k), l∈p x k p x k l ∀ k ∈ N, ∀ l ∈ L k∈N x k l .d k umax • c l ∀ l ∈ L, ∀ d ∈ D x k p , x k l 0 ∀ l ∈ L, ∀ p ∈ P (k), ∀ k ∈ N umax 1
programming techniques, applying a combined column and constraint generation method [START_REF] Ben-Ameur | Routing of Uncertain Traffic Demands[END_REF]. In a traditional robust routing application, the obtained routing configuration is applied during long-term periods of time (i.e. daily routing). In this sense, we refer to robust routing as Stable Robust Routing (SRR).

III. MULTI-HOUR ROBUST ROUTING

In [START_REF] Casas | An Adaptive Multi Temporal Approach for Robust Routing[END_REF] we present the advantages of the SRR with respect to traditional routing approaches: SRR offers stability guarantees against traffic uncertainty and traffic timevariations at a reasonable cost. However, considering a single routing scheme for long-time periods is conservative and results in sub-optimal performance. We propose a simple approach to shrink and adapt the uncertainty set along time that outperforms the SRR. Based on rough knowledge of traffic variations (i.e. considering expected traffic behavior), we propose to optimally divide the uncertainty set and build a multi-hour routing configuration, considering a single SRR configuration for each sub-set. Daily traffic changes can be seen as a time variation of the uncertainty set. At each time t, the routing matrix R and the link load values y(t) = y t define an instantaneous uncertainty set

D(t) = {d ∈ R m , R × d y t , d 0}.
The continuous union of infinite instantaneous uncertainty sets along time t defines the daily uncertainty set 1(a) explains this idea. Assuming this set is an union of polytopes, [START_REF] Ben-Ameur | Between Fully Dynamic Routing and Robust Stable Routing[END_REF] provides a theoretical study of the optimal partitioning of D t , using a partitioning hyper plane. [START_REF] Ben-Ameur | Between Fully Dynamic Routing and Robust Stable Routing[END_REF] proves that this is a NP-hard problem, except for the case where a partitioning direction is previously fixed. We define a partitioning hyper plane by its direction vector α and a value w: α.d = w. In the MHRR approach, we consider a particular direction for partitioning: the time direction. In that case, w represents the time of the day. We define h + 1 hyperplanes at times {w 1 , w 2 , .., w h+1 }. The intersection between D t and the half-spaces defined by these partitioning hyperplanes results in h uncertainty sub-sets

D t = (d, t) ∈ R m+1 , d ∈ ∪ t1 t tτ D(t), t 1 t t τ . Figure
D i = {D t ∩ {d, α.d w i } ∩ {d, α.d w i+1 }} , ∀i = 1, ., h.
Let D i be the smallest single-time set that contains all demands d(t) ∈ D i , w i t w i+1 (see figure 1(b)). A SRR configuration R i robust is computed for each sub-set D i . Each routing configuration is finally applied at each time interval. The optimal values of routing changes w * = {w * 2 , . . . , w * h } are the solution for the following optimization problem (w 1 and w h+1 are fixed a priori, as they define the considered time interval of analysis):

w * (D t ) = arg min w max i=1..h u max (D i ) (4)
where u max (D i ) is the solution for (3) for polytope D i . [START_REF] Ben-Ameur | Between Fully Dynamic Routing and Robust Stable Routing[END_REF] presents a simple algorithm to approximately solve (4) (within an arbitrary precision), using a generalization of a simple dichotomy methodology. The MHRR presents a trade-off between performance and routing stability. The more intervals we use, the more adapted the routing becomes. However, the number of intervals should be bounded as many routing changes may lead to instabilities and performance degradation. In a general case, 2 sub-sets are enough to handle the usual daily variation.

MHRR Evaluation

We present a comparative analysis between SRR and MHRR in Abilene, an Internet2 backbone network. Abilene consists of 12 router-level nodes and 30 OC192 links (2 OC48). The used router-level network topology and traffic demands are available at [START_REF] Zhang | Abilene Dataset 04[END_REF]. Traffic data consists in 6month traffic matrices collected every 5' via Netflow from the Abilene Observatory [START_REF]The Abilene Observatory[END_REF]. The time-variation of the polytope is not a simple homothety [START_REF] Casas | An Adaptive Multi Temporal Approach for Robust Routing[END_REF]; in this sense, we will show that a routing configuration change during the day improves routing performance. Let R o be the historical routing matrix of Abilene, not necessarily optimal (R o is available at [START_REF] Zhang | Abilene Dataset 04[END_REF]). We consider a single time partitioning (i.e. 2 routing intervals), w 1 = 20:00, w 2 = w * and w 3 = 21:00, where w * is the solution for (4). For each time interval, we consider the smallest polytope that includes all possible realizations over that period:

D A,B = {d ∈ R m , R o × d y A,B , d 0} (5)
where y A = y 20:00-w * max and y B = y w * -21:00 max (maximum values for each link). In this way, D A includes all traffic demands between 20:00 and w * and D B between w * and 21:00 (see figure 1(b)). For each polytope, we compute a SRR configuration, R A robust and R B robust . In order to compare stable and multi-hour approaches, we apply both routing configurations during the whole evaluation period. We include the routing performance obtained with R o (curve historical routing) to appreciate the time variation of traffic loads. Figure 22 2(a) compares the routing performance (MLU) between these two RR configurations. Polytope D A is well suited for smaller loads, so R A robust performs better during the first half of the day, when network load is lower. However, when traffic increases, demands that do not belong to D A produce higher link utilizations than those obtained with R B robust . The MHRR consists in computing the time when routing must be changed (w * ≈ 8:00 in this case), using the corresponding routing configuration depending on the time of the day (R A robust before w * and R B robust after). The MHRR approach presents a performance improvement of 15% with respect to the SRR approach before w * , reaching a near 20% of over-efficiency after w * . We repeat the same evaluation but considering a traffic demand that drastically changes (i.e. a large timevariation of the polytope, caused by a volume anomaly). Figure 2(b) presents an abrupt change in MLU (almost 14 times higher) at time 18:00. In this case, we assume that this change is known in advance (note that in the general case, it is not possible to predict these abrupt changes). The optimal moment for changing routing is w * ≈ 18:00. The MHRR approach definitely outperforms the SRR in this experience, presenting a MLU between 10% and 60% smaller during the whole evaluation period.

IV. DEALING WITH UNEXPECTED EVENTS

The proposed MHRR approach offers a robust and efficient routing configuration, given a rough knowledge of the daily uncertainty set. However, in the presence of volume anomalies it is no longer possible to apply the MHRR as the daily uncertainty set is unknown. For those cases, we propose a fast volume anomaly detection/isolation algorithm to quickly identify faulty traffic. This detection allows to decide as soon as possible the moment when routing configuration must be changed. The goal of the algorithm is to detect/isolate an additive change θ in the time series of traffic demand d(t) from a sequence of link load measurements y(t) = R × d(t). We use link loads as input to avoid relying on seldom available traffic demands. In this work, we focus on detecting and isolating a "localized" anomaly, θ = θ (δ 1,i , . . . , δ i,i , . . . , δ m,i )

T , where δ i,j = 0 if i = j and δ i,i = 1 (this corresponds to a change θ in OD flow i). If several OD flows are simultaneously corrupted, the detection/isolation algorithm produces an alarm and identifies only one faulty OD flow. The algorithm can be extended to detect/isolate simultaneous anomalies, but its' complexity (n o hypotheses, see IV-B) grows highly. The isolation of the anomalous traffic is possible since an anomaly in a given OD flow typically spans multiple links. Real traffic demands follow a non-observable model from link load measurements: since r < m, it is impossible to retrieve d(t) from y(t) without additional assumptions on the traffic demand. To overcome this difficulty, we propose a parsimonious linear model for non-anomalous traffic. This model renders traffic demands observable and therefore, it allows to separate usual from anomalous traffic.

A. Stochastic Traffic Model for Anomaly Detection

We assume that the stochastic process of the OD traffic demand d(t) obeys the following linear model:

d(t) = λ(t) + ξ(t) (6)
where λ(t) ∈ R m is the mean traffic demand and ξ(t) is a white Gaussian noise with covariance matrix Σ(t) that represents the model error together with the natural variability of the OD flows (based on the results obtained in [START_REF] Coates | Internet tomography[END_REF]). The process λ(t) represents the "regular" part of the OD TM which can be correctly modeled when the behavior of the network is anomaly-free. We propose to parameterize this vector by exploiting the stationarity of the spatial distribution of the TM. One of the few invariants of Internet traffic is that a small percentage of flows contribute to a large proportion of total traffic [START_REF] Johansson | Data-driven Traffic Engineering: techniques, experiences and challenges[END_REF], [START_REF] Medina | Traffic Matrix Estimation: Existing Techniques and New Directions[END_REF]. Hence, if we assume that the traffic distribution between the different OD couples is spatially stationary in the absence of an anomaly, the order of increasing OD flows remains constant during long time periods. The proposed traffic model takes advantage of the stationary property of this ordering. We propose to classify OD flows in three different classes, depending on their volume: large OD flows, small OD flows and medium-size OD flows. The sorted components can be interpreted as a discrete increasing signal. The curve obtained by interpolating this discrete signal is assumed to be a continuous curve, hence it can be parameterized by using a polynomial approximation. Figure 3 shows the OD flows, sorted in the increasing order of their volume of traffic, as a function of the time t. Since data are vectors of finite dimension, we propose to use the following method to design a discrete spline basis: (i) we choose a continuous spline basis; (ii) we discretize all these splines according to m points uniformly chosen in the interval [1; m] and (iii) we rearrange all these discrete signals according to previous sorting order.

We finally obtain the following linear model for the anomalyfree traffic demand:

d(t) = Sµ(t) + ξ(t) (7)
where S = (s 1 s 2 . . . s q ) is a m × q known matrix with columns s j and q is small with respect to m. The vectors s i , which correspond to the rearranged discrete spline, form a set of known basis vectors describing the spatial distribution of the traffic and µ(t) = (µ 1 (t) . . . µ q (t)) T is the unknown time varying parameter vector which describes the OD flow intensity distribution with respect to the set of vectors s i . The model for the anomaly-free link traffic is given by:

y(t) = Hµ(t) + ζ(t), (8)
where H = RS and ζ(t) = Rξ(t). In this way, we can describe the usual behavior of traffic demands from simple link measurements. The computation of the rank of H is not simple since it depends on the routing matrix R. In practice, since the number of columns of H is very small, the product RS and its rank can be computed very fast. Therefore, we will assume that H is full column rank. Finally, the covariance matrix Σ is unknown. The remedy consists in computing an estimate Σ of Σ. Results on the estimation of Σ can be found in [START_REF] Vardi | Network Tomography: Estimating Source-Destination Traffic Intensities from Link Data[END_REF].

B. Volume Anomaly Detection/Isolation

The detection/isolation of a volume anomaly at time t 0 can be treated as a hypothesis testing problem where the null hypothesis H 0 t0 = {the OD flows are anomaly-free at time t 0 } is tested against m alternatives H j t0 = {the j-th OD flow presents an anomalous additional amount of traffic θ from the time t 0 }. The change detection algorithm has to compute a pair (T, ν), where T is the alarm time at which a ν-type change (ν ∈ {1, 2, . . . , m}) is detected and isolated, based on link traffic observations y 1 , y 2 , ... The hypothesis testing can be written as

H 0 : y(t) ∼ N (H µ(t), RΣR T ) , t = 1, 2, . . . , (9)
H j t0 :    y(t) ∼ N (H µ(t), RΣR T ), t = 1, . . . , t 0 -1, y(t) ∼ N (H µ(t) + θ j r j , RΣR T ) , θ j,1 |θ j | θ j,2 , t = t 0 , t 0 + 1, . . . (10)
where r j, 1 j m denotes the normalized j-th column of R and 0 < θ j,1 < θ j,2 < +∞ are some known bounds on the change intensity of the j-th OD flow (these bounds are introduced for technical reasons but they can be chosen arbitrarily). As we show in the Appendix, we can simplify this problem by eliminating the non-anomalous traffic. In this case, hypothesis [START_REF] Brutlag | Aberrant Behavior Detection in Time Series for Network Monitoring[END_REF] can be rewritten as

H j t0 :
   z(t) ∼ N (0, I r-q ) , t = 1, . . . , t 0 -1, z(t) ∼ N (θ j v j , I r-q ) , θ j,0 |θ j | θ j,1 , t = t 0 , t 0 + 1, . . . [START_REF] Cheng | Use of Spectral Analysis in Defense Against DoS Attacks[END_REF] where v j is a known vector and z(t) are the normalized residuals obtained from y(t) after filtering the non-anomalous traffic. The vector v j corresponds to the signature in the residuals of a change in OD flow j. We use the optimal recursive algorithm (T r , ν r ) proposed in [START_REF] Nikiforov | A lower bound for the detection/isolation delay in a class of sequential tests[END_REF] to solve [START_REF] Cheng | Use of Spectral Analysis in Defense Against DoS Attacks[END_REF] :

T r = min 1 k m {T r (k)}, ν r = arg min 1 k m {T r (k)} T r (k) = inf t 1 : min 0 j =k m [g t (k, j) -h k,j ] 0 (12)
with g t (k, j) = g t (k, 0)g t (j, 0). The recursive functions g t (k, 0) are defined by

g t (k, 0) = (g t-1 (k, 0) + z t (k, 0)) + (13)
z t (k, 0) = log f k (z(t)) f 0 (z(t)) (14)
g 0 (k, 0) = 0 for every 1 k m and g t (0, 0) = 0 for all t. f 0 represents the probability density function of anomaly-free traffic measurements. f k is the probability density function of residuals z(t 0 ), z(t 0 + 1), .. after a change of type k. The thresholds h k,j are chosen by the following formula:

h k,j = h d if 1 k m and j = 0 h i if 1 k, j m and j = k
where h d is the detection threshold and h i is the isolation threshold. For given bounds γ and β, this algorithm is asymptotically optimal, i.e. it reaches the lower bound of the maximum mean delay for detection [START_REF] Nikiforov | A lower bound for the detection/isolation delay in a class of sequential tests[END_REF]. The choice of the detection and isolation thresholds h d and h i is discussed (with practical comments and simulations) in [START_REF] Nikiforov | A simple recursive algorithm for diagnosis of abrupt changes in random signals[END_REF].

C. Validation

We demonstrate the ability of the detection/isolation algorithm to detect and identify a volume anomaly in SNMP link flow data from two different networks (different not only in the topology but also in the behavior of traffic demands): a large Tier-2 network (50 nodes, 168 measured links and 2450 non-zero OD flows, sampled at a 10 minute rate) and Abilene (the Abilene dataset consists in Netflow traces, so we use the supplied routing matrix to retrieve link loads). Figure 4 shows the typical realizations of the decision functions g t (i, 0) and s t (i) = min 0 i =k m [g t (i, k)h i,k ] vs the elapsed time. The functions s t (i) are used to "monitor" the OD flows; when the function s t (i) exceeds 0, OD flow i is declared faulty. It is assumed that the anomaly in the Tier-2 network begins at time 3660, and at time 1070 in Abilene. Note that after this time, several decision functions g t (i, 0) rapidly grow. Each function g t (i, 0) is associated with OD flow i and when this function grows, it means that OD flow i is suspected of carrying an abnormal amount of traffic. Contrary to g t (i, 0), only decision function s t (159) (s t (87) in Abilene) associated to faulty OD flow 159 (87 respectively) grows and finally exceeds the threshold. Hence, the functions s t (i) permit us to isolate the faulty OD flow among all the OD flows associated to functions g t (i, 0) that have rapidly grown. At time 3660 (1070 respectively), an alarm is raised and the algorithm selects the faulty OD flow 159 (87 respectively). The decision function s t (i) needs only 1 observation (10 minutes in the Tier-2 network or 5 minutes in Abilene, but this is the smallest delay than can be achieved given these sampling-rates) to detect and isolate the faulty OD flow. An interesting observation is that the detection/isolation algorithm achieves good results in both networks, even though the respective traffic demand behaviors are completely different between these two networks.

V. CONCLUSIONS AND FUTURE WORK

In this paper, we address the routing under traffic uncertainty problem. We provide a solution that not only deals with current dynamic traffic demands in a robust and efficient way but also detects and isolates large-volume anomalous traffic, improving network operation. We extend the robust routing paradigm by introducing the notion of time-varying uncertainty set, setting up a multi-hour robust routing scheme. We show that this approach achieves better resource utilization than previous stable robust proposals in different scenarios. We introduce an original linear spline-based parsimonious model to parameterize normal traffic behavior from widely available link load measurements. Based on this model, we present a statistical algorithm to detect and isolate volume anomalies in network traffic. We apply this algorithm to cope with sudden and large traffic changes in current dynamic demands, complementing the multi-hour robust scheme. Many important issues remain open for further study. A deep evaluation of the challenges involved in changing routing configuration after the detection of the anomalous traffic should be performed. The isolation ability of the proposed algorithm can be highly exploited in order to recompute routing in a more efficient way. The impact of routing re-configuration on end-to-end traffic must be explored, especially considering the imposed QoS restrictions in the current end-user Internet-services scenario.

Fig. 1 .

 1 Fig. 1. (a) Daily variation of the polytope Dt, (b) time partitioning of Dt.

Fig. 2 .

 2 Fig. 2. Routing performance, stable vs. multi-hour robust routing.

Fig. 3 .

 3 Fig. 3. Approximation of real OD flows by the spline-based model.

Fig. 4 .

 4 Fig. 4. Typical realizations of decision functions for a Tier-2 network (a,b) and Abilene (c,d).

  1 , d 2 , . . . , d o }, TM time series d(t), etc.[START_REF] Ben-Ameur | Routing of Uncertain Traffic Demands[END_REF] defines this set as a polytope, based on the intersection of several half-spaces that result from linear constraints imposed to traffic demand. The Robust Routing Optimization Problem (RROP) consists of minimizing u max , considering all demands within D (3). The RROP can be efficiently solved by linear

	minimize	umax
	subject to:

APPENDIX -ELIMINATION OF NON-ANOMALOUS TRAFFIC

Non-anomalous traffic Hµ(t) is eliminated by projecting the measurement vector y(t) on the null space of H. By using the invariant properties of the Gaussian law, the general covariance matrix in [START_REF] Brutlag | Aberrant Behavior Detection in Time Series for Network Monitoring[END_REF] is reduced to the identity one. Let us define the matrix W = (w1, .., wr-q) of size r × (rq) composed of eigenvectors w1, .., wr-q of the projection matrix P ⊥ H = Ir -H(H T H) -1 H T corresponding to eigenvalue 1. The matrix W satisfies the following conditions: W T H = 0, W W T = P ⊥ H and W T W = Ir-q. The matrix W can be considered as a linear rejector that eliminates the non-anomalous traffic. Under hypothesis H j t 0 , the sequence W T y(t) can be modeled as W T y(t) = W T ζ(t) + θj W T rj , j = 1, .., m. Since W T ζ(t) is a correlated Gaussian vector with covariance matrix Σ = W T RΣR T W , each vector W T y(t) is normalized by using the square root matrix Σ 1 2 , z(t) = Σ -1 2 W T y(t) ∼ N ( θj vj, Ir-q), with vj = Σ -1 2 W T rj .