Method and apparatus for network-level anomaly inference

1. Field of the Invention

Embodiments of the present invention generally relate to network monitoring systems and, more particularly, to a method and apparatus for network-level anomaly inference.

2. Description of the Related Art

Today's networks are large, geographically distributed and constantly evolving. Currently anomaly detection is limited to anomaly detection at individual sites or links. However, network-wide views of traffic matrices are critical inputs for all manual or automated processes for large networks, including demand forecasting, capacity planning, traffic and network engineering, and operational maintenance. Going beyond anomaly detection at the link level to explanation at the network level via widely available network data would provide the necessary network-wide views. Therefore, there exists a need in the art for network-level anomaly inference.

Method and apparatus for network-level anomaly inference in a network is described. In one embodiment, link load measurements are obtained for multiple time intervals. Routing data for the network is obtained. Link level anomalies are extracted using temporal analysis on the link load measurements over the multiple time intervals. Network-level anomalies are inferred from the link-level anomalies.

So that the manner in which the above recited features of the present invention can be understood in detail, a more particular description of the invention, briefly summarized above, may be had by reference to embodiments, some of which are illustrated in the appended drawings. It is to be noted, however, that the appended drawings illustrate only typical embodiments of this invention and are therefore not to be considered limiting of its scope, for the invention may admit to other equally effective embodiments.

FIG. 1 is a block diagram depicting an exemplary embodiment of a network architecture in accordance with one or more aspects of the invention;

FIG. 2 is a flow diagram depicting an exemplary embodiment of a method for identifying network-level anomalies in accordance with one or more aspects of the invention; and

FIG. 3 is a block diagram depicting an exemplary embodiment of a computer suitable for implementing the processes and methods described herein.

FIG. 1 is a block diagram depicting an exemplary embodiment of a network architecture 100 in accordance with one or more aspects of the invention. The network architecture 100 illustratively comprises a network 102 , routers 104 - 1 through 104 - 5 (collectively referred to as routers 104 ), links 106 - 1 through 106 - 6 (collectively referred to as links 106 ), a database server 108 , a route monitor 112 , and an application server 110 . The network 102 comprises a packet network configured to propagate packets in accordance with a particular network protocol, such as internet protocol (IP), and transport protocol, such as transmission control protocol (TCP), user datagram protocol (UDP), and the like.

The routers 104 are interconnected by the links 106 . Each of the routers is configured to provide link load (i.e., traffic volume) measurements at predefined time intervals for storage in the database server 108 . For example, the routers 104 may employ the simple network management protocol (SNMP) to provide the link load measurements to a network management server 109 , which stores link measurement data 114 in the database server 108 . Operation of SNMP to collect such data is well known in the art.

The routers 104 employ routing protocols to facilitate the exchange of routing information. The routing protocols are used to build forwarding tables in the routers 104 . Exemplary routing protocols include interior gateway protocols (IGPs), such as open shortest path first (OSPF) or intermediate system to intermediate system (ISIS), and exterior gateway protocols (EGPs), such as the border gateway protocol (BGP). Operation of such routing protocols to exchange routing information is well known in the art. The route monitor 112 is configured to provide up-to-date measurements of network routing employed by the routers 104 (“routing data”).

The application server 110 is configured to implement a network-level anomaly analyzer 116 . The network-level anomaly analyzer 116 is configured to obtain link load measurements from the database server 108 and routing data from the route monitor 112 and to identify anomalous origin-destination (OD) flows network-wide. Operation of the network-level anomaly analyzer 116 is described below. Those skilled in the art will appreciate that the network architecture 100 is merely illustrative and that the network-level anomaly analyzer 116 may be employed in a myriad of network architectures. Notably, a typical large (e.g., Tier 1) network may include hundreds of routers and thousands of links. In general, the network-level anomaly analyzer 116 is configured to process link load measurements and routing data to identify anomalous OD flows network-wide.

FIG. 2 is a flow diagram depicting an exemplary embodiment of a method 200 for identifying network-level anomalies in accordance with one or more aspects of the invention. The method 200 may be employed by the network-level anomaly analyzer 116 of FIG. 1 . The method 200 begins at step 201 . At step 202 , link load measurements are obtained for multiple time intervals. At step 204 , routing data for the network is obtained. At step 206 , link level anomalies are extracted using temporal analysis of link load measurements over time. At step 208 , network level anomalies are inferred from the link level anomalies. The method 200 ends at step 299 . Aspects of the method 200 are described in more detail immediately below.

In one embodiment of the invention, anomalies are extracted at the link level and explanations for the anomalies are inferred at the network-level, using a tomograhic approach. This process is referred to herein as network anomography. Network tomography describes several problems: inferring topology, or link performance of a network from end-to-end measurements, or inferring Origin-Destination (OD) traffic demands from link traffic measurements. These problems can be written using the following linear relationship:

b=Ax   Eq. 1,

where b contains link load measurements, A is the routing matrix, and x includes the unknown traffic matrix elements written as a vector. For a network with n links, and m OD flows, we define the routing matrix to be the n×m matrix A=[a ij ] where a ij indicates the fraction of traffic from flow j to appear on link i.

The problem of inferring the OD traffic-matrix is massively under-constrained: there will be approximately N 2 OD flows to estimate with only O(N) link measurements. Hence tomography methods seek to introduce additional information, often in the form of some kind of traffic model. Anomography problems are also highly under-constrained, but the models used to describe traffic are inappropriate for anomalies—by definition these events are generated by completely different processes from normal network traffic. Furthermore, in anomography, detection is combined with inference. In contrast, in standard network tomography problems, only a set of traffic matrix elements is inferred. Hence there are important differences between anomography and network tomography.

It is also important to note that routing matrices change over time. While routing matrices may be taken to be constant, in some cases it is important to allow for the fact that routing is not constant, and neither is the traffic. In order to allow for variable routing, the traffic measurements and the routing matrix are indexed over time. Thus, the relationship between the link traffic and OD traffic matrix may be represented ast:

b j =A j x j   Eq. 2,

where A j is an n×m routing matrix, x j is a length-n vector of unknown OD flow traffic volumes, and b j is a length-m vector of link loads, at time interval j.

In one embodiment, it is assumed that the routing matrices A j are time-invariant and are denoted by A. Under this assumption, we can combine all t linear systems (Eq. 2) can be combined into a single equation using matrix notation:

B=AX   Eq. 3,

where B=[b 1 b 2 . . . b t ] is the matrix formed by having b j as its column vectors, and X=[x 1 x 2 . . . x t ] is the matrix with x j as its column vectors.

In the present embodiment, a late-inverse solution strategy to network anomography is employed. Whereas early-inverse solutions infer OD flow data from link load measurements and anomaly detection is applied, late-inverse approaches extract the anomalous traffic from the link load observation, then form and solve a new set of inference problems:

{tilde over (B)}=A{tilde over (X)}   Eq. 4,

where {tilde over (B)}=[{tilde over (b)} 1 {tilde over (b)} 2 . . . {tilde over (b)} t ] is the matrix of anomalous traffic in the observables, and {tilde over (X)}=[{tilde over (x)} 1 {tilde over (x)} 2 . . . {tilde over (x)} t ] is the matrix of OD flow anomalies to be diagnosed, over t time intervals. The idea of the late-inverse method is to defer “lossy” inference to the last step.

While the new inference problems (Eq. 4) share the same linear-inverse structure as in network tomography (Eq. 3), the characteristics of the unknowns are different, and so is the solution strategy. A framework for late-inverse anomography methods is described. In this framework, {tilde over (B)} is formed by multiplying B with a transformation matrix T. Depending on whether a left or right multiplying transformation matrix is used, the framework can be furthered divided into the following two classes: (1) spatial anomography, where a left multiplying transformation matrix T is used to form {tilde over (B)}, i.e., {tilde over (B)}=TB; and (2) temporal anomography, where a right multiplying transformation matrix T is used to form {tilde over (B)}, i.e., {tilde over (B)}=BT. It is worth noting a subtle distinction between spatial and temporal anomography. Temporal anomography has the nice property that it is self-consistent in that given B=AX and {tilde over (B)}=BT, a solution to Eq. 4 can always be obtained by setting {tilde over (X)}=XT. In other words, we can apply exactly the same analysis can be applied on B and on X to obtain the corresponding anomalies {tilde over (B)} and {tilde over (X)}, respectively (i.e., the same transformation matrix T is used). Spatial anomography, in contrast, does not have this property, because with B=AX and {tilde over (B)}=TB, the formulated problem instance (Eq. 4) is not guaranteed to have a solution.

Four embodiments of temporal anomography analysis techniques for extracting anomalous traffic {tilde over (B)} from link load observations B are disclosed, which include ARIMA modeling, Fourier analysis, Wavelet analysis, and temporal Principle Components Analysis (PCA). In general, the anomalous link traffic can be separated by performing temporal analysis on the time series of each link. Consider a set of link traffic data over time t: B=[b 1 b 2 . . . b t ]. The process of extracting anomalies by exploiting the temporal structure within the data points can be modeled as a linear transformation of the time series: {tilde over (B)}=[{tilde over (b)} 1 {tilde over (b)} 2 . . . {tilde over (b)} t ]=BT, where the transformation matrix T can be either explicit or implicit.

