Pedro Casas
email: pedro.casas@enst-bretagne.fr

Lionel Fillatre
email: lionel.fillatre@utt.fr

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

Robust and Reactive Traffic Engineering for Dynamic Traffic Demands

Keywords: Traffic Uncertainty, Multi-Hour Robust Routing, Anomaly Detection/Isolation, Reactive Robust Routing

Traffic Engineering (TE) has become a challenging mechanism for network management and resources optimization due to uncertain and difficult to predict 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 longtime 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 detect and identify faulty traffic to take the appropriate countermeasures. In this paper, we introduce the Reactive Robust Routing (RRR) for TE, an approach that combines both proactive and reactive techniques to tackle the problem. Based on expected traffic patterns, we adapt the uncertainty set and build a multi-hour yet robust routing scheme that outperforms the stable robust approach. For the case of anomalous and unexpected traffic, we propose a fast anomaly detection/isolation algorithm to detect and localize abrupt changes in traffic flows 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 probability. We validate these proposals using real data from two different backbone networks and we show how the RRR can handle uncertain and highly dynamic traffic in an automatic fashion, simplifying network operation.

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 usual 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 crowds 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 fashion 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. In this sense, RR provides performance guarantees (i.e. worst-case bounds) for all possible traffic variations within the uncertainty set. The RR approach can be used as a proactive technique to deal with dynamic traffic. It can handle changing demands at a reasonable cost (with respect to an ideal but illusive optimal adaptive routing) up to a certain limit (given by the size of the uncertainty set). However, applying a RR algorithm to address both traffic behaviors (usual traffic as well as volume anomalies) is an inefficient strategy: a single routing can not be suitable for both situations.

On the contrary, a reactive approach could be used as a complementary strategy to enhance RR performance, responding to abrupt and large traffic changes with an effective routing reconfiguration. Volume anomalies may have an important impact on the network performance, causing sudden situations of strong network congestion. The early detection and isolation of these anomalies allows to modify the routing as soon as possible, limiting their impact. In this work, we propose a signal processing algorithm for fast load change detection/isolation. Through out the paper, we use the term anomaly isolation to refer to the identification and localization of an anomalous flow among the network traffic.

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, applying 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, and discusses corresponding methods for robust OSPF optimization. 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. In [START_REF] Wang | COPE: Traffic Engineering in Dynamic Networks[END_REF], the authors introduce COPE, an approach to deal with this tradeoff in the size of the uncertainty set, combining traditional algorithms with oblivious routing. COPE optimizes routing for predicted demands and bounds worst-case performance to ensure acceptable efficiency under unexpected traffic events. Nevertheless, COPE proposes a long-term stable routing configuration as previous works do (24hs periods), losing the adaptability (and hence the performance efficiency) of reactive routing. Besides, it is possible not only to assure performance guarantees for unexpected events, but to obtain optimal routings for this traffic.

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 sub-optimal. 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

In the final remarks of [START_REF] Johansson | Data-driven Traffic Engineering: techniques, experiences and challenges[END_REF], the authors raised an interesting reflection: "it is not clear whether time-varying demands should be addressed using proactive (e.g. robust routing) or reactive (dynamic, adaptive) methods". In this work, we propose to use both proactive and reactive 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). The stable RR may be costly. However, it is easy to control its cost by shrinking the uncertainty set. 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 optimally calculate the exact times 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. This algorithm allows both to detect the volume anomaly and to identify and localize the anomalous flow. To overcome the stability problems of previous approaches, we propose a non data-driven traffic model which remains stable over time. The main contribution of this detection algorithm relies on the well established conditions of optimality that it presents, a feature generally absent in previous works.

Both proactive and reactive methods are combined into a novel TE approach for dynamic traffic demands: the Reactive Robust Routing (RRR). This approach uses the MHRR to handle typical changes in traffic demands and the detection/isolation algorithm to deal with unexpected volume anomalies. The RRR exploits the isolation ability of the detection/isolation algorithm to compute an adapted robust routing configuration after the anomalous traffic detection, reducing its impact on network performance during its prevalence. In addition, it also provides a simple yet effective method to automatically detect the end of the anomaly, returning to the MHRR configuration. Contrary to previous works in the field, our proposal optimizes routing in a robust and adaptive fashion for every possible traffic demand (and not only for the commoncase traffic). A key feature of the RRR approach relies on the fact that the whole routing configuration/reconfiguration algorithm is completely automatic, an interesting property that simplifies network operation by self-managing. All the proposed algorithms in this work 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. Section V presents the Reactive Robust Routing, showing the automatic interaction between the proactive and the reactive components through complete real and simulated examples. Finally, Section VI 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 which element 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 in 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. 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 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 in minimizing u max , considering all demands within D (3). The RROP can be efficiently solved by linear 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 routminimize umax subject to:

