Traffic Rate Network Tomography with Higher-Order Cumulants

Network tomography was formulated by Y. Vardi in the seminal paper [1]. The goal is to estimate the rates of traffic flows over source-destination pairs from traffic flows over links of the networks. Traffic flow rates are measured, for example, by the number of packets per second. We use $X$ to denote a vector of $L$ source-destination traffic flows, and $Y$ to denote a vector of $M$ link traffic flows. Normally, $M\ll L$. All source-destination traffic flows are assumed independent Poisson random variables. Thus, $X$ comprises a vector of $L$ independent Poisson random variables. Link traffic flow measurements are assumed passive and require no probes. For networks operating under a deterministic routing regime, traffic flows from each source node is routed to the destination node over a fixed path. Traffic flows over each link may originate from multiple source-destination traffic flows. Thus, we define a binary variable $a_{ij}$ such that $a_{ij}=1$ when traffic over source-destination $j$ passes through link $i$ and $a_{ij}=0$ otherwise. We also define an $M\times L$ binary routing matrix $A=\{a_{ij},i=1,\ldots,M;j=1,\ldots,L\}$ which is assumed known throughout the paper. It follows that $Y=AX$. The expected value of the $j$th Poisson source-destination traffic flow constitutes its traffic flow rate. Thus, the vector of source-destination traffic flow rates is the expected value of $X$ which we denote by $\lambda=E\{X\}$.

The network tomography problem is that of estimating the rate vector $\lambda$ from given realizations of $Y$. This is a rather challenging inverse problem since $Y=AX$ is an underdetermined set of equations, and the estimate of $\lambda$ must be non-negative. Vardi invoked the terminology “LININPOS” for (LINear INverse POSitive problem). Maximum likelihood estimation of $\lambda$ is not feasible since the components of $Y$ are dependent Poisson random variables with no explicitly known distribution. The expectation-maximization (EM) algorithm is also not useful for this problem as it requires calculation of the conditional mean $E\{X\mid Y\}$ which in turns requires knowledge of the joint distribution of $X$ and $Y$. Vanderbei and Iannone [2] developed an EM algorithm in which $E\{X\mid Y\}$ is simulated. A common approach is to rely on the explicit form of $E\{X\mid Y\}$ for jointly Gaussian $X$ and $Y$, see, e.g., [1], [3], [4].

Vardi resorted to moment matching in which $\lambda$ is estimated from a linear matrix equation relating the first two empirical moments of $Y$ to the corresponding first two theoretical moments of $AX$. The approach is applicable to networks operating under deterministic as well as random routing regimes. In the latter case, there are multiple alternative paths for each source-destination pair which can be selected according to some probability law. Vardi invoked an iterative procedure for estimating $\lambda$ which was previously developed for image deblurring [5]. This useful procedure will be discussed in Section II. Vardi’s work ignited extensive research in the areas of network inference and medical tomography.

In this paper we explore the benefits of supplementing Vardi’s second-order moment matching approach with higher-order empirical cumulants. Here too the cumulant matching equation is linear in $\lambda$, and the approach is applicable to networks operating under deterministic as well as random routing regimes. Higher-order cumulants introduce new useful information on the unavailable distribution of $Y$ which can be leveraged in estimating the source-destination flow rates. By using a sufficient number of empirical cumulants, the linear mapping involved in the cumulant matching approach can achieve full column rank. In an ideal scenario where the cumulant matching equations are consistent, and the cumulants of the link traffic flows are accurately known, the rate vector can be recovered without error. As we demonstrate in Section V, this ideal situation is approachable when a sufficiently large number of realizations of $Y$ are available.