In one embodiment, autoregressive integrated moving average (ARIMA) modeling is used. The ARIMA methodology is a class of linear time-series forecasting techniques that capture the linear dependency of the future values on the past. It is able to model a wide spectrum of time-series behavior. An ARIMA model includes three order parameters: the autoregressive parameter (p), the number of differencing passes (d), and the moving average parameter (q). Models can be summarized as ARIMA(p,d,q). A model described as ARIMA(0,1,2) means that it contains p=0 (zero) autoregressive parameters and q=2 moving-average parameters, which were computed for the time series after it was differenced once (d=1).

A general ARIMA model of order (p, d, q) can be expressed as:

z k - ∑ i = 1 p ⁢ ⁢ ϕ i · z k - i = e k - ∑ j = 1 q ⁢ ⁢ θ j · e k - i , Eq . ⁢ 5

where z k is obtained by differencing the original time series d times (when d≧1) or by subtracting the mean from the original time series (when d=0), e k is the forecast error at time k, φ i (i=1, . . . , p) and θ j (j=1, . . . , q) are the autoregression and moving average coefficients, respectively.

There are well known techniques for estimating the parameters p, d, q, φ i and θ j for a given time series, and given the parameters. The model is applied by predicting the terms z k and then estimating the prediction errors by e k+1 =z k+1 −{circumflex over (z)} k+1 , which then form the anomalous traffic (the traffic which does not fit the model). In practice, the parameters used in the ARIMA model are sometimes chosen to meet particular goals intended by the implementer, rather than being estimated from a particular data set. This is because the parameters of use in a particular data set may change over time. However, the present invention employs adaptive techniques to overcome this problem, as described below.

If the time series are considered to be vectors of length t, then the above results can be written in matrix form. Taking the measurements b=(b 1 , . . . , b t ) T , the errors e=(e 1 , . . . , e t ) T can be obtained via right-multiplication by a transformation matrix {tilde over (b)} T =e T =b T T. Specifically, let I denote the t×t identity matrix, ∇ denote the “back shift” matrix, and 1 denote the t×t unit matrix, i.e.,

I = [ 1 0 0 … 0 0 0 1 0 … 0 0 … ⁢ 0 0 0 … 1 0 0 0 0 … 0 1 ] , V _ = [ 0 1 0 … 0 0 0 0 1 … 0 0 … 0 0 0 … 0 1 0 0 0 … 0 0 ] , 1 = [ 1 1 1 … 1 1 1 1 1 … 1 1 … 1 1 1 … 1 1 1 1 1 … 1 1 ] .

The differencing result, z=[z 1 z 2 . . . z t ] T , can then be represented by:

z T = { b T ⁡ ( I - ∇ ) d , for ⁢ ⁢ d ≥ 1 , b T - 1 n ⁢ b T ⁢ 1 = b T ⁡ ( I - 1 n ⁢ 1 ) , for ⁢ ⁢ d = 0 Eq . ⁢ 6

Equation 5 can be written in matrix notation as

z T - ∑ i = 1 p ⁢ ⁢ ϕ i ⁢ z T ⁢ ∇ ⅈ = ⅇ T - ∑ j = 1 q ⁢ ⁢ θ j ⁢ ⅇ T ⁢ ∇ j ,

or equivalently as

Extending ARIMA based models to multivariate time series is straightforward. As noted earlier, the matrix B is constructed with the measurements at each time period b i as its columns. Using the transformations just described, the following is obtained:

Replacing Z by the matrix form of Equation 6, it is seen that E=BT is indeed a transformation given by right-multiplying B with a matrix T. In fact, any linear filtration of the elements of a time series can be modeled by a right multiplying matrix transformation. If the transformation is time-invariant, then the matrix in question will be Toeplitz (the values along diagonals will be constant). For anomaly detection, the forecast errors are identified as anomalous link traffic, {tilde over (B)}=E. That is, traffic behavior that cannot be well captured by the model is considered anomalous.

In another embodiment, Fourier analysis is used. Fourier analysis is the process of decomposing a complex periodic waveform into a set of sinusoids with different amplitudes, frequencies, and phases. The sum of these sinusoids can exactly match the original waveform. This lossless transform presents a new perspective of the signal under study (in the frequency domain), which has proved useful in very many applications.

For a discrete-time signal x 0 , x 1 , . . . , x N−1 , the Discrete Fourier Transform (DFT) is defined by:

f n = 1 N ⁢ ∑ k = 0 N - 1 ⁢ ⁢ x k ⁢ ⅇ - j ⁢ ⁢ k2π ⁢ ⁢ n / N , for ⁢ ⁢ 0 ≤ n ≤ N - 1 ,

where f n is a complex number that captures the amplitude and phase of the signal at the n-th harmonic frequency (with base frequency 1/N) and j denotes the √{square root over (−1)}. Note that for a real signal {f n } is symmetric, i.e., f n =f N−1−n . Lower n corresponds to a lower frequency component, with f 0 being the DC component, or the average of the input series, and f n with n close to N/2 corresponding to high frequencies.

The Inverse Discrete Fourier Transform (IDFT) is used to reconstruct the signal in the time domain by

x n = ∑ k = 0 N - 1 ⁢ ⁢ f k ⁢ ⅇ j ⁢ ⁢ k ⁢ ⁢ 2 ⁢ π ⁢ ⁢ n / N ⁢ ⁢ for ⁢ ⁢ 0 ≤ n ≤ N - 1.

An efficient way to implement the DFT and IDFT is through an algorithm known as the Fast Fourier Transform (FFT). The computational complexity of the FFT is O(N log(N)).

The FFT can be used to extract anomalous link traffic, {tilde over (B)}, by filtering out the low frequency components in the link traffic time series. In general, low frequency components capture the daily and weekly traffic patterns, while high frequency components represent the sudden changes in traffic behavior. Working in the frequency domain provides the opportunity to distinguish these two kinds of behaviors. FFT based anomography is as follows: First, transform link traffic B into the frequency domain: F=FFT(B): apply FFT on each row of B (recall that a row corresponds to the time series of traffic data on one link). The result is the corresponding frequency domain series, in each row of F. Second, remove low frequency components: i.e., set F i =0, for i∈[1,c]∪[N−c,N], where F i is the i-th column of F and c is a cut-off frequency (for example, 10-minute aggregated link traffic data of one week duration can be used with c=171, corresponding to a frequency of one cycle per hour). Third, transform back into the time domain: i.e., take {tilde over (B)}=IFFT(F). The result is the high frequency components in the traffic data, which is used as anomalous link traffic, {tilde over (B)}.

The DFT and IDFT may be represented as right-matrix products. In setting columns of F to zero, and performing the IDFT we are taking a linear combination of the columns of F, which in turn are a linear combination of those of X. Hence, the overall process above can be modeled as a right-multiplying matrix transformation {tilde over (B)}=BT. Note also that in thresholding at frequency c, the symmetry of F is preserved, and so although F may contain complex elements, the resulting transform will be real.

In another embodiment, wavelets are used to extract anomalous link traffic. Wavelets are mathematical functions that segregate data into different frequency components, and then study each component with a resolution matched to its scale. They provide a powerful means for isolating characteristics of signals via a combined time-frequency representation and are often considered superior to traditional Fourier methods, especially in situations where the signal contains transients, such as discontinuities and sharp spikes.

As described below, {tilde over (B)} can be computed as the high-frequency components of link traffic B. {tilde over (B)} can also be computed as the mid-frequency components of B in essentially the same way. First, use wavelets to decompose B into different frequency levels: W=WAVEDEC(B), by applying a multi-level 1-D wavelet decomposition on each row of B. The result is a wavelet decomposition vector, which we save as one row in matrix W. An exemplary wavelet that can be used is the Daubechies wavelet of order 6. Second, remove low- and mid-frequency components in W by setting all coefficients at frequency levels higher than w c to 0. Here, w c is a cut-off frequency level. For example, 10-minute aggregated link traffic data of one week duration can be used, and w c can be set at 3. That is, only keep coefficients at frequency levels 1, 2, and 3. Third, reconstruct the signal: {tilde over (B)}=WAVEREC(B). The result is the high-frequency components in the traffic data.

It can be verified that the process of WAVEDEC and WAVEREC only involves linear combinations of columns of B. As a result. the {tilde over (B)} derived through the wavelet based anomography can also be modeled as right multiplying matrix transformation.

In another embodiment, temporal PCA is employed to extract the link-level anomalies. Consider the link load matrix B=[b 1 b 2 . . . b t ]. We can think of each row as a t-dimensional vector. What we are looking for is a new coordinate system, v 1 , v 2 , . . . , v t , such that the projection of the m links (on v 1 , v 2 , . . . , v t ) has energy concentrated on the first several axes. This is exactly what PCA provides. PCA is applied on B T . Then, an anomalous subspace is defined and anomalies that have projections in the anomalous subspace are extracted. In this way, we obtain a left multiplying transformation matrix T, i.e., {tilde over (B)}=TB T . Taking transpose on both side of the equation, we have {tilde over (B)}=({tilde over (B)} T ) T =(TB T ) T =BT T , where T T is a right multiplying transformation matrix that extracts anomalies from B. Note that unlike ARIMA, FFT, or wavelet based tomography, temporal PCA have the property that random reordering of the series, b 1 , b 2 , . . . , b t , does not affect the outcome of the algorithm.

Once the matrix of link anomalies {tilde over (B)} is obtained (step 206 in FIG. 2 ), the next step is to reconstruct OD flow anomalies {tilde over (X)} by solving a series of ill-posed linear inverse problems {tilde over (b)} j =A{tilde over (x)} j . Three embodiments of inference algorithms for solving these problems are described below. All three algorithms deal with the under-constrained linear system by searching for a solution that minimizes some notions of vector norm, three examples of which are:

1) The l 2 norm of a vector v is defined as

 v  2 = ( Σ i ⁢ v i 2 ) 1 / 2 ,