(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
ing 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 continious 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 the usual daily variation.

MHRR Evaluation

We present a comparative analysis between SRR and MHRR in Abilene, an Internet2 backbone network. Abilene consists in 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 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 loading 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, provided 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" anomaly1 , θ = θ (δ 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). 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. The process λ(t) represents the "regular" part of the OD TM which can be correctly modeled when the behaviour 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 mediumsize 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 on 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 anomaly-free 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 , ... This algorithm is optimal in the sense that it minimizes the maximum mean delay for detection/isolation, for a given minimum mean time before a false alarm γ and maximum false isolation probability β, both defined by an expert user. 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 (10) 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, . . . ( 11
)
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] Roughan | Traffic Engineering with Estimated Traffic Matrices[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 of this evaluation 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. REACTIVE ROBUST ROUTING

Both proactive and reactive methods (the MHRR and the anomaly detection/isolation algorithm respectively) are combined into a single approach we refer to as the Reactive Robust Routing (RRR). This approach provides an automatic method for robust routing configuration/reconfiguration, based on the monitoring of the network state. The RRR exploits the isolation ability of the detection/isolation algorithm to compute a new robust routing configuration after the detection of an anomalous OD flow; at the same time, it detects the end of the anomaly (if there is any) and returns to the usual MHRR routing.

A. Routing Reconfiguration

We propose a simple method that exploits both the RR approach and the isolation ability of previous detection/isolation algorithm to compute the new routing scheme to apply after the detection step. The idea of this reconfiguration is to minimize the impacts of the detected anomaly on the network performance. We assume that before the detection of the anomalous traffic, a stable RR configuration R A robust is applied, computed on the basis of the historical routing R o and the link load y o that results from the MHRR algorithm (5) (R A robust is obtained from (3), using

D A = {d ∈ R p , R o × d y o , d 0}).
C = {d * ∈ R p , R o × d * y o + R o θ, d 0}.
The reader should bear in mind that the type of anomalies we deal with are generally caused by external factors (e.g. external routing changes, flash crowds, denial of service attacks); this justifies the relevance of the uncertainty set expansion with respect to R o . Said in other words, we detect and identify which is the anomalous OD flow, an then we consider a bigger uncertainty set that takes into account the abrupt change of this OD flow. The new RR scheme R C robust is the solution for (3), using D C . To avoid the estimation of the unknown anomalous volume θ, we can expand D A to the limits of links' capacities, in the direction of OD flow k :

D C = d * ∈ R p , R o × d * y o .ζ i,k + C.ζ i,k , d 0 , where ζ i,k = 1 if OD flow k is routed across link i = 1.
.r and ζ i,k = 0 otherwise, and † denotes the ones' complement of † (bitwise NOT, i.e. 0 → 1 and 1 → 0). While the outcome of this approach may result in routing inefficiency, it avoids the estimation errors of θ (i.e. we build a more robust routing). Figure 5 presents the evaluation of the RRR approach in the presence of a sudden and abrupt load change. We consider the same situation of figure 2(b), comparing the routing performance of the MHRR and the RRR respectively. As in section III, we assume the daily uncertainty set is completely known for the case of the MHRR (i.e. the abrupt change is known in advance). For RRR, the anomaly is automatically detected and the new routing configuration is computed and immediately applied, based on the expansion of the uncertainty set. We can appreciate that the routing performance of the RRR is slightly worst than the one obtained with the MHRR (less than 2%). Nevertheless, the RRR represents a real scenario, where the anomaly can not be forecasted and has to be detected to compute an accurate rerouting.

B. Back to the MHRR scenario

After the anomaly detection and the robust routing reconfiguration, we must provide a way to detect the end of the anomaly, in order to return to the MHRR situation. This detection can be easily achieved by using a simplified version of our detection algorithm: suppose that we detect and isolate a k type anomaly at time t 0 (i.e. OD flow k is declared as anomalous). For every time t > t 0 , we only monitor OD flow k until no anomaly alarms are raised, showing the end of this anomaly (remember that in this work we have only considered "localized" anomalies, i.e. anomalies in a single OD flow at a time). As we focus on a single OD flow, the multi-hypotheses test ( 9), ( 10) becomes a single hypothesis test, where the null hypothesis H 0 t = {the OD flow k is anomaly-free at time t} is tested against H 1 t = {the k-th OD flow presents an anomalous additional amount of traffic at time t} (note that this hypothesis test is slightly modified, as both hypothesis are exchanged). Instead of using a sequential probability test (e.g. single CUSUM algorithm [START_REF] Basseville | Detection of abrupt changes: theory and applications[END_REF]), we propose to apply a simple LLR test (Log Likelihood Ratio Test) at each time t. For this purpose, we consider the decision function τ t (k)

τ t (k) = log f k (z(t)) f 0 (z(t)) + -h d (15)
In the presence of an anomalous behavior in OD flow k, τ t (k) > 0 and the anomaly alarm keeps raising. When τ t (k) becomes negative, the anomaly alarm ceases, pointing out its end. To conclude this section, we present in figure 6 an evaluation of the complete RRR approach under the presence of a volume attack (e.g. single DoS attack). We introduce an artificial sudden and large volume change in OD flow 63 of the Abilene dataset. This artificial traffic is put on top of the usual daily traffic between times 1125 and 1350. The first step of the RRR consists in computing the MHRR, using an expected daily uncertainty set. The optimal division (4) results in w * = 1230. The evaluation begins at time 1020, when the MHRR decides to apply the SRR R A robust (SRR A in fig. 6). The detection/isolation algorithm continuously monitor the network state, and at time t 0 = 1125 detects and localizes an anomalous behavior in OD flow 63 (figure 6(a)). After the detection (and before the new sampling of link loads y(t 0 + 1), i.e. a 5' time-window) the new routing configuration is computed, according to V-A. At time t = t 0 + 1 the new routing configuration is deployed and the anomaly-end detection phase begins. It is important to note that the matrix H = RS as well as the usual-traffic rejector (see the Appendix) must be recomputed after the change of the routing matrix R (in fact, the same re-computation must be conducted every time the routing matrix changes, restarting the detection algorithm to avoid transient effects). The decision function τ t (63) remains positive for every time t > t 0 , until time t = 1350, when the negative value of τ t (63) shows the end of the anomalous behavior in OD flow 63. At this time, the RRR compares t with w * in order to decide which routing to apply, whether SRR R A robust or R B robust (R A robust if t < w * or R B robust if t > w * ). Once the new routing configuration is established, the anomaly detection/isolation algorithm starts again to search for anomalous behaviors. The performance improvements of the RRR are evident, up to a 40% wrt the MHRR and near 50% wrt the traditional SRR approach.

VI. CONCLUSIONS AND FUTURE WORK

In this paper, we address the routing optimization 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 usual traffic behavior from widely available link load measurements. Compared to many other traffic models, ours remains stable along time, a necessary condition to achieve reliable results. Based on this model, we present a statistical algorithm to detect and isolate volume anomalies in network traffic. This algorithm presents wellestablished conditions of optimality, unavailable in previous proposals in the field. We apply this algorithm to cope with sudden and large traffic changes in current dynamic demands, complementing the multi-hour robust scheme. We propose a simple method that exploits both the RR approach and the isolation ability of previous detection/isolation algorithm to compute the new routing scheme to apply after the detection step. The idea of this reconfiguration is to minimize the impacts of the detected anomaly on the network performance. All these algorithms are merged into a new proposal for robust and reactive routing optimization, the Reactive Robust Routing. The RRR approach deals with traffic uncertainty in a completely automatic fashion, simplifying network management.

We believe that the RRR represents a first step towards a dynamic and robust routing policy, but many important issues remain open for further study. A deep evaluation of the impact of routing re-configuration on end-to-end traffic must be conducted, especially considering the imposed QoS restrictions in the actual 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).

  After the detection and isolation of the faulty OD flow k, the anomalous-free traffic demand d takes the value d * = d + θ, with θ = θ.δ k , where δ k = (δ 1,k , . . . , δ k,k , . . . , δ p,k ) T , δ i,k = 0 if i = k and δ k,k = 1. We can expand the uncertainty set D A in the directions of the routed OD flow k (with respect to R o , i.e. the routing configuration that defined D A ), obtaining an expanded uncertainty set D

Fig. 5 .

 5 Fig. 5. Reactive Robust Routing -routing reconfiguration after the detection af a large and abrupt traffic change.

  End of the ktype Anomaly Level of detection τt(OD flow k=63) (a) Anomaly detected in OD flow 63. (b) Anomaly-end detected.

Fig. 6 .

 6 Fig. 6. Reactive Robust Routing performance under a simulated DoS attack.

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 simultaneous anomalies, but the complexity (n o operations) grows like m 4 , where m is the n o of OD flows.

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] Ben-Ameur | Between Fully Dynamic Routing and Robust Stable Routing[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.


Lionel Fillatre

Igor Nikiforov

Pedro Casas

Sandrine Vaton

OPTIMAL VOLUME ANOMALY DETECTION IN NETWORK TRAFFIC FLOWS

Optimal detection 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 simple link load measurements, the anomaly detection is an ill-posed decision-making problem. The method proposed in this paper consists of finding a linear parsimonious model of ambient traffic (nuisance parameter) and detecting anomalies by using an invariant detection algorithm based on a separation of the measurement space into disjoint subspaces corresponding to normal and anomalous network traffic. The method's ability to detect anomalies is evaluated in real traffic from Abilene, a United States backbone network. The theoretically expected results are confirmed.

INTRODUCTION

Network management becomes very complex as networks increase in size and complexity. The traffic demand in 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 a denial-of-service attack, a flash crowd event, a virus/worm propagation, etc.) that spans multiple physical links of the network. The reliable detection of these unusual and significant changes in the OD traffic matrix is an important issue for network operation.