A vast literature exists on network tomography. Several forms of network tomography have been studied in the literature depending on the type of measurement data and the network performance parameters (e.g., rate, delay) of interest: (i) source-destination path-level traffic rate estimation based on link-level traffic measurements [1, 6, 3, 4, 7, 8, 9]; (ii) link-level parameter estimation based on end-to-end path-level traffic measurements [2], [10, 11, 12, 13, 14]; (iii) network topology discovery from traffic measurements [15, 16, 17, 18]. Approaches to network tomography can also be characterized as active (e.g., [10, 11, 19, 20]), whereby explicit control traffic is injected into the network to extract measurements, or passive (e.g., [1, 3, 21, 14]), whereby the observation data is obtained from existing network traffic. A somewhat outdated survey circa 2003 can be found in [22]. A more recent survey [23] discusses network tomography in conjunction with network coding.

Closest to Vardi’s work is the contemporary work of Vanderbei and Iannone [2] which relies on a Poisson model for incoming traffic. The goal is to estimate the rate of traffic on each link connecting input and output nodes from traffic counts at the input and output nodes. Vanderbei and Iannone did not resort to moment matching but rather developed a simulation based EM algorithm for maximum likelihood estimation of the rates. A thorough Bayesian approach to Vardi’s problem was developed by Tebaldi and West [6] using Markov Chain Monte Carlo simulation. See also [9]. Another closely related work to Vardi’s problem appeared in [3], [4], [22] where maximum likelihood estimation of the source-destination rates from link data was implemented under a Gaussian rather than a Poisson traffic model. In [7], source-destination rates were estimated by utilizing spatiotemporal correlation of nominal traffic, and the fact that traffic anomalies are sparse. In [8], conditions for identifiability of higher order cumulants in estimation of source-destination traffic from link measurements were established. In [13] an algorithm was developed for choosing a set of linearly independent source-destination measurement paths from which additive link metrics are individually estimated.

The plan for this paper is as follows. In Section II we present the minimum I-divergence iterative procedure which plays an important role in this paper. In Section III we address rate estimation in networks with deterministic routing. We develop a new systematic set of cumulant matching equations for estimating the rate vector from up to the fourth-order empirical cumulant. We also develop an upper bound on the mean squared error (MSE) in least squares estimation of the rate vector from empirical cumulants. We conclude this section with a discussion on some implementation issues of the proposed cumulant matching approach. Complexity of the approach is discussed in Section III-D. In Section IV we discuss rate estimation in networks operating under a random routing regime. Numerical results for the NSFnet [24] are presented in Section V. Concluding remarks are given in Section VI. Technical details of the derivations are deferred to the Appendix.

Under the Poisson model for $X$, the moment (and cumulant) matching approach results in a set of linear equations in $\lambda$. Generally speaking, the set of equations has the form of $\hat{\eta}(Y)=B\lambda$ where $\hat{\eta}(Y)$ is a vector of $M_{b}$ empirical moments (first and second-order in [1]) of $Y$, $B=\{b_{ij},i=1,\ldots,M_{b};j=1,\ldots,L\}$ is a constant zero-one matrix that depends on $A$ but not on $\lambda$, and $B\lambda$ represents the corresponding theoretical moments.

To estimate $\lambda$ that satisfies $\hat{\eta}(Y)=B\lambda$, Vardi proposed to use an iterative estimation approach which originated in the study of another inverse problem concerning image deblurring [5], [25]. Let $\lambda_{j}^{\textrm{old}}$ denote a current estimate of the $j$th component of $\lambda$, and let $\lambda_{j}^{\textrm{new}}$ denote the new estimate of that component at the conclusion of the iteration. Let $(B\lambda^{\textrm{old}})_{i}$ denote the $i$th component of $B\lambda^{\textrm{old}}$. Similarly, let $\hat{\eta}_{i}(Y)$ denote the $i$th component of $\hat{\eta}(Y)$. The iteration is given by