where v i is the i-th element of vector v.

2) The l 1 norm of a vector v is defined as ∥v∥ 1 =Σ i |v i |, i.e., the sum of the absolute value of each element of v. 3) The l 0 norm of a vector V is defined as

 v  0 = Σ i ⁢ v i 0 ,

i.e., the number of non-zero elements of v.

One solution to {tilde over (b)}=A{tilde over (x)} is the pseudoinverse solution {tilde over (x)}=A + {tilde over (b)}, where A + is the pseudoinverse (or Moore-Penrose inverse) of matrix A. It is known that {tilde over (x)}=A + {tilde over (b)} is the solution to the problem {tilde over (b)}=A{tilde over (x)} that minimizes the l 2 norm of the anomaly vector, i.e., it solves:

minimize ∥{tilde over (x)}∥ 2 subject to ∥{tilde over (b)}−A{tilde over (x)}∥ 2 is minimal  Eq. 8.

In practice, we expect only a few anomalies at anyone time, so {tilde over (x)} typically has only a small number of large values. Hence it is natural to proceed by maximizing the sparsity of {tilde over (x)}, i.e., solving the following l 0 norm minimization problem:

minimize ∥{tilde over (x)}∥ 0 subject to {tilde over (b)}=A{tilde over (x)}   Eq. 9.

The l 0 norm is not convex and is notoriously difficult to minimize, so in practice one needs to either approximate the l 0 norm with a convex function or use heuristics, for example, a greedy algorithm.

Another approach to approximate l 0 norm minimization is to convexify Equation 9 by replacing the l 0 norm with an l 1 norm. It can be shown that l 1 norm minimization results in the sparsest solution for many large under-determined linear systems.

minimize ∥{tilde over (x)}∥ 1 subject to {tilde over (b)}=A{tilde over (x)}   Eq. 10.

In the presence of measurement noise, the constraints {tilde over (b)}=A{tilde over (x)} may not always be satisfiable. In this case, we can add a penalty term ∥{tilde over (b)}−A{tilde over (x)}∥ 1 to the objective and reformulate Equation 10 as:

minimize λ∥{tilde over (x)}∥ 1 +∥{tilde over (b)}−A{tilde over (x)}∥ 1   Eq. 11,

where λ∈[0,1] controls the degree to which the constraints {tilde over (b)}=A{tilde over (x)} are satisfied. For example, λ=0.001.

Equation 11 can be cast into the following equivalent Linear Programming (LP) problem, for which solutions are available even when A is very large, owing to modern interior-point linear programming methods:

minimize λΣ i u i +Σ j v j

{tilde over (b)}=A{tilde over (x)}+z

subject to u≧{tilde over (x)}, u≧−{tilde over (x)}

v≧z, v≧−z   Eq. 12.

Note that it is common to use ∥{tilde over (b)}−A{tilde over (x)}∥ 2 2 (instead of ∥{tilde over (b)}−A{tilde over (x)}∥ 1 ) as the penalty term in Equation 11. This alternative formulation can be efficiently solved using methods like Iterative Reweighted Least Squares and has been successfully applied to recover sparse overcomplete representations in the presence of noise.

Another heuristic solution to l 0 norm minimization is to apply the greedy algorithm. For example, the greedy heuristic has been successfully applied to wavelet decomposition, where it goes by the name of Orthogonal Matching Pursuit (OMP). Following a similar technique, a greedy solution can be used to maximize the sparsity of {tilde over (x)}. The algorithm starts with an empty set I of non-zero positions for {tilde over (x)} and then iteratively add new non-zero positions to I. During each iteration, for each position p∉I, the algorithm tests how much it can reduce the residual {tilde over (b)}−A{tilde over (x)} by including p as a non-zero position. More specifically, let J=I∪{p}. The algorithm estimates the values for the non-zero elements of {tilde over (x)} (denoted as {tilde over (x)} J ) by solving the following least squares problem:

minimize ∥{tilde over (b)}−A J {tilde over (x)} J ∥ 2   Eq. 13,

where A J =A[.,J] is a submatrix of A formed by the column vectors of A corresponding to positions in J. The residual is then computed as e J =∥{tilde over (b)}−A J {tilde over (x)} J ∥ 2 . The algorithm then greedily chooses the position p that gives the smallest e J and add it to I. The algorithm stops whenever either the residual energy falls below some tolerable inaccuracy e max or the number of non-zero positions exceeds some threshold l max 0 .

An exemplary pseudo-code representation of the greedy algorithm described above is as follows:

In the above embodiments, it is assumed that the routing matrices are constant. In another embodiment, we wish to allow for dynamic routing changes, and so we must allow A j to vary over time. In addition, if some measurements are missing (say at time j), we may still form a consistent problem by setting the appropriate row of A j to zero. Thus, for realistic SNMP measurements where missing data are often an issue, we still wish to vary A j even for static routing. Note that we assume that the routing measurements are obtained via a route monitor which gives us accurate, up-to-date measurements of routing (at least at the time scale of SNMP measurements, e.g. minutes).

When the routing matrix is non-constant, there is no reason to believe that the measurements B should follow a simple model such as an ARIMA model. Even where the traffic itself follows such a model, a simple routing change may change a link load measurement by 100%, for instance by routing traffic completely away from a particular link. If we were to apply the ARIMA model to the measurements B, we would see such a change in routing as a level-shift anomaly. However, its cause is not an unknown change in X (to be discovered), but rather a known change in the routing matrices A j . Likewise, it no longer makes sense to try to exploit spatial correlations which arose from a particular routing, to the case of another routing.

However, it is no less reasonable to approximate the traffic matrix X by an ARIMA model (than B when the routing is constant), even when routing may change. Under such a modeling assumption, we can write {tilde over (X)}=XT. We know also that the measurements are given by Equation 2. An approach to the solution is therefore to seek a solution {tilde over (X)} that is consistent with these equations, but also minimizes one of the norms (described above) at each time step. We choose to minimize the l 1 norm ∥{tilde over (x)} j ∥ 1 here because: (i) it allows us to naturally incorporate link load constraints at multiple time intervals; and (ii) it is more accurate than both the pseudoinverse and the greedy algorithms for static anomography.

Unfortunately, for transform based methods (the Fourier, wavelet, and PCA methods) the number of constraints becomes very large (as t grows). On the other hand, the set of constraints for the ARIMA model can be written in a form such that it does not grow with t. Hence, in the following we concentrate on generalizing the ARiMA approach. We first present the basic algorithm for ARIMA(p, d, q) models with d≧1. To improve its efficiency, we develop two simple techniques that significantly reduce the problem size. Those skilled in the art will appreciate that the algorithm can also be extended to handle ARIMA models with d=0. We will also discuss model selection and parameter estimation, two important issues for applying ARIMA-based anomography.

We are going to seek solutions that are consistent with the measurements b j =A j x j , for j=1, . . . , t, and an ARIMA model that gives {tilde over (X)}=XT, where T is the same transformation matrix implicitly defined by Equation 6 and Equation 7. Importantly, we do not wish to have to estimate X (or we may as well use an early-inverse method). The advantage of the ARIMA model is we do not need to know X, but only linear combinations of X.

Let L be the backshift operator, whose effect on a process z can be summarized as (Lz) k =z k−1 . Let the AR polynomial Φ(L) be:

Φ ⁡ ( L ) = ∑ i = 0 d + p ⁢ γ i ⁢ L i ⁢ = def ⁢ ( 1 - ∑ i = 1 p ⁢ ϕ i ⁢ L i ) ⁢ ( 1 - L ) d .

Let y k−i =γ i x k−i . We now identify e={tilde over (x)} in the ARIMA model described in Equation 5 (or rather its multivariate extension). By definition the sum

∑ i = 0 d + p ⁢ y k - i = z k - ∑ i = 1 p ⁢ ϕ i ⁢ z k - i

and so, for d≧1, the ARIMA model (Equation 5) can be rewritten:

∑ i = 0 d + p ⁢ y k - i = x ~ k - ∑ j = 1 q ⁢ θ j ⁢ x ~ k - j . Eq . ⁢ 14

Define c k−i =γ i b k−i , then as y k−i =γ i x k−i , the measurement equation (2) implies

A k−i y k−i =c k−i , i= 0,1 , . . . , d+p   Eq. 15.

We can compute {tilde over (x)} 1 ,{tilde over (x)} 2 , . . . , {tilde over (x)} t iteratively by solving a series of l 1 norm minimization problems P k (k=1,2, . . . , t):