High hardware requirements are necessary to networkwide collect and process the direct OD flow measurements [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. For this reason detecting an anomaly in the traffic matrix from SNMP measurements is a difficult task.

Several approaches are proposed to remedy this problem. The first group of methods consists in detecting anomalies in SNMP measurements without taking into account the traffic matrix. Such methods typically use time series (AutoRegressive Integrated Moving Average or ARIMA models among others) [START_REF] Thottan | Anomaly detection in IP networks[END_REF][START_REF] Tartakovsky | A novel approach to detection of intrusions in computer networks via adaptive sequential and batch-sequential change-point detection methods[END_REF][START_REF] Krishnamurthy | Sketch-based change detection: methods, evaluation, and applications[END_REF][START_REF] Ahmed | Multivariate online anomaly detection using kernel recursive least squares[END_REF] to model the SNMP measurements' evolution in time and 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 linear mathematical relation between the OD traffic matrix and the SNMP measurements, represented by the routing matrix. Moreover, they cannot be applied when the routing matrix varies in time (dynamic routing) and they cannot be used to estimate the OD traffic matrix. The second group of methods [START_REF] Coates | Internet tomography[END_REF][START_REF] Tebaldi | Bayesian inference on network traffic using link count data[END_REF][START_REF] Cao | Time-varying network tomography: router link data[END_REF] exploits this linear mathematical relation. 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 different models for SNMP measurements (wavelets, PCA among others) and OD flows (ARIMA time series) to detect anomalies. A major drawback of these methods is the lack of theoretical results on their optimality properties (maximization of the probability to detect an anomaly for example) of the studied methods. Finally, the last group of methods consists in using the Kalman filtering technique [START_REF] Soule | Traffic matrices: balancing measurements, inference and modelling[END_REF] to model the time evolution of the traffic matrix and to detect changes in the OD flows. Unfortunately, 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 main contributions of this study are the following: firstly, 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 highly non-stationary in time measurements and to dynamic routing. Finally, an optimal invariant detection algorithm is proposed to detect anomalies directly from SNMP measurements (no need of direct OD flow measurements). This algorithm is optimal in the sense that it maximises the probability of detecting the anomalies under a constrained false alarm probability.

PROBLEM STATEMENT

This section briefly presents the SNMP measurement model and the anomaly detection problem.

SNMP Measurements

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]. The volume of traffic y(ℓ), 2) y [START_REF] Tartakovsky | A novel approach to detection of intrusions in computer networks via adaptive sequential and batch-sequential change-point detection methods[END_REF] x [START_REF] Coates | Internet tomography[END_REF][START_REF] Thottan | Anomaly detection in IP networks[END_REF] x(1, 2)