for $j=1,\ldots,L$. This iteration is particularly suitable for solving positive inverse problems. It reaches a fixed point when the moment matching equation $\hat{\eta}(Y)=B\lambda^{\textrm{old}}$ is satisfied. The iteration was studied by Snyder, Schulz and O’Sullivan [25] in a similarly formulated application of image deblurring. It was shown to monotonically decrease Csiszár’s I-divergence [26] between the original image convolved with the kernel, and the observed blurred image. The procedure turned out to be an EM iteration in the positron emission tomography problem, which follows a similar formulation as network tomography, but with the crucially facilitating difference that the Poisson components of $Y$ are now independent random variables [27], [28]. This iteration will play a central role in our cumulant matching approach.

In this section we present our cumulant matching approach and its theoretical performance bound for networks operating under a deterministic routing regime. The goal is to estimate the source-destination rate vector $\lambda$ from $N$ realizations of the link traffic flow $Y$ where $Y=AX$ and the routing matrix $A$ is known. Conditions for identifiability of $\lambda$ were given in [1]. Specifically, the parameter $\lambda$ is identifiable if all columns $\{a_{1},\ldots,a_{L}\}$ of $A$ are distinct and none is null.

In a typical network, there are multiple paths connecting every source node with every destination node. When the network operates under a deterministic routing regime, a single fixed path is used for every source-destination pair. When the network operates under a random routing regime, a path from the source to destination nodes is selected according to some probability law. Vardi attributed Markovian routing for traffic on each source-destination pair. The accessible nodes and links for the given pair are represented by the states and transition probabilities of the Markov chain, respectively. The transition probabilities of the Markov chain determine the probability of each path with the same source-destination address.

Tebaldi and West argued that random routing can be viewed as deterministic routing in a super-network in which all possible paths for each source-destination address are listed [6]. This approach results in an expanded zero-one routing matrix with each column in the original routing matrix replaced by multiple columns representing the feasible paths for the source-destination pair. The Poisson traffic flow on a source-destination pair is thinned into multiple Poisson traffic flows with the multinomial thinning probabilities being the probabilities of the paths with the same source-destination address. Thus, the super network and the original network operate under similar statistical (Poisson) models. The thinned Poisson rates can now be estimated as was originally done, e.g., using cumulant matching, and each source-destination rate estimate can be obtained from the sum of the thinned rate estimates in that source-destination pair.