P k :minimize ∥{tilde over (x)} k ∥ 1 subject to Equations 14 and 15  Eq. 16.

As an illustrative example, consider the simplest ARIMA model, ARIMA(0,1,0). In this case, p=q=0, so we can write:

Φ ⁡ ( L ) = ∑ i = 0 1 ⁢ γ i ⁢ L i = ( 1 - L ) ,

so γ 0 =1, γ 1 =−1, and Equation 14 becomes

x ~ k = ∑ i = 0 1 ⁢ y k - i ,

thus problem P k is:

minimize ∥{tilde over (x)} k ∥ 1

We apply zero-padding to handle the initial condition when k≦q or k≦d+p. Specifically, for the MA part of the model, we set {tilde over (x)} k−j to 0 whenever k≦j. For the AR part of the model, we apply zero-padding on the differenced series {(1−L) d x i }, which can be achieved by redefining the AR polynomial Φ(L)=Σ i γ i L i as

Similar to above, we can accommodate measurement noise by incorporating penalty terms into the objective to penalize against violation of constraints (Eq. 14) and (Eq. 15). We can then solve the resulting l 1 norm minimization problem by reformulating it as an equivalent LP problem.

One potential problem with the above algorithm is its high computational cost. Even though the computational cost is fixed relative to t, it is still highly dependent on the number of traffic matrix elements n, and the order parameters (p, d, q), and so the l 1 minimization problem P k can be very large even when A k−i stays constant. In contrast, the static anomography algorithm had quite good computation properties for constant A j . Below, two techniques are described to significantly reduce the size of P k . These techniques are result by the following observations: (i) the routing matrices are often quite stable and tend not to change in every time interval; and (ii) when the changes occur, they tend to be local changes and most rows of the routing matrix will remain the same.

The first technique reduces the problem size by merging constraints for intervals with the same routing matrix. Specifically, if there exist i 1 <i 2 such that A k−i 1 =A k−i 2 , we can use a single unknown vector y′ k−i 1 to represent y k−i 2 +y k−i 2 in Equation 14 and then replace the two sets of constraints on y k−i 1 and y k−i 2 in Equation 14 with a single set of constraints on y′ k−i 1 :

A k - i 1 ⁢ y k - i 1 = c k - i 1 ⁢ = def ⁢ c k - i 1 + c k - i 2

We can repeat this process until all A k−i are distinct.

The second technique exploits the fact that there is often a large subset of rows common to all A k−i . Before describing the technique, we first introduce some notations. Let S be a set of integers. Given a matrix M, let M S be the sub matrix of M that consists of all rows with indices in S. Similarly, given a vector v, let v S be the sub-vector of v that consists of elements with indices in S.

Using the above notations, let C be the set of row indices such that all A k−i C are equal (denoted by A C ). Let C be the set of row indices not in C. We can then decompose Equation 15 into

Summing up all the constraints in Equation 18 over i, we obtain:

A C ⁢ ∑ i = 0 d + p ⁢ y k - i = ∑ i = 0 d + p ⁢ c k - i C . Eq . ⁢ 20

Combining Equation 20 and Equation 14, we get a single set of constraints on {tilde over (x)} k :

∑ i = 0 d + p ⁢ c k - i C = A C ⁢ x ~ k - A C ⁢ ∑ j = 1 q ⁢ θ j ⁢ x ~ k - j . Eq . ⁢ 21

We can then replace Equation 18 with Equation 21, reducing P k to:

minimize ∥{tilde over (x)} k ∥ 1 subject to Eqs. 14, 19, and 21  Eq. 22.

In the special case when A j is time invariant, we have C=0. So the constraints y k−i (Eq. 19) become empty, causing y k−i to become free variables in Equation 22. In this case, we can further simplify Equation 22 by eliminating y k−i and the corresponding constraint Equation 14, resulting in

minimize ∥{tilde over (x)} k ∥ 1 subject to Equation 21  Eq. 23.

It can be verified that Equation 23 is equivalent to our original algorithm for time invariant A j , which first computes {tilde over (B)} by applying ARIMA modeling on B and then estimates {tilde over (x)} k by minimizing ∥{tilde over (x)} k ∥ 1 subject to {tilde over (b)} k =A{tilde over (x)} k . This is appealing in that we now have one unified algorithm for ARIMA-based anomography.

Below, we address two important issues on applying ARIMA type models: model selection (i.e., determining p, d, q) and parameter estimation (i.e., determining coefficients φ i ,θ j ). We do so by answering the following three questions. What data to use; How to estimate φ i ,θ j given (p, d, q); and How to select the model order (p,d, q).