x(1, 3)

x(1, 3)

(3, Network t 0 t 0 t 0 t t t y(1) = x(1, 2) + x(1, 3) + x(3, 2)
Figure 1: Detection of unusual changes in OD traffic volumes. A step-wise anomaly appears at time t 0 in the OD flow x(1, 3) and it is routed on links y(1) and y [START_REF] Thottan | Anomaly detection in IP networks[END_REF].

typically in bytes, on the link ℓ at time t (to simplify the notations, the subscript t is omitted) is provided by SNMP link load measurements [START_REF] Zhang | Network anomography[END_REF]. Let x(i, j) be the OD traffic demand from node i to node j at time t. This situation is shown in Here, X T denotes the transpose of the matrix X. Without loss of generality, the known matrix A is assumed to be full row rank, rank (A) = n.

Clearly, since m ≫ n, it is impossible to infer X from Y .

Volume anomaly detection problem

The detection problem consists in detecting 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) and y(3) (typical anomalies to be detected in SNMP measurements are presented in [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF]). As it was mentioned above, the main problem with gathering the traffic matrix from SNMP measurements is that n ≪ m.

To overcome this problem a parsimonious linear model of non-anomalous traffic has been used. The idea of this model is that the non-anomalous (ambient) traffic X can be represented at each time t by using a known family of q basis functions B = (b 1 , b 2 , . . . , b q ) such that q < n. Therefore, the ambient traffic can be expressed as X ≈ Bµ where the m × q matrix B is assumed to be known and µ ∈ R q is a vector of unknown coefficients which describes the OD flows' variation with respect to the set of vectors b i . The advantages of such a parametric model are the following: firstly, a nonparametric 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 [START_REF] Ringberg | Sensitivity of PCA for traffic anomaly detection[END_REF]. Secondly, the parametric detection methods' performance is better than the non-parametric one provided the adopted model is accurate enough. Therefore, the parametric approach is used to design a reliable detection algorithm derived from a spline-based OD flow model.