In this section we demonstrate the performance of our approach and the gain realized by using higher-order empirical cumulants²²2Our experiments were run on ARGO, a research computing cluster provided by the Office of Research Computing at George Mason University, VA. (URL: http://orc.gmu.edu).. We study the NSFnet [24] whose topology is shown in Fig. 1. The network consists of 14 nodes and 21 bidirectional links. Hence, it contains $L=14\cdot 13/2=91$ source-destination pairs.

This size network may represent a private network, a transportation network or a subset of interest of a large-scale network. The link weights in Fig. 1 are exclusively used to determine $k\geq 1$ shortest paths for each source-destination pair. Otherwise, they play no role in the traffic rate estimation problem. To determine the $k$ shortest paths between a given source-destination pair, we used the shortest simple paths function from the NetworkX Python library, which is based on the algorithm of Yen [34]. When $k=1$, the number of source-destination paths equals the number of source-destination pairs and the routing matrix $A$ is a $21\times 91$ matrix. The augmented routing matrix $\mathcal{A}_{r}$ achieves full column rank when $r=2$.

With $k\geq 2$, we can assign multiple paths to each source-destination pair and treat them as distinguishable new source-destination pairs. The routing matrix $A$ thus becomes fatter and using higher-order empirical cumulants becomes beneficial. For example, when $k=2$, we have $L=182$ source-destination paths, the column rank of $\mathcal{A}_{2}$ is $162$ and $\mathcal{A}_{3}$ has full column rank. Thus, using this example we focus on third-order cumulant matching.

The network with $k=2$ and a $21\times 182$ routing matrix $A$ could also be seen as a super-network in the Tebaldi-West sense [6] for a network with $L=91$ source-destination pairs operating under a random routing regime with two possible paths per each source-destination pair. The accuracy of the rate estimation for the random routing network is determined by the accuracy of the rate estimation in the deterministic routing super-network. Thus, it suffices to focus on rate estimation in the deterministic routing network with the $21\times 182$ routing matrix $A$.

The arrival rates $\{\lambda_{j}\}$ in our experiment were generated randomly from the interval $[0,4]$. In each of $T=500$ simulation runs, a rate vector $\lambda$ was generated, and was subsequently used in generating $N$ statistically independent identically distributed Poisson vectors $\{X_{n}\}$ which were transformed into the vectors $\{Y_{n}=AX_{n}\}$ using the assumed known routing matrix $A$. The $N$ statistically independent identically distributed Poisson vectors $\{Y_{n}\}$ were used to generate the empirical cumulants using (29) when $r=1$ and (30) when $r=2,3$. We experimented with $N$ in the range of $10\,000$ to $500\,000$. The MSE in estimating the cumulants for the NSFnet with $M=21$ are given in Table II. Clearly, the MSE decreases monotonically with $N$ for $r=1,2,3$.

The empirical cumulants were used to estimate the rate vector in the current run. For the cumulants regularization we have used $\epsilon_{2}=1$ and $\epsilon_{3}=.01$. Let $\lambda_{t}(i)$ and $\hat{\lambda}_{t}(i)$ denote, respectively, the $i$th component of $\lambda$ and its estimate at the $t$th run where $i=1,\ldots,L$ and $t=1,\ldots,T$. For each estimate we evaluated the normalized MSE defined by

and the averaged normalized MSE defined by

The MSE in estimating $\lambda_{i}$ is approximately $\xi_{i}^{2}\cdot E\{\lambda^{2}(i)\}$ when $T$ is sufficiently large.

Two rate estimators were used, the iterative estimator (1) and the least squares estimator (60). The iteration was initialized uniformly with all rates set to $.1$. It was terminated after $300$ iterations. The least squares regularization factor was set to $\gamma=.0005$.

The performance of the iterative estimator (1) is demonstrated in Table III. The table shows the averaged normalized MSE $\overline{\xi^{2}}$ in estimating the path rates using $r=1,2,3$ empirical cumulants estimated from $N$ vectors $\{Y_{n}\}$ with $N$ ranging from $N=10\,000$ to $N=500\,000$. The averaged normalized MSE $\overline{\xi^{2}}$ in estimating the path rates using $r=1,2,3$ theoretical cumulants is shown as well. As can be observed in the table, $\overline{\xi^{2}}$ monotonically decreases with $N$ for $r=2,3$. For a given $N$, the performance using $r=3$ cumulants is better than using $r=2$ moments, except for $N=10\,000$ where the difference is only $.0005$. Similar behavior can be seen in Table IV when the least squares estimator (60) is used for the path rate estimation. The percentage of negative rate estimates observed with the least squares estimator is shown in Table V. This percentage is negligible for all values of $r$ and $N$ and does not exceed $5\%$. Comparing the two estimators (1) and (60), the least squares estimator appears superior as it consistently provides lower $\overline{\xi^{2}}$ at the expense of negligible percentage of negative rate estimates.

To assess the effectiveness of the third-order cumulant matching approach and compare it with Vardi’s second-order moment matching approach, we contrast the performance of the former approach with the best performance of the latter approach as obtained by using the theoretical rather than the empirical moments. Let $\overline{\xi^{2}}(r,N)$ denote the averaged normalized MSE from (62) obtained when $\lambda$ is estimated using $r$ empirical cumulants that were estimated from $N$ realizations of the vector $Y$. Let $\overline{\xi^{2}}(r,\infty)$ denote the averaged normalized MSE $\overline{\xi^{2}}$ from (62) obtained by using $r$ theoretical cumulants, or effectively by using $N=\infty$. Table VI shows $\overline{\xi^{2}}(2,N)/\overline{\xi^{2}}(2,\infty)$ and $\overline{\xi^{2}}(3,N)/\overline{\xi^{2}}(2,\infty)$ for $N=200\,000$ and the estimators (1) and (60) discussed in Tables III–IV.