The first question is what data to use for model selection and parameter estimation. This is important because the data determines the model and the parameters that we can learn. In the context of network anomography, we would like the model and parameters to capture the normal behavior of x i . For this purpose, we propose to select our models and parameters based on the traffic aggregated at different ingress points. More specifically, let I be the set of indices corresponding to all the ingress points in the link load vectors b i . We will use the series of sub-vectors b i I as the input data for model selection and parameter estimation.

This has several advantages. First, b i I is readily available and does not require any inference. More importantly, ingress traffic is largely invariant to internal topology and routing changes in the local domain under consideration, making our algorithms applicable even in the presence of topology and routing changes. In addition, each element of b i I aggregates a number of OD flows. So the effect of anomalies or missing data on individual OD flows is less significant, making our results more robust.

To estimate φ i ,θ j , given (p, d, q) and input vector series {b k I }, we can estimate the autoregression and moving-average coefficients φ i ,θ j by constructing a state-space model and then applying the standard Kalman filter adaptation.

To select the model order, we first select the degree of differencing (d). The optimal degree of differencing is often the one at which the standard deviation of the differenced series is the lowest. We can apply this rule to determine d. More specifically, for each d∈{0,1,2,3,4}, we compute the differenced series

Z d =[z d,1 z d,2 . . . z d,t ]=(1− L ) d [b 1 I b 2 I . . . b t I ]

Let

E ⁡ [ Z d ] = 1 t ⁢ ∑ i = 1 t ⁢ z d , i

and

V ⁡ [ Z d ] = 1 t ⁢ ∑ i = 1 t ⁢  z d , i - E ⁡ [ Z d ]  2 2 .

We then pick the d that results in the minimum variance V[Z d ].

Once we have d, we can automate the choice of p and q by applying an information based criterion such as the AIC (Akaike's Information Criterion) or AICC. Information based criteria are designed to achieve a good balance between model parsimony and low prediction error.

FIG. 3 is a block diagram depicting an exemplary embodiment of a computer 300 suitable for implementing the processes and methods described herein. Notably, computer 300 may be used to implement the network-level anomaly analyzer 116 of FIG. 1 and the method 200 of FIG. 2 . The computer 300 includes a central processing unit (CPU) 301 , a memory 303 , various support circuits 304 , and an I/O interface 302 . The CPU 301 may be any type of microprocessor known in the art. The support circuits 304 for the CPU 301 include conventional cache, power supplies, clock circuits, data registers, I/O interfaces, and the like. The I/O interface 302 may be directly coupled to the memory 303 or coupled through the CPU 301 . The I/O interface 302 may be coupled to various input devices 312 and output devices 311 , such as a conventional keyboard, mouse, printer, and the like.

The memory 303 may store all or portions of one or more programs and/or data to implement the processes and methods described herein. Notably, the memory 303 may store program code to be executed by the CPU 301 for implementing the network-level anomaly analyzer 116 of FIG. 1 and performing the method 200 of FIG. 2 . Although one or more aspects of the invention are disclosed as being implemented as a computer executing a software program, those skilled in the art will appreciate that the invention may be implemented in hardware, software, or a combination of hardware and software. Such implementations may include a number of processors independently executing various programs and dedicated hardware, such as ASICs.

The computer 300 may be programmed with an operating system, which may be OS/2, Java Virtual Machine, Linux, Solaris, Unix, Windows, Windows95, Windows98, Windows NT, and Windows2000, WindowsME, and WindowsXP, among other known platforms. At least a portion of an operating system may be disposed in the memory 303 . The memory 303 may include one or more of the following random access memory, read only memory, magneto-resistive read/write memory, optical read/write memory, cache memory, magnetic read/write memory, and the like, as well as signal-bearing media as described below.

An aspect of the invention is implemented as a program product for use with a computer system. Program(s) of the program product defines functions of embodiments and can be contained on a variety of signal-bearing media, which include, but are not limited to: (i) information permanently stored on non-writable storage media (e.g., read-only memory devices within a computer such as CD-ROM or DVD-ROM disks readable by a CD-ROM drive or a DVD drive); (ii) alterable information stored on writable storage media (e.g., floppy disks within a diskette drive or hard-disk drive or read/writable CD or read/writable DVD); or (iii) information conveyed to a computer by a communications medium, such as through a computer or telephone network, including wireless communications. The latter embodiment specifically includes information downloaded from the Internet and other networks. Such signal-bearing media, when carrying computer-readable instructions that direct functions of the invention, represent embodiments of the invention.

While the foregoing is directed to embodiments of the present invention, other and further embodiments of the invention may be devised without departing from the basic scope thereof, and the scope thereof is determined by the claims that follow.