OD FLOW MODEL

The derivation of the ambient traffic matrix model consists of two different steps: i) description of the ambient traffic by using a spatial stationary gravity model and ii) linear approximation of the gravity model by using polynomial splines.

General gravity model

Gravity models [START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF][START_REF] Gunnar | Traffic matrix estimation on a large IP backbone -A comparison on real data[END_REF] assume that the traffic volume between two nodes is directly proportional to the relative attraction of each node and inversely proportional to some function of the separation (typically the distance) between these nodes. Let us assume that, at time t, the node i of the network is characterized by an attractive factor A(i) related to the incoming flow at node i and a repulsion factor R(i) associated to the outgoing flow from node i. A network-adapted formulation of the gravity model [START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF] may be given by the following equation:

x(i, j) ≈ h (R(i) • A( j)) • d(i, j) -1 , (1)
where d(i, j) is the deterrence function representing the separation between nodes i and j at time t and h : R → R is a weighting function relating the repulsion and attraction factors. Typically, the function h is a non-decreasing smooth function, which is a reasonable assumption in the absence of another physical model on the network. The function h is unknown and makes it possible to model more complex interactions between nodes and to give more flexibility to the model as suggested (but not used) in [START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF][START_REF] Gunnar | Traffic matrix estimation on a large IP backbone -A comparison on real data[END_REF]. Following [START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF], it is assumed that d(i, j) = d for all 1 ≤ i, j ≤ r where d is an unknown real value. This assumption is quite reasonable in the case of an autonomous system which is geographically limited, as the Abilene network [START_REF] Zhang | Fast accurate computation of large-scale IP traffic matrices from link loads[END_REF] for example, since the distance between nodes are not significant for end-users. Hence, it follows that (1) can be rewritten as:

x(k) = x(i k , j k ) ≈ d -1 • h(ω(k)) (2)
where ω(k) = R(i k ) • A( j k ) represents the importance of OD flow k in terms of traffic volume outgoing from i k and incoming to j k . Without any loss of generality, it is assumed that 0 ≤ ω(k) ≤ 1 for all k.

Spatial stationarity

Firstly, it is assumed that the order of importance between OD flows remains stable in time, i.e. ω(k

1 ) ≤ ω(k 2 ) ≤ . . . ≤ ω(k m ), (3)
for all time t, where k j is the index of the OD flow with the j-th rank of importance (the relevance of this assumption is shown in section 5). Secondly, one of the few invariants of Internet traffic is "the elephants and mice phenomenon" [START_REF] Medina | Traffic matrix estimation : existing techniques and new directions[END_REF]: a few percentage of flows contributes to a large proportion of total traffic. Hence, OD flows can typically be classified in three classes of traffic depending of their volume of traffic: dominating OD flows, negligible OD flows and medium-size OD flows. This last class of traffic is defined to add more flexibility to the classification of OD flows. These classes induce a segmentation of the interval [0; 1] in three parts [0; π 1 ], [π 1 ; π 2 ] and [π 2 ; 1], with 0 < π 1 < π 2 < 1, corresponding respectively to negligible, medium and dominating OD flows.

As it has been experimentally confirmed, the segmentation is time-stable, at least during several days (see section 5).

OD flow spline-based model

The function h is assumed to be non-decreasing with a certain smoothness and defined piecewise, with respect to knots π 1 and π 2 , on the interval [0; 1]. Therefore, it is linearly approximated by using polynomial splines (basic definitions and results on polynomial splines can be found in [START_REF] Nürnberger | Approximation by spline functions[END_REF]). Let {ω → b 1 (ω), . . . , ω → b q (ω)} be a basis of q = p + 3 functions for the space of splines of degree p with p -1 continuous derivatives and 2 knots. By using (2), there exists a unique time-dependent vector µ = (µ(1), . . . , µ(q)) T such as X ≈ Bµ where B is the m × q matrix whose element at position (i, j) is b j (ω(i)). Here, since the function h is unknown, sampling points ω(i) can be arbitrarily chosen in the interval [0; 1] provided condition (3) is verified. The existence of dominant OD flows is sometimes modelled by using α-stable laws [START_REF] Doukhan | Theory and Applications of Long-Range Dependence[END_REF]. In such a context, the spline-based approximation is naturally justified by the necessity to have a piecewise approximation of such a heavy-tailed power law distribution. Finally, it is assumed that model errors together with the natural variability of the OD flows follow a spatial Gaussian distribution [START_REF] Coates | Internet tomography[END_REF], which leads to the model:

X = Bµ + ξ (4)
where ξ ∼ N (0, γ 2 Σ) is a Gaussian noise with the m × m spatial diagonal covariance matrix Σ = diag(σ 2 1 , . . . , σ 2 m ). The matrix Σ is assumed to be known and stable in time. On the contrary, the scalar γ 2 serves to model the mean level of the variance (due to the natural OD flow time variability) and it may depend on the time. In practice, Σ and γ 2 are estimated from a few anomaly-free SNMP measurements.

ANOMALY DETECTION PROBLEM

The goal is to detect an anomaly with the highest probability of detection for a given probability of false alarm, i.e. activating an alarm in absence of anomaly, which is an undesirable event.

Hypotheses testing: problem statement

According to the previous section, the non-anomalous link load measurement model is given by the following linear model :

Y = A Bµ + Aξ = Gµ + ζ , (5)
where Y = (y(1), . . . , y(n)) T and ζ ∼ N (0, γ 2 AΣA T ). Without any loss of generality, the resulting matrix G = A B is assumed to be full column rank. Since the matrix Φ = AΣA T is known, the testing problem consists of choosing between the two alternatives:

H 0 = {Z ∼ N (θ + Hµ, γ 2 I n ); θ = 0, µ ∈ R q } (6) H 1 = {Z ∼ N (θ + Hµ, γ 2 I n ); θ = 0, µ ∈ R q }, (7)
with

Z = Φ -1 2 Y , H = Φ -1 2 G, Φ -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 µ is considered as a nuisance vector parameter since i) it is completely unknown, ii) it is not necessary for the anomaly detection and iii) it can mask the anomalies. Typically, when an anomaly occurs in OD flow j, the vector θ has the form θ = ε Φ -1 2 a( j) where a( j) is the j-th normalized column of A and ε is the intensity of the anomaly.

The aim is to detect the presence of an anomalous vector θ not explicable by the ambient traffic model (4).

Let K α = {φ : sup µ∈R q Pr θ =0,µ (φ (Z) = H 1 ) ≤ α} be the class of tests φ : R n → {H 0 , H 1 } with upper-bounded maximum false alarm probability, where the probability Pr θ ,µ stands for the vector of observations Z being generated by the distribution N (θ + Hµ, γ 2 I n ) and α is the prescribed probability of false alarm. The power function β is defined with the probability of correct detection: β (θ ; µ) = Pr θ =0,µ (φ (Z) = 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 µ 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 } (see an introduction to the principle of invariance in [START_REF] Lehman | Testing Statistical Hypotheses[END_REF]). The maximal invariant statistics (also called "parity vector") U = W Z is the transformation of the measured output Z into a set of nq linearly independent variables by projection onto the left null space of the matrix H. The matrix W T = (w 1 , . . . , w n-q ) of size n × (nq) is composed of the eigenvectors w 1 , . . . , w n-q of the projection matrix P ⊥ H = I n -H(H T H) -1 H T corresponding to eigenvalue 1. The matrix W satisfies the following conditions: W H = 0, W T W = P ⊥ H , WW T = I n-q . Let S be the family of surfaces

S = {S c : c > 0} with S c = {θ : P ⊥ H θ 2 2 /γ 2 = c 2 }. Then, it is shown [19] that the test φ * (Z) = H 0 if Λ(Z) = P ⊥ H Z 2 2 /γ 2 < λ α H 1 else , (8)
where the threshold λ α is chosen to satisfy the false alarm bound α, Pr θ =0,µ (Λ(Z) ≥ λ α ) = α, 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 θ /γ 2 .

NUMERICAL RESULTS

This section shows the relevance of the model and the performance of the detection algorithm.

Description of the data set

The evaluation of the proposed methods requires the knowledge of the real OD traffic flows. Such measurements are quite difficult to obtain in a commercial network but are available for the Abilene network. The Abilene backbone is composed of r = 12 core routers and m = 144 OD flows. For these numerical experiments, n = 42 backbone links are measured. More details on this network are given in [START_REF] Zhang | Fast accurate computation of large-scale IP traffic matrices from link loads[END_REF] and real data are available in [START_REF] Zhang | 6 months of Abilene traffic matrices[END_REF]. The primary data inputs are the time series of link loads (bytes across interfaces) gathered through SNMP. The sampling rate is one measurement per 10 minutes, i.e. each measurement corresponds to the total of volume of traffic (in bytes) which has passed through a given link during 10 minutes. Two sets of measurements are used: the first one, the anomaly-free data set, is composed of 6 anomaly-free SNMP measurements (one hour measurement period) and the second one, the testing data set, is composed of 720 SNMP measurements (five days measurement period). Let T a (respectively T b ) be the set of time index associated to SNMP measurements of the anomaly-free (resp. testing) data set. The anomaly-free data set is measured one hour before the testing one.

To identify the set of "true" anomalies in the testing data set (as a precursor to the validation step), unusual deviations from the mean in each OD flow are manually detected. Manual inspection declares an anomaly if the unusual deviation intensity of the guilty OD flow leads to an augmentation of traffic 1) larger than 1.5% of the total amount of traffic on the network and 2) larger than 1% of the amount of the traffic carried by the links routing this guilty OD flow. Hence, only significant volume anomalies are considered as "true anomalies" (small volume anomalies have little influence on link utilisations). Let T • b ⊂ T b be the set of time index t associated to the 680 non consecutive SNMP measurements of the testing data set manually declared as anomaly-free (40 measurements of the testing data set are affected by at least one significant volume anomaly).

Numerical validation of the model

Although many aspects could potentially be included in the evaluation, the focus is on the potential impact of performance errors on traffic engineering tasks. Hence, the root mean square error (RMSE) is used:

RMSE label (t) = m ∑ i=1 x label t (i) -x t (i) 2 , ∀t ∈ T • b .
Here, x t (i) denotes the true traffic volume of OD flow i at time t ∈ T • b and x label t (i) denotes the corresponding estimate for the method entitled 'label'. Three estimates are compared : 1) simple gravity estimate [START_REF] Zhang | Fast accurate computation of large-scale IP traffic matrices from link loads[END_REF] with label 'SG', 2) tomogravity estimate [START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF] with label 'TG' and 3) spline-based Maximum Likelihood (ML) estimate with the label 'SML'. Since the measurement model ( 5) is a Gaussian linear one, the optimal estimate of X t is the ML estimate [START_REF] Rao | Linear statistical inference and its applications[END_REF] X SML t given by:

X SML t = B(H T H) -1 H T Z t .
The statistical properties of the ML estimate are well known [START_REF] Rao | Linear statistical inference and its applications[END_REF] contrary to the simple gravity and tomogravity estimates [START_REF] Zhang | Fast accurate computation of large-scale IP traffic matrices from link loads[END_REF][START_REF] Zhang | Estimating point-to-point and point-to-multipoint traffic matrices: an information-theoretic approach[END_REF]. The temporal correlation of the noise sequence (ζ t ) t≥1 is ignored for the following reason: it can be theoretically shown that the integration of the ARMA (Au-toRegressive and Moving Average) model of (ζ t ) t≥1 does not change the covariance matrix (and the bias) of X SML t . The spline-based model is computed by using SNMP measurements of the short anomaly-free data set : 1) the tomogravity estimate x T G t (k) is computed for all OD flow k and all t ∈ T a , 2) the mean flow values x T G (k) = ∑ t∈T a

x T G t (k) are computed and 3) sorted in ascending order to obtain a rough estimate of the OD flow ranks. The splinebased model is designed with cubic splines (p = 3) and knots π 1 = 0.8507 and π 2 = 0.9830 with sampling points ω(k) uniformly distributed in the interval [0; 1]. Small variations on the values π 1 and π 2 have no serious effect on the results.

Method

SG TG SML Total RMSE 9 336.9 3 934.9 3 765.6 Table 1: Total RMSE (in kilobytes) for 680 anomalyfree measurements for gravity (SG), tomogravity (TG) and spline-based (SML) models.

The mean value x T G (k) is also used as an estimate σ 2 k of σ 2 k , which leads to an estimate Φ of Φ (quite efficient and sufficient in practice). In this estimation step, it is not necessary to know γ 2 t as explained in [START_REF] Rao | Linear statistical inference and its applications[END_REF]. b RMSE(t) 2 on this time period is computed in Table 1 as a global indicator to compare the methods. Clearly, the spline-based estimate gives better results than the others. To verify the spatial Gaussian assumption, residuals U t are computed for each t ∈ T • b . The Kolmogorov-Smirnov test [START_REF] Lehman | Testing Statistical Hypotheses[END_REF] at the level 5% accepts the Gaussian hypothesis for 670 of these measurements (acceptation 98.5% of the time). also plotted on the same figure with dashed lines. It shows that the "shape" of the OD flows is well modeled by the proposed spline-based model and is stable over the time. Hence, the spatial stationarity assumption is quite reasonable.

Numerical validation of the detection algorithm

The detection algorithm is applied to the SNMP measurements of the testing data set. For the detection purpose, it is crucially important to have a good estimate of the noise level γ 2 t . This parameter is estimated from the short anomaly-free data set by using the ML estimate of noise variance [START_REF] Rao | Linear statistical inference and its applications[END_REF] The results are presented in Table 2. The second column shows that the proposed test obtains a false alarm rate of 1.18% comparable with the prescribed level of false alarm α = 0.01. The probability to detect a volume anomaly is about 77.5%. The third column presents the results obtained by the PCA test described in [START_REF] Lakhina | Diagnosing network-wide traffic anomalies[END_REF]. The threshold of this test is chosen to obtain a similar false alarm rate of 1.03%. Clearly, the PCA test is not as sensitive (correct detection rate about 12.50%) as the proposed test. Indeed, the PCA decomposition of SNMP measurements is too rough to detect small (but significant) anomalies. Finally, Fig. 4 shows the correct detection rates of the spline-based test and the PCA test for different false alarm rates varying between 0 and 1. For example, for a correct detection rate of 80%, the false alarm rate of the spline-based test (≈ 2%) is clearly better than that of the PCA test (≈ 70%). Hence it can be concluded that the numerical results confirm the theoretical properties of the detection method and shows that the proposed test outperforms the PCA approach, at least for the Abilene data set.

CONCLUSION

The problem of anomaly detection in OD traffic volume from SNMP measurements has been considered as a statistical hypotheses testing problem with nuisance parameters (nonanomalous traffic). Since the number of SNMP measurements is significantly lower than the number of OD flows, an original linear spline-based parsimonious model is proposed to describe the non-anomalous traffic and to overcome the ill-posed nature of the SNMP measurement model. Results obtained with real data traffic from a United States backbone network show that both the OD traffic matrix estimation and the volume anomaly detection approaches outperform the methods previously applied in the field.

Fig. 1 .

 1 The link loads and the traffic matrix are simply related by a linear equation Y = A X where Y = (y(1), . . . , y(n)) T , X = (x(1), . . . , x(m)) T contains the m (m ≫ n) unknown traffic matrix elements x(k) = x(i k , j k ) written as a vector and A = {a(ℓ, k)} is the n × m routing matrix where 0 ≤ a(ℓ, k) ≤ 1 represents the fraction of OD flow k volume routed through link ℓ.

Figure 2 :

 2 Figure 2: Comparison between the SG, TG and SML RMSE for 680 anomaly-free measurements.

Fig. 2

 2 Fig.2shows the error RMSE(t) over the set T • b . The xaxis of the figure has no time meaning and it corresponds to the index of each measurement. The sum ∑ t∈T • b RMSE(t) 2 on this time period is computed in Table1as a global indicator to compare the methods. Clearly, the spline-based estimate gives better results than the others. To verify the spatial Gaussian assumption, residuals U t are computed for each t ∈ T • b . The Kolmogorov-Smirnov test[START_REF] Lehman | Testing Statistical Hypotheses[END_REF] at the level 5% accepts the Gaussian hypothesis for 670 of these measurements (acceptation 98.5% of the time).

Figure 3 :

 3 Figure 3: Approximation of real OD flows (full lines) by the spline-based model (dashed lines). Typical non-anomalous OD flows, sorted in the increasing order of their volume of traffic, are shown as functions of time t in Fig. 3. Since the "shape" of sorted OD flows is almost constant over the time, only a few sorted OD flows are plotted. The SML estimated profiles of the OD flows are

Figure 4 :

 4 Figure 4: Correct detection rate versus false alarm rate for the spline-based test (solid line) and the PCA test (dotted line).

Table 2 :

 2 Results of the detection for 720 measurements composed of 680 anomaly-safe measurements and 40 anomalous measurements for the spline-based and PCA tests.

	Type of situation	Spline-based	PCA
	Normal working	672 (98.82 %) 673 (98.97 %)
	False alarms	8 (1.18 %)	7 (1.03 %)
	Missed detections	9 (22.50 %)	35 (87.50 %)
	Correct detections	31 (77.50 %)	5 (12.50 %)

  in residuals U t . Since this parameter can slowly vary in time, its value is updated during the test. During the test, at time t, if no anomaly has been declared one hour before, γ 2 t is estimated by its value one hour before.

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

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