# Estimating Network Flow Length Distributions via Bayesian Nonnegative Tensor Factorization

Authors: Barış Kurt; Ali Taylan Cemgil; Güneş Karabulut Kurt; Engin Zeydan
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

In this paper, we develop a framework to estimate network flow length distributions in terms of the number of packets. We model the network flow length data as a three-way array with day-of-week, hour-of-day, and flow length as entities where we observe a count. In a high-speed network, only a sampled version of such an array can be observed and reconstructing the true flow statistics from fewer observations becomes a computational problem. We formulate the sampling process as matrix multiplication so that any sampling method can be used in our framework as long as its sampling probabilities are written in matrix form. We demonstrate our framework on a high-volume real-world data set collected from a mobile network provider with a random packet sampling and a flow-based packet sampling methods. We show that modeling the network data as a tensor improves estimations of the true flow length histogram in both sampling methods.


---

## Body

## 1. Introduction

Monitoring network statistics is crucial for the maintenance and infrastructure planning for network service providers. Statistical information about traffic patterns helps a service provider to characterize its network resource usage and user behavior, to infer future traffic demands, to detect traffic and usage anomalies, and to provide insights to improve the performance of the network [1]. However, network monitoring has become a difficult task due to increasingly high-volume and high-speed data over modern networks, and in most cases, it requires special hardware. For this reason, sampling [2] becomes a viable approach for extracting statistics from such high-speed networks. In this work, we are interested in one of the most important network statistics, the flow length distribution (FLD).A network flow is defined as a set of Internet protocol (IP) packets with the same signature observed within a limited time period. The flow signature is composed of the IP and port pairs of both the source and destination nodes together with level-3 protocol types such as transport control protocol (TCP) or user datagram protocol (UDP). A flow starts with the arrival of the first packet and terminated when the interpacket timeout is exceeded. The total number of packets in a flow is referred to as the flow length and the length distribution of a set of flows that are terminated in a time window is called flow length distribution.In this work, we are using one of the most popular methods for collecting per-flow information, i.e., passive measurement. In this method, network packets are processed as they pass through a passive measurement beacon connected to the network, e.g, router. The beacon keeps a look-up table for flow identification. The beacon processes a packet by searching its corresponding flow inside the look-up table using its signature. If such a flow is found, its statistics are updated. Otherwise, the packet is treated as the first packet of a new flow, and the new flow is inserted into the table. Once a flow is terminated, its statistics are transferred to a storage.The flow length histogram can be calculated exactly by processing every packet that passes through the measurement beacon. In order to implement such a direct method, the monitoring beacon needs to maintain a table to hold information for all active flows on the network. However, substantial amount of concurrent flows with very short packet interarrival times of current high-speed networks (on the order of 10 Gbps to 100 Gbps inside carrier’s network today) make this brute-force counting method very costly to implement. First of all, this method would require a large amount of memory to record the flow table. Secondly, in a high-speed link, the interarrival times between packets, which may be as small as 8 nanoseconds in an OC-768 link, may be smaller than the time required to process flow hash operations such as identifying the packet and updating the flow statistics.The characteristics of the network traffic data inevitably lead to the development of alternative methods for measurement such as random sampling, where a fraction of the network traffic is randomly selected and processed. The simplest sampling method is the uniform packet sampling [3–6], used in commercial systems [7, 8]. In uniform sampling, each packet is selected with a predefined constant probability. This approach is easy to implement since it does not require the flow identification of each packet. However, recovering the true flow length distribution from the random packet sampled traffic is a challenging problem. The unbiased estimator of the original flow length n for sampling probability π is n^m=m/π, where m is observed flow length. The relative error of this estimator, calculated as 1/π−1 [3], grows unboundedly for short flows as the sampling rate gets smaller. The high error on the small flow lengths comes from the fact that most of the samples are collected from longer flows.Flow-based adaptive sampling methods [9–14] were proposed for more accurate flow length estimation. In these methods, each incoming packet is processed and then sampled with a probability that is a function of the current sampled length of the flow that the packet belongs to. Here, the main idea is to use a smaller memory by compressing the flow statistics counters in the router. However, these methods need to be implemented on specialized and expensive hardware due to the mandatory packet identification and look-up step.Both packet-based and flow-based adaptive sampling methods rely on numerical methods to recover the true FLD. In this work, we propose a framework that can be used to recover the true FLD from the sampled observation obtained by any sampling method. This framework uses a variant of the nonnegative tensor factorization NTF model, which we call the thin nonnegative tensor factorization (ThinNTF), where the “thin” prefix emphasizes that the factorization is applied directly to the samples, or namely “thinned” data.In our framework, the network traffic data is modeled as a 3-way array, containing the number of flow length observations, with dimensions interpreted as (1) flow length, (2) hour-of-day, and (3) day-of-week to capture the hourly and daily periodicity in the data. The nonnegative factorization of this tensor basically gives us estimates in the form of a nonparametric mixture model. Therefore, our model is an improvement of the nonparametric flow length models in [3, 6] by having the capability of modeling data with an arbitrary amount of mixture components and use the periodicity.While the ordinary NTF model [15] factorizes an observation tensor, the ThinNTF directly factorizes its sampled version and recovers the original tensor from the estimated factors. We take a fully Bayesian approach here and provide a generative model for the TNTF and a variational Bayes algorithm for inference. The contributions in this paper can be listed as follows:(i)
We model one week of flow length observations as a 3-dimensional tensor and observe the periodic behavior.(ii)
We propose a novel tensor factorization scheme, ThinNTF, which is able to find the factors of a latent tensor from its sampled counterpart. By doing so, we also solve the reconstruction problem.(iii)
We apply ThinNTF to real-world data sampled with two different sampling methods: uniform random packet sampling and flow-based adaptive sampling.The structure of the paper is as follows. In Section2, we provide the related works on network sampling and tensor factorization. In Section 3, we describe our real-world data and how we visualize it as a tensor. Additionally, we describe the sampling methods that we used to sample the data. In Section 4, we describe our ThinNTF model and the variational Bayes algorithm for estimating the factors. In Section 5, we describe our real-world data collection architecture. In Section 6, we present our synthetic and real-world experiments and results. Finally, in Section 7, we draw our conclusions.


## 2. Related Work

Sampling methods have long been applied to network traffic monitoring. A survey on fundamental network sampling strategies is given in [2]. Uniform packet sampling is extensively studied by various authors. Duffield et al. [3] propose the first nonparametric model for flow length distribution and provide a maximum likelihood estimation to recover the flow lengths. Riberio et al. [4] show that using protocol specific information gives better flow length distribution estimates in TCP flows. Yang and Michailidis [6] adopt the maximum likelihood approach to estimate both flow length and flow volume (number of bytes in a flow) distributions. Additionally, they model the data with a nonparametric mixture model of two components, where the first component models small flows and the second models large ones.Flow-based sampling methods are proposed as alternatives to the uniform packet sampling since packet sampling has theoretical limitations when recovering true flow statistics [5]. These methods process every incoming packet and apply sampling conditionally. Kumar et al. [13, 16] propose two different algorithms where the flow size counters are compressed statistically. They also propose a nonuniform packet sampling algorithm based on sketch counting [12]. Hu et al. [10, 14] propose another nonuniform packet sampling algorithm, called adaptive nonlinear sampling (ANLS), for estimating flow lengths per each flow and then adopt this method to flow volume [11]. In our experiments, we are going to use ANLS as an example of flow-based sampling methods since it is the current state-of-the-art nonuniform sampling method.Nonnegative tensor factorization is the generalization of the nonnegative matrix factorization (NMF) [17] to multiway arrays. In NMF, a nonnegative matrix is approximated with a multiplication of two nonnegative matrices. Minimizing the Kullback–Leibler divergence between the initial matrix and multiplied factors is a popular formulation of this method and can be solved with fixed-point iterations [18] or full Bayesian methods [19]. NMF has been used in many applications such as spectral data analysis [20], face recognition [17], and document clustering [21].Modeling the flow length distribution as a mixture of distributions is first proposed by [6]. However, according to our best knowledge, there is no previous work that models a large volume of flow length data as a tensor. This work fills a gap in the literature by introducing tensor factorization methodology to network monitoring.


## 3. Problem Description

We describe our problem as a tensor thinning problem, where the count entities of the original flow lengths are stored in a tensor. We formulate the sampling process as a matrix multiplication operated on this data tensor. In order to do that, each sampling model should be represented as a matrix that transforms the original data tensor to a sampled one. We provide matrices for two sampling models: uniform packet sampling and ANLS flow-based packet sampling.


### 3.1. Notation and Indexes

For a clear notation, the scalar values are denoted by lightface letters, such as the index variablej and its maximum value J. The vectors are represented by boldface lower case letters, such as vector x. Boldface upper case letters represent matrices, such as F,H, and D, and the tensors are represented with calligraphic upper case letters, i.e, X. The individual entries in matrices and tensors are written like scalars, i.e., fi,r and xi,j,k. The index denotes all the entries in the given dimension. For example si,: is the ith row of the S matrix and Xi,:,: is the ith slice of the tensor X in the first dimension.The index parameters are also fixed for clarity. The list of indexes and their ranges and semantic descriptions are given in Table1. For example, the i index always represents an original flow length, while ν presents a sampled flow length. The range of ν starts from 0, since all of the packets of a flow may be discarded during the sampling process yielding a zero-length sampled flow, which is never observed.Table 1
Indexes in the model.
Index
Range
Description
I
1
,
I
Original flow lengths
N
0
,
I
Sampled flow lengths
J
1
,
J
Hours of day
K
1
,
K
Days
R
1
,
R
Components


### 3.2. Data Tensor

The original flow length data is represented in anI×J×K tensor X, with individual elements xi,j,k, regarded as the number of flows that has length i measured at the hour j of the day k. In this setup, I is the maximum flow length value, J is the hours of the day, and K is the days of the week. For our real-world data, collected continuously for 1 week, these values are I=2,000,000,J=24, and K=7.Working with large maximum flow size is not feasible for two reasons. First, learning a mixture model where each flow component has 2 million parameters is not a good formulation for this problem. Secondly, 99.9% of flows in our data have less than 100 packets, which means the tensorX will be very sparse for i>100. The clamping process can be defined as(1)X¯i,j,k=Xi,j,k,fori<Imax,∑l≥ImaxXl,j,k,fori=Imax,where X¯ is the clamped tensor. The clamping does not require any change in the model and inference equations that are given in Section 4. Therefore, for notational clarity, we only use X as the generic original data tensor.Figure1 shows the vertical slices of our unsampled real-world data tensor X, collected at the backbone of a mobile operator during a one week period, from Monday to Sunday, and clamped at Imax=50. The intensity images are generated from the logarithm of the flow length counts. The daily and hourly patterns are easily recognizable in the original FLD data.Figure 1
Slices of the original flow length tensorX.


### 3.3. Sampling Methods

Independent of the sampling method, we can define anI×I+1 size S matrix, where I is the maximum flow length with entries si,ν interpreted as the probability of sampling ν packets from an original flow of length i. Naturally, S is a lower diagonal matrix of the form:(2)S=s1,0s1,100⋯0s2,0s2,1s2,20⋯0⋮⋮⋮⋮⋯⋮sI,0sI,1sI,2sI,3⋯sI,I,where its lth row defines a probability distribution for the sampled flow length of a flow of size l. Given a flow size distribution x∈ℤ≥I, where ℤ≥ is the set of nonnegative integers, the expected sampled flow length distribution would be given by y^=STx. It immediately follows that the sampled flow size distribution y has length I+1, with y0 being the proportion of sampled flows with none of their packets sampled. During the sampling process, this value will never be observed naturally since the flow identification is performed only on selected packets. In all experiments throughout this paper, the y vector (or Y tensor, which will be described later on) will be element-wise multiplied with a binary mask vector m (or binary mask tensor M), whose entries are set to 1 except the ones corresponding to zero-sampled flow lengths, in order to simulate the real-life scenario.For any given sampling method, we can calculate theS matrix directly if a closed-form expression is available. Otherwise, it can be approximated by simulating the sampling process and counting the sampling statistics. In this paper, both uniform random sampling and the ANLS provide closed-form expressions for the calculation of S matrix.An important practical issue is that, if the original tensorX is clamped at Imax, the S matrix must also be clamped. In that case, a last row entry sImax,ν must present the probability of selecting ν packets from a flow of length greater or equal to Imax. This clamping operation can be done by calculating a full-size S matrix first and setting s¯Imax,ν∝∑i=ImaxIsi,ν with a naive assumption that after Imax, the original flow sizes are uniformly distributed.


#### 3.3.1. Uniform Sampling Method

In uniform sampling, each packet is processed with a fixed probability ofπ, irrespective of the flow it belongs to. In this method, the sampling matrix entries si,ν are calculated directly through Binomial distribution with i trials and π success probability, i.e.,(3)si,ν=iνπν1−πi−ν,forν≤i,0,otherwise.Algorithm1 describes how the flow table is updated with uniform sampling upon the arrival of a new packet. The algorithm uniformly draws a random number in interval 0,1, and if it is less than π, the packet is processed; otherwise, the packet is discarded. For processing the packet, a look-up operation is performed on the flow table to find and update the flow that the packet belongs to. If no such flow is found, a new flow is created using the packet’s signature.Algorithm 1:Uniform packet sampling algorithm.
(1)
function Sample Uniformly Random π, flow_table, packet
(2)
if π > rand_double (0, 1) then
(3)
flow = flow_table·look-up (packet)
(4)
if flow is null then
(5)
flow= new Flow (packet)
(6)
else
(7)
flow·length+ = 1
(8)
flow_table·insert_or_update (flow)Figure2(a) shows the top 10×11 entries of the lower diagonal S matrices with different sampling probabilities. As the sampling probability π gets smaller, fewer packets from a flow gets observed, and the flow may even be missed when none of its packets are observed. The rightmost sampling matrix shows the case when π=1/64, where the matrix has a very high concentration of zero-length sampled flows.Figure 2
Sampling matrices for two different sampling schemes.
(a)
(b)Figure2(b) shows the original Monday data (the leftmost matrix) and its sampled versions under uniform sampling with the probabilities shown on the top sampling matrices. Here, we see that for π=1/64, the observed flow lengths are mostly less than 10, while the majority are not observed at all.


#### 3.3.2. ANLS Sampling Method

The ANLS [10] will be used as the representative of the flow-based adaptive sampling methods. In ANLS, each packet is sampled according to the number of packets previously sampled from its corresponding flow. If a sampled flow has length x, the probability of its next packet to be sampled (px;u) is calculated as(4)fx;u=1+ux−1u,(5)px;u=1fx−1;u−fx;u.Here,fx;u is a monotonically increasing function of flow length x, parametrized with u, which makes px;u monotonically decreasing. u determines the tendency of sampling packets. As u gets smaller, more packets are sampled and estimating original flow lengths gets easier.The ANLS method is described in details in Algorithm2. Prior to the sampling, f and p are calculated in the SetupANLS function, according to equations (4) and (5). During sampling, for each incoming packet, a look-up operation is performed unconditionally. If the corresponding flow is found, it is updated with probability relative to its current observed size. Otherwise, a new flow is created, ensuring that every flow is observed with at least one packet.Algorithm 2:ANLS sampling algorithm.
(1)
function Setup ANLSu
(2)
f[0] = 0
(3)
for i∈1,I do
(4)
fi=1+ui−1/u
(5)
pi=1/fi−1−fi
(6)
return f, p
(7)
function SampleWithANLS (p, flow_table, packet)
(8)
flow = flow_table·look-up (packet)
(9)
if flow is null then
(10)
flow= new Flow (packet)
(11)
else if p[flow·length] > rand_double (0, 1) then
(12)
flow·length+ = 1
(13)
flow_table·insert_or_update (flow)We calculate the sampling matrixS recursively for ANLS. In this method, the first packet is always sampled since p1,u=1 independent of u. We start by assigning all zero sampling probabilities as s:,0=0 and s1,1=1. The recursive equation for calculating the sampling matrix can be deduced as(6)si,ν∝si−1,ν−1pi−1,u+si−1,ν1−pi−1,u.Figure3 shows the ANLS sampling matrices for first 10 flow lengths and the Monday data sampled with them, respectively, similarly to Figure 2. First, we can see that when u is small, the S matrix looks like identity and as u gets larger, the sampling probability of large flows decreases. Secondly, compared to uniform sampling, the ANLS method has much higher sampling ratios than uniform sampling. However, operating with such high sampling ratios would require specialized hardware in real time.Figure 3
Visualization of the Monday slice with uniform and ANLS sampling methods.
(a)
(b)


## 3.1. Notation and Indexes

For a clear notation, the scalar values are denoted by lightface letters, such as the index variablej and its maximum value J. The vectors are represented by boldface lower case letters, such as vector x. Boldface upper case letters represent matrices, such as F,H, and D, and the tensors are represented with calligraphic upper case letters, i.e, X. The individual entries in matrices and tensors are written like scalars, i.e., fi,r and xi,j,k. The index denotes all the entries in the given dimension. For example si,: is the ith row of the S matrix and Xi,:,: is the ith slice of the tensor X in the first dimension.The index parameters are also fixed for clarity. The list of indexes and their ranges and semantic descriptions are given in Table1. For example, the i index always represents an original flow length, while ν presents a sampled flow length. The range of ν starts from 0, since all of the packets of a flow may be discarded during the sampling process yielding a zero-length sampled flow, which is never observed.Table 1
Indexes in the model.
Index
Range
Description
I
1
,
I
Original flow lengths
N
0
,
I
Sampled flow lengths
J
1
,
J
Hours of day
K
1
,
K
Days
R
1
,
R
Components


## 3.2. Data Tensor

The original flow length data is represented in anI×J×K tensor X, with individual elements xi,j,k, regarded as the number of flows that has length i measured at the hour j of the day k. In this setup, I is the maximum flow length value, J is the hours of the day, and K is the days of the week. For our real-world data, collected continuously for 1 week, these values are I=2,000,000,J=24, and K=7.Working with large maximum flow size is not feasible for two reasons. First, learning a mixture model where each flow component has 2 million parameters is not a good formulation for this problem. Secondly, 99.9% of flows in our data have less than 100 packets, which means the tensorX will be very sparse for i>100. The clamping process can be defined as(1)X¯i,j,k=Xi,j,k,fori<Imax,∑l≥ImaxXl,j,k,fori=Imax,where X¯ is the clamped tensor. The clamping does not require any change in the model and inference equations that are given in Section 4. Therefore, for notational clarity, we only use X as the generic original data tensor.Figure1 shows the vertical slices of our unsampled real-world data tensor X, collected at the backbone of a mobile operator during a one week period, from Monday to Sunday, and clamped at Imax=50. The intensity images are generated from the logarithm of the flow length counts. The daily and hourly patterns are easily recognizable in the original FLD data.Figure 1
Slices of the original flow length tensorX.


## 3.3. Sampling Methods

Independent of the sampling method, we can define anI×I+1 size S matrix, where I is the maximum flow length with entries si,ν interpreted as the probability of sampling ν packets from an original flow of length i. Naturally, S is a lower diagonal matrix of the form:(2)S=s1,0s1,100⋯0s2,0s2,1s2,20⋯0⋮⋮⋮⋮⋯⋮sI,0sI,1sI,2sI,3⋯sI,I,where its lth row defines a probability distribution for the sampled flow length of a flow of size l. Given a flow size distribution x∈ℤ≥I, where ℤ≥ is the set of nonnegative integers, the expected sampled flow length distribution would be given by y^=STx. It immediately follows that the sampled flow size distribution y has length I+1, with y0 being the proportion of sampled flows with none of their packets sampled. During the sampling process, this value will never be observed naturally since the flow identification is performed only on selected packets. In all experiments throughout this paper, the y vector (or Y tensor, which will be described later on) will be element-wise multiplied with a binary mask vector m (or binary mask tensor M), whose entries are set to 1 except the ones corresponding to zero-sampled flow lengths, in order to simulate the real-life scenario.For any given sampling method, we can calculate theS matrix directly if a closed-form expression is available. Otherwise, it can be approximated by simulating the sampling process and counting the sampling statistics. In this paper, both uniform random sampling and the ANLS provide closed-form expressions for the calculation of S matrix.An important practical issue is that, if the original tensorX is clamped at Imax, the S matrix must also be clamped. In that case, a last row entry sImax,ν must present the probability of selecting ν packets from a flow of length greater or equal to Imax. This clamping operation can be done by calculating a full-size S matrix first and setting s¯Imax,ν∝∑i=ImaxIsi,ν with a naive assumption that after Imax, the original flow sizes are uniformly distributed.


### 3.3.1. Uniform Sampling Method

In uniform sampling, each packet is processed with a fixed probability ofπ, irrespective of the flow it belongs to. In this method, the sampling matrix entries si,ν are calculated directly through Binomial distribution with i trials and π success probability, i.e.,(3)si,ν=iνπν1−πi−ν,forν≤i,0,otherwise.Algorithm1 describes how the flow table is updated with uniform sampling upon the arrival of a new packet. The algorithm uniformly draws a random number in interval 0,1, and if it is less than π, the packet is processed; otherwise, the packet is discarded. For processing the packet, a look-up operation is performed on the flow table to find and update the flow that the packet belongs to. If no such flow is found, a new flow is created using the packet’s signature.Algorithm 1:Uniform packet sampling algorithm.
(1)
function Sample Uniformly Random π, flow_table, packet
(2)
if π > rand_double (0, 1) then
(3)
flow = flow_table·look-up (packet)
(4)
if flow is null then
(5)
flow= new Flow (packet)
(6)
else
(7)
flow·length+ = 1
(8)
flow_table·insert_or_update (flow)Figure2(a) shows the top 10×11 entries of the lower diagonal S matrices with different sampling probabilities. As the sampling probability π gets smaller, fewer packets from a flow gets observed, and the flow may even be missed when none of its packets are observed. The rightmost sampling matrix shows the case when π=1/64, where the matrix has a very high concentration of zero-length sampled flows.Figure 2
Sampling matrices for two different sampling schemes.
(a)
(b)Figure2(b) shows the original Monday data (the leftmost matrix) and its sampled versions under uniform sampling with the probabilities shown on the top sampling matrices. Here, we see that for π=1/64, the observed flow lengths are mostly less than 10, while the majority are not observed at all.


### 3.3.2. ANLS Sampling Method

The ANLS [10] will be used as the representative of the flow-based adaptive sampling methods. In ANLS, each packet is sampled according to the number of packets previously sampled from its corresponding flow. If a sampled flow has length x, the probability of its next packet to be sampled (px;u) is calculated as(4)fx;u=1+ux−1u,(5)px;u=1fx−1;u−fx;u.Here,fx;u is a monotonically increasing function of flow length x, parametrized with u, which makes px;u monotonically decreasing. u determines the tendency of sampling packets. As u gets smaller, more packets are sampled and estimating original flow lengths gets easier.The ANLS method is described in details in Algorithm2. Prior to the sampling, f and p are calculated in the SetupANLS function, according to equations (4) and (5). During sampling, for each incoming packet, a look-up operation is performed unconditionally. If the corresponding flow is found, it is updated with probability relative to its current observed size. Otherwise, a new flow is created, ensuring that every flow is observed with at least one packet.Algorithm 2:ANLS sampling algorithm.
(1)
function Setup ANLSu
(2)
f[0] = 0
(3)
for i∈1,I do
(4)
fi=1+ui−1/u
(5)
pi=1/fi−1−fi
(6)
return f, p
(7)
function SampleWithANLS (p, flow_table, packet)
(8)
flow = flow_table·look-up (packet)
(9)
if flow is null then
(10)
flow= new Flow (packet)
(11)
else if p[flow·length] > rand_double (0, 1) then
(12)
flow·length+ = 1
(13)
flow_table·insert_or_update (flow)We calculate the sampling matrixS recursively for ANLS. In this method, the first packet is always sampled since p1,u=1 independent of u. We start by assigning all zero sampling probabilities as s:,0=0 and s1,1=1. The recursive equation for calculating the sampling matrix can be deduced as(6)si,ν∝si−1,ν−1pi−1,u+si−1,ν1−pi−1,u.Figure3 shows the ANLS sampling matrices for first 10 flow lengths and the Monday data sampled with them, respectively, similarly to Figure 2. First, we can see that when u is small, the S matrix looks like identity and as u gets larger, the sampling probability of large flows decreases. Secondly, compared to uniform sampling, the ANLS method has much higher sampling ratios than uniform sampling. However, operating with such high sampling ratios would require specialized hardware in real time.Figure 3
Visualization of the Monday slice with uniform and ANLS sampling methods.
(a)
(b)


## 3.3.1. Uniform Sampling Method

In uniform sampling, each packet is processed with a fixed probability ofπ, irrespective of the flow it belongs to. In this method, the sampling matrix entries si,ν are calculated directly through Binomial distribution with i trials and π success probability, i.e.,(3)si,ν=iνπν1−πi−ν,forν≤i,0,otherwise.Algorithm1 describes how the flow table is updated with uniform sampling upon the arrival of a new packet. The algorithm uniformly draws a random number in interval 0,1, and if it is less than π, the packet is processed; otherwise, the packet is discarded. For processing the packet, a look-up operation is performed on the flow table to find and update the flow that the packet belongs to. If no such flow is found, a new flow is created using the packet’s signature.Algorithm 1:Uniform packet sampling algorithm.
(1)
function Sample Uniformly Random π, flow_table, packet
(2)
if π > rand_double (0, 1) then
(3)
flow = flow_table·look-up (packet)
(4)
if flow is null then
(5)
flow= new Flow (packet)
(6)
else
(7)
flow·length+ = 1
(8)
flow_table·insert_or_update (flow)Figure2(a) shows the top 10×11 entries of the lower diagonal S matrices with different sampling probabilities. As the sampling probability π gets smaller, fewer packets from a flow gets observed, and the flow may even be missed when none of its packets are observed. The rightmost sampling matrix shows the case when π=1/64, where the matrix has a very high concentration of zero-length sampled flows.Figure 2
Sampling matrices for two different sampling schemes.
(a)
(b)Figure2(b) shows the original Monday data (the leftmost matrix) and its sampled versions under uniform sampling with the probabilities shown on the top sampling matrices. Here, we see that for π=1/64, the observed flow lengths are mostly less than 10, while the majority are not observed at all.


## 3.3.2. ANLS Sampling Method

The ANLS [10] will be used as the representative of the flow-based adaptive sampling methods. In ANLS, each packet is sampled according to the number of packets previously sampled from its corresponding flow. If a sampled flow has length x, the probability of its next packet to be sampled (px;u) is calculated as(4)fx;u=1+ux−1u,(5)px;u=1fx−1;u−fx;u.Here,fx;u is a monotonically increasing function of flow length x, parametrized with u, which makes px;u monotonically decreasing. u determines the tendency of sampling packets. As u gets smaller, more packets are sampled and estimating original flow lengths gets easier.The ANLS method is described in details in Algorithm2. Prior to the sampling, f and p are calculated in the SetupANLS function, according to equations (4) and (5). During sampling, for each incoming packet, a look-up operation is performed unconditionally. If the corresponding flow is found, it is updated with probability relative to its current observed size. Otherwise, a new flow is created, ensuring that every flow is observed with at least one packet.Algorithm 2:ANLS sampling algorithm.
(1)
function Setup ANLSu
(2)
f[0] = 0
(3)
for i∈1,I do
(4)
fi=1+ui−1/u
(5)
pi=1/fi−1−fi
(6)
return f, p
(7)
function SampleWithANLS (p, flow_table, packet)
(8)
flow = flow_table·look-up (packet)
(9)
if flow is null then
(10)
flow= new Flow (packet)
(11)
else if p[flow·length] > rand_double (0, 1) then
(12)
flow·length+ = 1
(13)
flow_table·insert_or_update (flow)We calculate the sampling matrixS recursively for ANLS. In this method, the first packet is always sampled since p1,u=1 independent of u. We start by assigning all zero sampling probabilities as s:,0=0 and s1,1=1. The recursive equation for calculating the sampling matrix can be deduced as(6)si,ν∝si−1,ν−1pi−1,u+si−1,ν1−pi−1,u.Figure3 shows the ANLS sampling matrices for first 10 flow lengths and the Monday data sampled with them, respectively, similarly to Figure 2. First, we can see that when u is small, the S matrix looks like identity and as u gets larger, the sampling probability of large flows decreases. Secondly, compared to uniform sampling, the ANLS method has much higher sampling ratios than uniform sampling. However, operating with such high sampling ratios would require specialized hardware in real time.Figure 3
Visualization of the Monday slice with uniform and ANLS sampling methods.
(a)
(b)


## 4. Methodology

Our methodology is based on the nonnegative factorization of the data tensor. Our model, which we call ThinNTF, introduces the sampling matrix as a constant factor to the original NTF with the Poisson–Gamma observation model. The rationale for using factorization for recovering true flow sizes is that the flow size distributions have daily periodic behavior, as we show in Section3. Inferring the factors, instead of the original matrix, requires less parameter estimation and results in smoother estimates compared to the standard maximum likelihood estimation described in [3]. In this section, we first describe the original tensor factorization model and we provide our novel version: the ThinNTF. At the end of the section, we present the full-Bayesian variational Bayes algorithm for the inference.


### 4.1. Nonnegative Tensor Factorization

NTF is the generalization of the 2-dimensional NMF model to multiple dimensions. In NTF, an N-dimensional tensor is approximated by the multiplication of lower dimensional factors. Unlike NMF, tensor factorization can be done in multiple ways. In this work, we are going to use the PARAFAC [22–24] factorization scheme. In PARAFAC, an I1×I2×⋯×IN tensor is approximated by In×R matrices for n∈1,N. Here, R is the number of components, i.e., the number of clusters in the data. Figure 4 shows the PARAFAC factorization of our FLD tensor X, into 3 factors: an I×R factor F for representing the flow length clusters, a J×R factor H for representing hourly behavior, and a K×R factor D for representing the daily behavior of the data. Every single entry of the X tensor is approximated by(7)xi,j,k≈x^i,j,k=∑rfi,rhj,rdk,r.Figure 4
PARAFAC factorization.Bro [25] explains that the PARAFAC factorization is unique under certain circumstances, where uniqueness is defined as begin unable to rotate the factorization without loss of fit. NMF and NTF are statistical models that impose nonnegativity constraint without uniqueness property. The uniqueness may be important if individual factors are of special interest. In our case, we are concerned with the estimation of the original data tensor X from sampled tensor Y, but not the individual factors for any interpretation. Our problem is more close to a missing value imputation problem; hence, uniqueness is not a requirement.


### 4.2. Thin Nonnegative Tensor Factorization

ThinNTF is basically an NTF with an additional constant factor, which in our case is the sampling matrixS. Figure 5 shows the graphical models of the NTF and the ThinNTF models for factorizing original and sampled flow length observations. In the graphical models, the shaded nodes are the observed entities and the unshaded ones are the latent entities.Figure 5
Graphical models representing the dependency structure of NTF and ThinNTF models in PARAFAC scheme.
(a)
(b)In Section3.3, we have described the sampling process as a matrix multiplication operation with a sampling matrix S. In ThinNTF, this sampling matrix operates on the original tensor X and creates a thinned version of it, which we call Y, by down-sampling its entries according to a sampling scheme, as shown in Figure 6. The entries of Y tensor yν,j,k presents the number of flows of sampled-length ν, at hour j at day k. The ⊗1 operation denotes the 1-mode product of matrix ST and tensor X, which corresponds to the set of matrix multiplications Y:,:,k=STX:,:,k for k∈1,K.Figure 6
ThinNTF model.In this scheme, one can immediately suspect thatx can be estimated by ST−1⊗1Y. However, this solution is not feasible for several reasons. First, the S matrix is not square and hence not invertible. Instead, its pseudoinverse can be calculated, but this does not impose nonnegativity. Moreover, the top slice of the Y tensor, which stores the number of flows with zero-sampled size, is never observed and hence must be estimated. Therefore, we need a solid statistical model and an inference method to estimate X under this model.In ThinNTF, we observe theY tensor, but try to factorize the X tensor, which is latent (Figure 5(b)). In the end, the factors of X are going to provide us an approximation X^ which solves the original flow length distribution reconstruction problem. We mathematically express this approximation as(8)yν,j,k≈y^ν,j,k=∑i,rsν,ifi,rhj,rdk,r,where F,H, and D are described in exactly the same way in the original NTF case. In Subsections 3.3.1 and 3.3.2, we described two different S matrices for two different schemes. ThinNTF model can be employed with any sampling method as long as it is described with a sampling matrix.


### 4.3. Generative Model

Taking the Bayesian approach, we first provide a generative model for the ThinNTF and then describe how we can estimate the posterior probabilities of model parameters (in this case, the factor matrices) conditioned on the sampled flow length observationsY and the sampling matrix S using the well-known Bayes rule. Table 2 contains all tensors and matrices used in the model together with their index sets.Table 2
Tensors in the model and their corresponding index sets.
Tensor
Index set
Description
X
i
,
j
,
k
Original flow length tensor
Y
ν
,
j
,
k
Sampled flow length tensor
M
ν
,
j
,
k
Mask tensor
W
ν
,
i
,
j
,
k
,
r
Latent variable tensor
F
i
,
r
Flow length factor
H
j
,
r
Hour of day factor
D
k
,
r
Day of week factor
S
i
,
ν
Sampling matrix
A
F
,
B
F
i
,
r
Gamma priors forF
A
H
,
B
H
j
,
r
Gamma priors forH
A
D
,
B
D
k
,
r
Gamma priors forDThe original and latent data tensorX and the sampled and observed data tensor Y have nonnegative integer entries. The natural probability distribution for this type of count data is the Poisson distribution. We assume that each entry of a latent 5-dimensional tensor W is drawn from a Poisson distribution whose parameters are functions of sampling matrix S and factors F,H, and D, such as(9)wν,i,j,k,r∼POwν,i,j,k,r;sν,ifi,rhj,rdk,r,where the Poisson distribution is defined as(10)POw;λ=expwlogλ−λ−logΓw+1.We choose the prior distributions for the factor entries as the Gamma distribution since it is the conjugate prior of Poisson distribution [26]. For each entry of factor F, we write(11)fi,r∼Gfi,r;ai,rf,bi,rfai,rf,where the Gamma distribution is described as(12)Gf;κ,Θ=expκ−1logf−fΘ−κlogΘ−logΓκ,with shape parameter κ and scale parameter θ. In our generative model, the parameters for Gamma distributions are κ=ai,rf and Θ=bi,rf/ai,rf respectively. This means that the mean of fi,r is κΘ=bi,rf, which is independent of ai,rf. The variance of fi,r becomes κΘ2=bi,rf2/ai,rf, which means that as ai,rf gets smaller, the factors gets sparser.In order to avoid repetition, we are going to omit the equations regarding the factorsH and D throughout the paper. These factors behave exactly like factor F, and it is easy to derive equations related to these factors once their corresponding equation for F is given.Finally, we generateX and Y tensor from W. Each entry wν,i,j,k,r of W can be interpreted as the number of original flows of length i, generated on hour j, day k, by cluster r, and observed as length ν. By summing W over dimensions cluster r and original lengths i, we get the sampled observations tensor Y. Similarly, by summing W over dimensions cluster r and sampled lengths ν, we get the original flow length tensor X. The whole generative process is summarized in Algorithm 3. The set of all indexes and tensors in the model are summarized in Tables 1 and 2, respectively.Algorithm 3:TNTF generative model.
(1)
function RandInit (S,AF,BF,AH,BH,AD,BD)
//Sample factorF from Gamma (AF,BF)
(2)
for all i∈1,I,r∈1,R do
(3)
fi,r∼Gfi,r;ai,rf,bi,rf/ai,rf
//Sample factorH from Gamma (AH,BH)
(4)
for all k∈1,K,r∈1,R do
(5)
hj,r∼Ghj,r;aj,rh,bj,rh/aj,rh
//Sample factorD from Gamma (AD,BD)
(6)
for all j∈1,J,r∈1,R do
(7)
dk,r∼Gdk,r;ak,rd,bk,rd/ak,rd
//Sample latent tensorW from Poisson distributions
(8)
for all ν∈1,I+1,i∈1,I,j∈1,J,k∈1,K,r∈1,R do
(9)
wν,i,j,k,r∼POwν,i,j,k,r;sν,ifi,rhj,rdk,r
(10)
return {F, H, D, W}
(11)
function GenerateData (S,AF,BF,AH,BH,AD,BD)
//Randomly initialize factors and latent tensor
(12)
F,H,D,W ⟵ RandInit (S,AF,BF,AH,BH,AD,BD)
//Generate original tensorx
(13)
for all i∈1,I,j∈1,J,k∈1,K do
(14)
xi,j,k=∑ν,rWν,i,j,k,r
//Generate sampled tensorY
(15)
for all ν∈1,I+1,j∈1,J,k∈1,K do
(16)
yν,j,k=∑i,rwν,i,j,k,r
(17)
return {F, H, D, W, X, Y}


### 4.4. Variational Bayes

After defining the generative model, we can inter the factorsF,H, and D of a sampled flow length observation tensor Y. In the original NMF paper, Lee, and Seung [17] provide fixed-point update equations for inferring the factors. Bro [25] gives similar fixed-point equations for updating the factors in PARAFAC factorization. Cemgil [19] shows that these updates correspond to the Kullback–Leibler minimization between the original matrix (or tensor X) and the approximated one (X^) and also provides a full Bayesian variational algorithm for the matrix factorization. Ermis et al. [15] provide a similar variational algorithm for the Gamma–Poisson tensor factorization.We start our Bayesian inference by calculating the posterior distributions over the factorsF,H, and D conditioned on observed tensor Y. For notational clarity, we introduce θ=AF,BF,AH,BH,AD,BD as the list of model hyperparameters. The log-likelihood observing Y under the model parameters θ is written as(13)logpY∣θ,S=log∫F,H,DdFdHdD∑WpY,W,F,H,D∣θ,S.This log-likelihood is intractable due to the integration over the latent factors, but it is lower bounded as(14)logpY∣θ,S≤Lθ=logpY,W,F,H,Dθ,SqW,F,H,D+HqW,F,H,D,where q is an auxiliary joint distribution of latent factors. This bound is tight when qW,F,H,D=pW,F,H,DY,θ,S. However, this is also intractable to calculate. Instead, we use a variational approximation [27] for q such that(15)qW∝explogpY,W,F,H,D∣θqF,H,D,qF∝explogpY,W,F,H,D∣θqW,H,D,qH∝explogpY,W,F,H,D∣θqW,F,D,qD∝explogpY,W,F,H,D∣θqW,F,H,where we iteratively update the posterior distribution of each factor by calculating the expectation of the logarithm of the full joint likelihood pY,W,F,H,D under the posteriors of all other latent factors.


### 4.5. Update Equations

Here, we provide the update equations forqW and qF. The updates of qH and qD can be easily deduced from the update equations of qF. The full joint likelihood whose expectation will be calculated at each step is(16)Jθ=logpY,W,F,H,Dθ=logpY∣W+logpW∣F,H,D+logpFθ+logpH∣θ+logpD∣θ,where YW is a degenerate distribution (δ) to make sure the summation of ∑i,rW equals Y. By inserting the necessary Poisson and Gamma distributions given in the generative model in the equation (16), we get the following expression:(17)Jθ=∑ν,j,kmν,j,klogδyν,j,k−∑i,rwν,i,j,k,r+∑i,r∑ν,j,kmν,j,kwν,i,j,k,rlogsν,ifi,rhj,rdk,r−sν,ifi,rhj,rdk,r−logΓwν,i,j,k,r+1+∑i,rai,rf−1logfi,r−fi,rai,rfbi,rf−ai,rflogbi,rfai,rf−logΓai,rf+∑j,raj,rh−1loghj,r−hj,raj,rhbj,rh−aj,rhlogbj,rhaj,rh−logΓaj,rh+∑k,rak,rd−1logdk,r−dk,rak,rdbk,rd−ak,rdlogbk,rdak,rd−logΓak,rd.


#### 4.5.1. Update Rule forqW

Considering the terms in the log-likelihood expression in equation (17), that only includes wν,i,j,k,r, we find that(18)qwv,:,j,k,:∝expmν,j,klogδyν,j,k−∑i,rwν,i,j,k,r+∑i,rmν,j,kwν,i,j,k,rlogsν,ifi,rhj,rdk,r−logΓwν,i,j,k,r+1∝multinomialwν,j,k,:,xi,j,k,pν,i,j,k,rmν,j,k,where, wν,j,k,:,: becomes multinomial distributed. The expectation of W is calculated as(19)pν,i,j,k,r=expsν,i+logfi,r+loghj,r+logdk,r∑i,rexpsν,i+logfi,r+loghj,r+logdk,r,wν,i,j,k,r=yν,j,kpν,i,j,k,r.


#### 4.5.2. Update Rule forqF

Similarly, considering the terms in log-likelihood equation (17) that only includes fi,r, we find that(20)qfi,r∝∑ν,j,kmν,j,kwν,i,j,k,r+ai,rf−1logfi,r−∑ν,j,kmν,j,ksν,ihj,rdk,r+ai,rfbi,rffi,r∝Gammafi,r;αi,rf,βi,rf,where fi,r becomes Gamma distributed with shape and scale parameters(21)αi,rf=ai,rf+∑ν,j,kmν,j,kwν,i,j,k,r,(22)βi,rf=ai,rfbi,rf+∑ν,j,kmν,j,ksν,ihj,rdk,r−1.We calculate the expectation offi,r and the logarithm of fi,r as(23)fi,r=αi,rfβi,rf,(24)logfi,r=Ψαi,rf+logβi,rf.The variational Bayes algorithm that uses the above equations is given in Algorithm4. The calculation of the lower bound is given in Appendix. The exact derivations of all equations can be found in [28].Algorithm 4:Variational Bayes algorithm.
(1)
function ThinNTF_VB (Y, S, AF,BF,AH,BH,AD,BD)
//Randomly initialize factors and latent tensor
(2)
F,H,D,W ⟵ RandInit (S,AF,BF,AH,BH,AD,BD)
(3)
repeat
(4)
Calculateαi,rf,βi,rf and fi,r as in equations (21)–(23).
(5)
Calculateαj,rh,βk,rh and hj,r similarly.
(6)
Calculateαk,rd,βk,rd and dk,r similarly.
(7)
Calculatelogfi,r as in equation (24).
(8)
Calculatelogfi,r similarly.
(9)
Calculatelogfi,r similarly.
(10)
Calculatewν,i,j,k,r as in equation (19).
(11)
Calculate lower bound
(12)
until Max iterations are reached or lower bound converged
(13)
return F, H, D, X


### 4.6. Computational Complexity

The nonnegative tensor factorization is an NP-hard problem [29]. The variational Bayes algorithm we introduced in Algorithm 4 is an iterative solution that converges to a local maximum solution. The complexity of each iteration is determined by the leading term, which is the equation (19). In general, calculating a ThinNTF model with R components for a κ-dimensional tensor with all dimensions of length N has OκNκ+1R complexity for a single iteration.


## 4.1. Nonnegative Tensor Factorization

NTF is the generalization of the 2-dimensional NMF model to multiple dimensions. In NTF, an N-dimensional tensor is approximated by the multiplication of lower dimensional factors. Unlike NMF, tensor factorization can be done in multiple ways. In this work, we are going to use the PARAFAC [22–24] factorization scheme. In PARAFAC, an I1×I2×⋯×IN tensor is approximated by In×R matrices for n∈1,N. Here, R is the number of components, i.e., the number of clusters in the data. Figure 4 shows the PARAFAC factorization of our FLD tensor X, into 3 factors: an I×R factor F for representing the flow length clusters, a J×R factor H for representing hourly behavior, and a K×R factor D for representing the daily behavior of the data. Every single entry of the X tensor is approximated by(7)xi,j,k≈x^i,j,k=∑rfi,rhj,rdk,r.Figure 4
PARAFAC factorization.Bro [25] explains that the PARAFAC factorization is unique under certain circumstances, where uniqueness is defined as begin unable to rotate the factorization without loss of fit. NMF and NTF are statistical models that impose nonnegativity constraint without uniqueness property. The uniqueness may be important if individual factors are of special interest. In our case, we are concerned with the estimation of the original data tensor X from sampled tensor Y, but not the individual factors for any interpretation. Our problem is more close to a missing value imputation problem; hence, uniqueness is not a requirement.


## 4.2. Thin Nonnegative Tensor Factorization

ThinNTF is basically an NTF with an additional constant factor, which in our case is the sampling matrixS. Figure 5 shows the graphical models of the NTF and the ThinNTF models for factorizing original and sampled flow length observations. In the graphical models, the shaded nodes are the observed entities and the unshaded ones are the latent entities.Figure 5
Graphical models representing the dependency structure of NTF and ThinNTF models in PARAFAC scheme.
(a)
(b)In Section3.3, we have described the sampling process as a matrix multiplication operation with a sampling matrix S. In ThinNTF, this sampling matrix operates on the original tensor X and creates a thinned version of it, which we call Y, by down-sampling its entries according to a sampling scheme, as shown in Figure 6. The entries of Y tensor yν,j,k presents the number of flows of sampled-length ν, at hour j at day k. The ⊗1 operation denotes the 1-mode product of matrix ST and tensor X, which corresponds to the set of matrix multiplications Y:,:,k=STX:,:,k for k∈1,K.Figure 6
ThinNTF model.In this scheme, one can immediately suspect thatx can be estimated by ST−1⊗1Y. However, this solution is not feasible for several reasons. First, the S matrix is not square and hence not invertible. Instead, its pseudoinverse can be calculated, but this does not impose nonnegativity. Moreover, the top slice of the Y tensor, which stores the number of flows with zero-sampled size, is never observed and hence must be estimated. Therefore, we need a solid statistical model and an inference method to estimate X under this model.In ThinNTF, we observe theY tensor, but try to factorize the X tensor, which is latent (Figure 5(b)). In the end, the factors of X are going to provide us an approximation X^ which solves the original flow length distribution reconstruction problem. We mathematically express this approximation as(8)yν,j,k≈y^ν,j,k=∑i,rsν,ifi,rhj,rdk,r,where F,H, and D are described in exactly the same way in the original NTF case. In Subsections 3.3.1 and 3.3.2, we described two different S matrices for two different schemes. ThinNTF model can be employed with any sampling method as long as it is described with a sampling matrix.


## 4.3. Generative Model

Taking the Bayesian approach, we first provide a generative model for the ThinNTF and then describe how we can estimate the posterior probabilities of model parameters (in this case, the factor matrices) conditioned on the sampled flow length observationsY and the sampling matrix S using the well-known Bayes rule. Table 2 contains all tensors and matrices used in the model together with their index sets.Table 2
Tensors in the model and their corresponding index sets.
Tensor
Index set
Description
X
i
,
j
,
k
Original flow length tensor
Y
ν
,
j
,
k
Sampled flow length tensor
M
ν
,
j
,
k
Mask tensor
W
ν
,
i
,
j
,
k
,
r
Latent variable tensor
F
i
,
r
Flow length factor
H
j
,
r
Hour of day factor
D
k
,
r
Day of week factor
S
i
,
ν
Sampling matrix
A
F
,
B
F
i
,
r
Gamma priors forF
A
H
,
B
H
j
,
r
Gamma priors forH
A
D
,
B
D
k
,
r
Gamma priors forDThe original and latent data tensorX and the sampled and observed data tensor Y have nonnegative integer entries. The natural probability distribution for this type of count data is the Poisson distribution. We assume that each entry of a latent 5-dimensional tensor W is drawn from a Poisson distribution whose parameters are functions of sampling matrix S and factors F,H, and D, such as(9)wν,i,j,k,r∼POwν,i,j,k,r;sν,ifi,rhj,rdk,r,where the Poisson distribution is defined as(10)POw;λ=expwlogλ−λ−logΓw+1.We choose the prior distributions for the factor entries as the Gamma distribution since it is the conjugate prior of Poisson distribution [26]. For each entry of factor F, we write(11)fi,r∼Gfi,r;ai,rf,bi,rfai,rf,where the Gamma distribution is described as(12)Gf;κ,Θ=expκ−1logf−fΘ−κlogΘ−logΓκ,with shape parameter κ and scale parameter θ. In our generative model, the parameters for Gamma distributions are κ=ai,rf and Θ=bi,rf/ai,rf respectively. This means that the mean of fi,r is κΘ=bi,rf, which is independent of ai,rf. The variance of fi,r becomes κΘ2=bi,rf2/ai,rf, which means that as ai,rf gets smaller, the factors gets sparser.In order to avoid repetition, we are going to omit the equations regarding the factorsH and D throughout the paper. These factors behave exactly like factor F, and it is easy to derive equations related to these factors once their corresponding equation for F is given.Finally, we generateX and Y tensor from W. Each entry wν,i,j,k,r of W can be interpreted as the number of original flows of length i, generated on hour j, day k, by cluster r, and observed as length ν. By summing W over dimensions cluster r and original lengths i, we get the sampled observations tensor Y. Similarly, by summing W over dimensions cluster r and sampled lengths ν, we get the original flow length tensor X. The whole generative process is summarized in Algorithm 3. The set of all indexes and tensors in the model are summarized in Tables 1 and 2, respectively.Algorithm 3:TNTF generative model.
(1)
function RandInit (S,AF,BF,AH,BH,AD,BD)
//Sample factorF from Gamma (AF,BF)
(2)
for all i∈1,I,r∈1,R do
(3)
fi,r∼Gfi,r;ai,rf,bi,rf/ai,rf
//Sample factorH from Gamma (AH,BH)
(4)
for all k∈1,K,r∈1,R do
(5)
hj,r∼Ghj,r;aj,rh,bj,rh/aj,rh
//Sample factorD from Gamma (AD,BD)
(6)
for all j∈1,J,r∈1,R do
(7)
dk,r∼Gdk,r;ak,rd,bk,rd/ak,rd
//Sample latent tensorW from Poisson distributions
(8)
for all ν∈1,I+1,i∈1,I,j∈1,J,k∈1,K,r∈1,R do
(9)
wν,i,j,k,r∼POwν,i,j,k,r;sν,ifi,rhj,rdk,r
(10)
return {F, H, D, W}
(11)
function GenerateData (S,AF,BF,AH,BH,AD,BD)
//Randomly initialize factors and latent tensor
(12)
F,H,D,W ⟵ RandInit (S,AF,BF,AH,BH,AD,BD)
//Generate original tensorx
(13)
for all i∈1,I,j∈1,J,k∈1,K do
(14)
xi,j,k=∑ν,rWν,i,j,k,r
//Generate sampled tensorY
(15)
for all ν∈1,I+1,j∈1,J,k∈1,K do
(16)
yν,j,k=∑i,rwν,i,j,k,r
(17)
return {F, H, D, W, X, Y}


## 4.4. Variational Bayes

After defining the generative model, we can inter the factorsF,H, and D of a sampled flow length observation tensor Y. In the original NMF paper, Lee, and Seung [17] provide fixed-point update equations for inferring the factors. Bro [25] gives similar fixed-point equations for updating the factors in PARAFAC factorization. Cemgil [19] shows that these updates correspond to the Kullback–Leibler minimization between the original matrix (or tensor X) and the approximated one (X^) and also provides a full Bayesian variational algorithm for the matrix factorization. Ermis et al. [15] provide a similar variational algorithm for the Gamma–Poisson tensor factorization.We start our Bayesian inference by calculating the posterior distributions over the factorsF,H, and D conditioned on observed tensor Y. For notational clarity, we introduce θ=AF,BF,AH,BH,AD,BD as the list of model hyperparameters. The log-likelihood observing Y under the model parameters θ is written as(13)logpY∣θ,S=log∫F,H,DdFdHdD∑WpY,W,F,H,D∣θ,S.This log-likelihood is intractable due to the integration over the latent factors, but it is lower bounded as(14)logpY∣θ,S≤Lθ=logpY,W,F,H,Dθ,SqW,F,H,D+HqW,F,H,D,where q is an auxiliary joint distribution of latent factors. This bound is tight when qW,F,H,D=pW,F,H,DY,θ,S. However, this is also intractable to calculate. Instead, we use a variational approximation [27] for q such that(15)qW∝explogpY,W,F,H,D∣θqF,H,D,qF∝explogpY,W,F,H,D∣θqW,H,D,qH∝explogpY,W,F,H,D∣θqW,F,D,qD∝explogpY,W,F,H,D∣θqW,F,H,where we iteratively update the posterior distribution of each factor by calculating the expectation of the logarithm of the full joint likelihood pY,W,F,H,D under the posteriors of all other latent factors.


## 4.5. Update Equations

Here, we provide the update equations forqW and qF. The updates of qH and qD can be easily deduced from the update equations of qF. The full joint likelihood whose expectation will be calculated at each step is(16)Jθ=logpY,W,F,H,Dθ=logpY∣W+logpW∣F,H,D+logpFθ+logpH∣θ+logpD∣θ,where YW is a degenerate distribution (δ) to make sure the summation of ∑i,rW equals Y. By inserting the necessary Poisson and Gamma distributions given in the generative model in the equation (16), we get the following expression:(17)Jθ=∑ν,j,kmν,j,klogδyν,j,k−∑i,rwν,i,j,k,r+∑i,r∑ν,j,kmν,j,kwν,i,j,k,rlogsν,ifi,rhj,rdk,r−sν,ifi,rhj,rdk,r−logΓwν,i,j,k,r+1+∑i,rai,rf−1logfi,r−fi,rai,rfbi,rf−ai,rflogbi,rfai,rf−logΓai,rf+∑j,raj,rh−1loghj,r−hj,raj,rhbj,rh−aj,rhlogbj,rhaj,rh−logΓaj,rh+∑k,rak,rd−1logdk,r−dk,rak,rdbk,rd−ak,rdlogbk,rdak,rd−logΓak,rd.


### 4.5.1. Update Rule forqW

Considering the terms in the log-likelihood expression in equation (17), that only includes wν,i,j,k,r, we find that(18)qwv,:,j,k,:∝expmν,j,klogδyν,j,k−∑i,rwν,i,j,k,r+∑i,rmν,j,kwν,i,j,k,rlogsν,ifi,rhj,rdk,r−logΓwν,i,j,k,r+1∝multinomialwν,j,k,:,xi,j,k,pν,i,j,k,rmν,j,k,where, wν,j,k,:,: becomes multinomial distributed. The expectation of W is calculated as(19)pν,i,j,k,r=expsν,i+logfi,r+loghj,r+logdk,r∑i,rexpsν,i+logfi,r+loghj,r+logdk,r,wν,i,j,k,r=yν,j,kpν,i,j,k,r.


### 4.5.2. Update Rule forqF

Similarly, considering the terms in log-likelihood equation (17) that only includes fi,r, we find that(20)qfi,r∝∑ν,j,kmν,j,kwν,i,j,k,r+ai,rf−1logfi,r−∑ν,j,kmν,j,ksν,ihj,rdk,r+ai,rfbi,rffi,r∝Gammafi,r;αi,rf,βi,rf,where fi,r becomes Gamma distributed with shape and scale parameters(21)αi,rf=ai,rf+∑ν,j,kmν,j,kwν,i,j,k,r,(22)βi,rf=ai,rfbi,rf+∑ν,j,kmν,j,ksν,ihj,rdk,r−1.We calculate the expectation offi,r and the logarithm of fi,r as(23)fi,r=αi,rfβi,rf,(24)logfi,r=Ψαi,rf+logβi,rf.The variational Bayes algorithm that uses the above equations is given in Algorithm4. The calculation of the lower bound is given in Appendix. The exact derivations of all equations can be found in [28].Algorithm 4:Variational Bayes algorithm.
(1)
function ThinNTF_VB (Y, S, AF,BF,AH,BH,AD,BD)
//Randomly initialize factors and latent tensor
(2)
F,H,D,W ⟵ RandInit (S,AF,BF,AH,BH,AD,BD)
(3)
repeat
(4)
Calculateαi,rf,βi,rf and fi,r as in equations (21)–(23).
(5)
Calculateαj,rh,βk,rh and hj,r similarly.
(6)
Calculateαk,rd,βk,rd and dk,r similarly.
(7)
Calculatelogfi,r as in equation (24).
(8)
Calculatelogfi,r similarly.
(9)
Calculatelogfi,r similarly.
(10)
Calculatewν,i,j,k,r as in equation (19).
(11)
Calculate lower bound
(12)
until Max iterations are reached or lower bound converged
(13)
return F, H, D, X


## 4.5.1. Update Rule forqW

Considering the terms in the log-likelihood expression in equation (17), that only includes wν,i,j,k,r, we find that(18)qwv,:,j,k,:∝expmν,j,klogδyν,j,k−∑i,rwν,i,j,k,r+∑i,rmν,j,kwν,i,j,k,rlogsν,ifi,rhj,rdk,r−logΓwν,i,j,k,r+1∝multinomialwν,j,k,:,xi,j,k,pν,i,j,k,rmν,j,k,where, wν,j,k,:,: becomes multinomial distributed. The expectation of W is calculated as(19)pν,i,j,k,r=expsν,i+logfi,r+loghj,r+logdk,r∑i,rexpsν,i+logfi,r+loghj,r+logdk,r,wν,i,j,k,r=yν,j,kpν,i,j,k,r.


## 4.5.2. Update Rule forqF

Similarly, considering the terms in log-likelihood equation (17) that only includes fi,r, we find that(20)qfi,r∝∑ν,j,kmν,j,kwν,i,j,k,r+ai,rf−1logfi,r−∑ν,j,kmν,j,ksν,ihj,rdk,r+ai,rfbi,rffi,r∝Gammafi,r;αi,rf,βi,rf,where fi,r becomes Gamma distributed with shape and scale parameters(21)αi,rf=ai,rf+∑ν,j,kmν,j,kwν,i,j,k,r,(22)βi,rf=ai,rfbi,rf+∑ν,j,kmν,j,ksν,ihj,rdk,r−1.We calculate the expectation offi,r and the logarithm of fi,r as(23)fi,r=αi,rfβi,rf,(24)logfi,r=Ψαi,rf+logβi,rf.The variational Bayes algorithm that uses the above equations is given in Algorithm4. The calculation of the lower bound is given in Appendix. The exact derivations of all equations can be found in [28].Algorithm 4:Variational Bayes algorithm.
(1)
function ThinNTF_VB (Y, S, AF,BF,AH,BH,AD,BD)
//Randomly initialize factors and latent tensor
(2)
F,H,D,W ⟵ RandInit (S,AF,BF,AH,BH,AD,BD)
(3)
repeat
(4)
Calculateαi,rf,βi,rf and fi,r as in equations (21)–(23).
(5)
Calculateαj,rh,βk,rh and hj,r similarly.
(6)
Calculateαk,rd,βk,rd and dk,r similarly.
(7)
Calculatelogfi,r as in equation (24).
(8)
Calculatelogfi,r similarly.
(9)
Calculatelogfi,r similarly.
(10)
Calculatewν,i,j,k,r as in equation (19).
(11)
Calculate lower bound
(12)
until Max iterations are reached or lower bound converged
(13)
return F, H, D, X


## 4.6. Computational Complexity

The nonnegative tensor factorization is an NP-hard problem [29]. The variational Bayes algorithm we introduced in Algorithm 4 is an iterative solution that converges to a local maximum solution. The complexity of each iteration is determined by the leading term, which is the equation (19). In general, calculating a ThinNTF model with R components for a κ-dimensional tensor with all dimensions of length N has OκNκ+1R complexity for a single iteration.


## 5. Data Collection

Applying a tensor model to the flow length estimation problem requires high-volume data collected over a long period, to capture the timely behavior of the network. The already available online data sets do not fulfill this requirement. Therefore, we collected our own real-world data from a mobile network service provider in Turkey [30]. The architecture of our system and the description of the data we collected are presented as follows.


### 5.1. System Architecture

The system architecture of a mobile operator’s general packet radio service (GPRS) network infrastructure, including radio access and core network elements, is illustrated in Figure7. IP traffic generated or received by mobile devices between mobile station (MS) and packet data network PDN, e.g., IP Multimedia Subsystem (IMS), is tunneled in the core network of mobile operators through serving GPRS support node (SGSN) and gateway GPRS support node (GGSN) via the user data part of the GPRS tunneling protocol (GTP) [31]. The GTP message exchanges include information such as the size of the traffic, IP session start and end time, and device and user identifiers.Figure 7
Placement of our monitoring server inside the premises of the mobile operator running a commercial cellular network.The Gn interface (Gn is an interface between SGSN and GGSN where GTP is the main protocol for network packets flowing through) carries user packets to be transferred between the mobile users and the Internet together with control packets necessary for the universal mobile telecommunications service (UMTS) core network [32]. All packets in this channel are carried by the GTP, which is an IP-based protocol for carrying GPRS data within UMTS networks, used for data encapsulation in order to keep the core network independent of the protocols that are being used between MS and the packet-switched network.The Gn interface carries mainly two types of GTP message structures or packets: GTP-C and GTP-U. GTP-C is used for signaling between SGSN and GGSN in core network which carries packet data protocol (PDP) context messages such as activating and deactivating user session, configuring service parameters or updating the session. GTP-U is used for transmitting user data between the radio access network and core network. In our experiments, we filtered out GTP-C packets (since they are not considered to be part of a flow due to flow definition), which makes 10% of the total Gn traffic. Therefore, the sampling is applied to GTP-U packets only. GTP is carried mainly over UDP.


### 5.2. Data Extraction Process

The mobile operator network consists of several districts with more than 10 regional core areas throughout Turkey. The average total traffic in all regional areas consists of approximately over 15 billion packets in the uplink direction and over 20 billion packets in the downlink direction daily. This corresponds to approximately 80 terabytes of total data flowing in uplink and downlink daily inside the mobile operator’s core network as a whole. In this work, the Gn interface which connects the SGSN and GGSN nodes are mirrored, and the network traffic is transferred into a FLD server located at mobile operator’s technology center in the core network. A speed of 200 Mbit/sec at peak hours can be observed through one of the mirrored interfaces in the core network.We monitored the network traffic in one of the servers of a mobile operator continuously for 10 days. We developed a packet extraction tool inside the monitoring server shown in Figure7, which parses each GTP-U packet and stores the packet signature together with packet length and arrival time, discarding their payload. The stored network data is processed offline for extracting the true flow lengths.After the data collection and flow extraction, the total number of packets collected is found to be4×1011, which makes up around 2.5×1010 flows. Figure 8 shows the cumulative flow length distribution of the data. We see that most of the flows have less than 5 packets, and 99.9% of them have less than 100 packets.Figure 8
Cumulative flow lengths in the real-world data.


## 5.1. System Architecture

The system architecture of a mobile operator’s general packet radio service (GPRS) network infrastructure, including radio access and core network elements, is illustrated in Figure7. IP traffic generated or received by mobile devices between mobile station (MS) and packet data network PDN, e.g., IP Multimedia Subsystem (IMS), is tunneled in the core network of mobile operators through serving GPRS support node (SGSN) and gateway GPRS support node (GGSN) via the user data part of the GPRS tunneling protocol (GTP) [31]. The GTP message exchanges include information such as the size of the traffic, IP session start and end time, and device and user identifiers.Figure 7
Placement of our monitoring server inside the premises of the mobile operator running a commercial cellular network.The Gn interface (Gn is an interface between SGSN and GGSN where GTP is the main protocol for network packets flowing through) carries user packets to be transferred between the mobile users and the Internet together with control packets necessary for the universal mobile telecommunications service (UMTS) core network [32]. All packets in this channel are carried by the GTP, which is an IP-based protocol for carrying GPRS data within UMTS networks, used for data encapsulation in order to keep the core network independent of the protocols that are being used between MS and the packet-switched network.The Gn interface carries mainly two types of GTP message structures or packets: GTP-C and GTP-U. GTP-C is used for signaling between SGSN and GGSN in core network which carries packet data protocol (PDP) context messages such as activating and deactivating user session, configuring service parameters or updating the session. GTP-U is used for transmitting user data between the radio access network and core network. In our experiments, we filtered out GTP-C packets (since they are not considered to be part of a flow due to flow definition), which makes 10% of the total Gn traffic. Therefore, the sampling is applied to GTP-U packets only. GTP is carried mainly over UDP.


## 5.2. Data Extraction Process

The mobile operator network consists of several districts with more than 10 regional core areas throughout Turkey. The average total traffic in all regional areas consists of approximately over 15 billion packets in the uplink direction and over 20 billion packets in the downlink direction daily. This corresponds to approximately 80 terabytes of total data flowing in uplink and downlink daily inside the mobile operator’s core network as a whole. In this work, the Gn interface which connects the SGSN and GGSN nodes are mirrored, and the network traffic is transferred into a FLD server located at mobile operator’s technology center in the core network. A speed of 200 Mbit/sec at peak hours can be observed through one of the mirrored interfaces in the core network.We monitored the network traffic in one of the servers of a mobile operator continuously for 10 days. We developed a packet extraction tool inside the monitoring server shown in Figure7, which parses each GTP-U packet and stores the packet signature together with packet length and arrival time, discarding their payload. The stored network data is processed offline for extracting the true flow lengths.After the data collection and flow extraction, the total number of packets collected is found to be4×1011, which makes up around 2.5×1010 flows. Figure 8 shows the cumulative flow length distribution of the data. We see that most of the flows have less than 5 packets, and 99.9% of them have less than 100 packets.Figure 8
Cumulative flow lengths in the real-world data.


## 6. Experiments and Results

We designed two sets of experiments in order to verify our model: synthetic and real-world experiments. In each set, we sampled the original data with both uniform and ANLS models with different sampling parameters. Then, we tried to recover the original tensor with ThinNTF models. The ThinNTF model takes a single parameterR which is the number of components in each factor. Additionally, we also represented data as I×JK matrix by unfolding the X tensor in the first dimension as described in [28] and applied the 2-dimensional version of ThinNTF, which we simply call thin nonnegative matrix factorization (ThinNMF). For Uniform sampling, we used the maximum likelihood estimation (MLE) defined in [3] as the baseline. For ANLS, we used both MLE and its own unbiased estimator of the model as baselines.Both ThinNMF and ThinNTF models explain the data as a linear combination ofR flow length distributions, stored in the columns of F matrix. In the ThinNMF model, we have JK coefficient sets for this combination. On the other hand, in ThinNTF, we have J coefficients for hour-of-day and K coefficients for day-of-week. The Cartesian products of these coefficient sets make a total of JK coefficients and create a dependency between the hour and day attributes. Therefore, we expect that ThinNTF captures the periodicity and give better estimates.During the experiments, we always run the stochastic algorithms, i.e., ThinNMF, ThinNTF, and MLE, for 10 times and keep the parameters of the model with the highest lower bound value. Then, we reported the success of our algorithm with the weighted mean relative distance (WMRD) metric as this was used in all previous flow size estimation works. The WMRD is a metric which gives more weights to the relative differences that occur with larger frequency. It is formulated as(25)wmrdx,X^=∑ixi−x^i∑ixi+x^i/2,where x is an original flow size distribution measured at the end of the hour and x^ is its estimated version. For the whole tensor X, we calculate the average WMRD value.Additionally, we report the Kullback–Leibler divergence between the original and the estimated tensors, since this is the metric minimized during the variational Bayes algorithm. The KL divergence between two distributionsx and X^ is calculated as(26)KLx,x^=∑ixilogx^ixi.


### 6.1. Experiments on Synthetic Data

We prepared our synthetic experiments to test the validity of our models. In this experiment set, we used the generative model of the ThinNTF model as described in Algorithm3 to generate a small network with maximum flow size I=10,J=24 and K=7. The original synthetic flow length distribution X is generated by a generative model with 3 components, where each component is a column in the F factor. We selected these 3 components as exponential, inverted exponential, and uniform random distributions. Therefore, in experiments, we used ThinNMF-R3 and ThinNTF-R3 models, where the suffix R3 shows that the model has 3 components.We sampled the synthetic data with uniform and ANLS sampling methods with different sampling parameters. The sampling was done simply by randomly drawing a sampled size for each flow according to the sampling probabilities in theS matrix. By this way, we ignored the flow splitting problem, and this gave us an ideal data for the ThinNTF model. We report and compare the mean standard deviation of these WMRD values for all experiments.The ThinNTF model always performed best with the uniform sampling model, as shown in Table3 as expected. On ANLS sampling, the MLE and the ANLS estimators performed better with high sampling probabilities, when u∈0.01,0.02,0.05, as shown in Table 4. On the other hand, when the sampling probability of the ANLS model decreases, the ThinNMF helped with better estimations. From the initial results, we conclude that the factorization is definitely helpful for more difficult uniform sampling method and helps lower the sampling probabilities in flow-based packet sampling. The results are also visible in Figure 9.Table 3
Uniform sampling results on synthetic data.
Period
ThinNMF-R3
ThinNTF-R3
MLE
2
0.53
0.49
0.88
4
0.63
0.59
1.20
8
0.65
0.61
1.29
16
0.68
0.61
1.41
32
0.74
0.61
1.37
64
0.85
0.61
1.50Table 4
ANLS sampling results on synthetic data.
U
ThinNMF-R3
ThinNTF-R3
MLE
ANLS
0.01
0.29
0.27
0.09
0.15
0.02
0.31
0.29
0.17
0.27
0.05
0.33
0.32
0.36
0.28
0.1
0.35
0.34
0.51
0.38
0.2
0.38
0.36
0.59
0.67
0.5
0.48
0.47
0.72
0.70Figure 9
Synthetic experiment results.
(a)
(b)


### 6.2. Experiments on Real-World Data

The original data collected from a mobile network provider as we described in Section5 are sampled with both sampling methods. However, this time, we simulated the real network offline by feeding the packets one by one to the monitoring server, as described in Section 5, This way, we were able to create the actual conditions on a sampler installed at a network provider’s backbone. This also created the flow splitting problem, since we applied a 30 seconds timeout in our flow hash. We set the maximum flow length as I=100, meaning that X100,:,: entries show the count of flows that have more than 99 packets. This clamping decision was made according to the cumulative distribution of flow lengths as shown in Figure 8. We also clamped the sampling matrices S so that they exactly match the model.Since the number of components in the original flow distribution is unknown, we run our experiments withR∈2,3,4 components for ThinNMF and ThinNTF. The rest of the experiment is similar to the synthetic one. The sampled Y matrix with shape 100×24×7 is factorized and X^ is reconstructed with the estimated factors. We reported and compared the mean and standard deviation of 24×7 WMRD and KL values.The factorization models, both ThinNMF and ThinNTF helped lower the WMRD score in both uniform and ANLS sampling methods. ThinNTF-R4 model consistently gave lower error than the MLE baseline for uniform model as shown in Table5 and Figure 10. Indeed, our factorization framework improved results overall for uniform sampling. However, since recovering true estimates in uniform sampling is quite difficult, we see less impact of the factorization as the sampling ratio increases.Table 5
Uniform sampling results on real data.
Period
ThinNMF
ThinNTF
MLE
R = 2
R = 3
R = 4
R = 2
R = 3
R = 4
2
0.23
0.24
0.23
0.21
0.25
0.22
0.41
4
0.55
0.52
0.53
0.50
0.48
0.49
0.69
8
0.94
0.93
0.94
0.91
0.90
0.87
0.97
16
1.15
1.11
1.11
1.09
1.05
1.04
1.05
32
1.25
1.24
1.24
1.16
1.13
1.10
1.22
64
1.31
1.29
1.30
1.09
1.06
1.04
1.22Figure 10
Real-world data results with Uniform sampler.
(a)
(b)
(c)
(d)Figure10 also gives the KL values between the true and estimated flow length distributions. While the scale of this metric is different, it gives consistent results with the WMRD. This shows that our model, which minimizes the KL metric, also minimizes the commonly used WMRD metric; hence, the model is suitable for this problem.Another important issue is that for uniform sampling, 3-way factorization is more successful than the 2-way factorization. The periodicity information which is captured by the ThinNTF model helps improve the estimates and makes it a more successful model for this sampling method.In ANLS, all our factorization models gave lower error values than the MLE and unbiased estimator of ANLS as shown in Table6 and Figure 11. Since ANLS is a more powerful sampling method than uniform sampling, the effect framework is slightly less for small sampling parameter u. However, both ThinNMF and ThinNTF gave better result while sampling smaller number of packets (when u is large). Furthermore, since we are trying to recover the same original data in both experiments, we can compare our ThinNMF and ThinNTF models under two sampling methods. We see that, in both methods, as the number of components increases, the models gave lower error rates. However, with the uniform sampling method, 3-dimensional methods give better results, while with ANLS, 2-dimensional models perform slightly better.Table 6
ANLS sampling results on real data.
U
ThinNMF
ThinNTF
MLE
ANLS
R = 2
R = 3
R = 4
R = 2
R = 3
R = 4
0.01
0.03
0.02
0.01
0.05
0.04
0.03
0.05
0.12
0.02
0.04
0.03
0.02
0.06
0.04
0.03
0.08
0.21
0.05
0.04
0.03
0.02
0.07
0.05
0.04
0.13
0.39
0.1
0.06
0.05
0.04
0.08
0.07
0.05
0.17
0.61
0.2
0.08
0.08
0.08
0.10
0.09
0.07
0.21
0.70
0.5
0.13
0.13
0.11
0.16
0.15
0.13
0.33
0.94Figure 11
Real-world data results with ANLS sampler.
(a)
(b)
(c)
(d)


### 6.3. Effect of Clamping

The choice of where to clamp the data can be given by multiple factors. First of all, one can set the clamping valueImax according to a value of special interest. Otherwise, we would like to choose a small Imax so that we deal with a dense tensor and less parameters. On the other hand, we would like to set Imax as high as possible so that the clamped portion of the data is as small as possible.We run the best algorithms found in previous section for uniform and ANLS sampling methods withImax∈25,50,75,100. The WMRD values are given in Figure 12. In both methods, Imax=25 gave relatively poor performance and Imax=100 was generally the best choice. Also the results with Imax≥50 are closer to each other. This is consistent with the graph in Figure 8, where the cumulative flow lengths do not change much after Imax=50. A final remark from this experiment is that as the clamping value increases, estimation becomes harder with small sampling rates. This explains the results in uniform sampling with sampling rate 1/64.Figure 12
Clamping experiments.
(a)
(b)


## 6.1. Experiments on Synthetic Data

We prepared our synthetic experiments to test the validity of our models. In this experiment set, we used the generative model of the ThinNTF model as described in Algorithm3 to generate a small network with maximum flow size I=10,J=24 and K=7. The original synthetic flow length distribution X is generated by a generative model with 3 components, where each component is a column in the F factor. We selected these 3 components as exponential, inverted exponential, and uniform random distributions. Therefore, in experiments, we used ThinNMF-R3 and ThinNTF-R3 models, where the suffix R3 shows that the model has 3 components.We sampled the synthetic data with uniform and ANLS sampling methods with different sampling parameters. The sampling was done simply by randomly drawing a sampled size for each flow according to the sampling probabilities in theS matrix. By this way, we ignored the flow splitting problem, and this gave us an ideal data for the ThinNTF model. We report and compare the mean standard deviation of these WMRD values for all experiments.The ThinNTF model always performed best with the uniform sampling model, as shown in Table3 as expected. On ANLS sampling, the MLE and the ANLS estimators performed better with high sampling probabilities, when u∈0.01,0.02,0.05, as shown in Table 4. On the other hand, when the sampling probability of the ANLS model decreases, the ThinNMF helped with better estimations. From the initial results, we conclude that the factorization is definitely helpful for more difficult uniform sampling method and helps lower the sampling probabilities in flow-based packet sampling. The results are also visible in Figure 9.Table 3
Uniform sampling results on synthetic data.
Period
ThinNMF-R3
ThinNTF-R3
MLE
2
0.53
0.49
0.88
4
0.63
0.59
1.20
8
0.65
0.61
1.29
16
0.68
0.61
1.41
32
0.74
0.61
1.37
64
0.85
0.61
1.50Table 4
ANLS sampling results on synthetic data.
U
ThinNMF-R3
ThinNTF-R3
MLE
ANLS
0.01
0.29
0.27
0.09
0.15
0.02
0.31
0.29
0.17
0.27
0.05
0.33
0.32
0.36
0.28
0.1
0.35
0.34
0.51
0.38
0.2
0.38
0.36
0.59
0.67
0.5
0.48
0.47
0.72
0.70Figure 9
Synthetic experiment results.
(a)
(b)


## 6.2. Experiments on Real-World Data

The original data collected from a mobile network provider as we described in Section5 are sampled with both sampling methods. However, this time, we simulated the real network offline by feeding the packets one by one to the monitoring server, as described in Section 5, This way, we were able to create the actual conditions on a sampler installed at a network provider’s backbone. This also created the flow splitting problem, since we applied a 30 seconds timeout in our flow hash. We set the maximum flow length as I=100, meaning that X100,:,: entries show the count of flows that have more than 99 packets. This clamping decision was made according to the cumulative distribution of flow lengths as shown in Figure 8. We also clamped the sampling matrices S so that they exactly match the model.Since the number of components in the original flow distribution is unknown, we run our experiments withR∈2,3,4 components for ThinNMF and ThinNTF. The rest of the experiment is similar to the synthetic one. The sampled Y matrix with shape 100×24×7 is factorized and X^ is reconstructed with the estimated factors. We reported and compared the mean and standard deviation of 24×7 WMRD and KL values.The factorization models, both ThinNMF and ThinNTF helped lower the WMRD score in both uniform and ANLS sampling methods. ThinNTF-R4 model consistently gave lower error than the MLE baseline for uniform model as shown in Table5 and Figure 10. Indeed, our factorization framework improved results overall for uniform sampling. However, since recovering true estimates in uniform sampling is quite difficult, we see less impact of the factorization as the sampling ratio increases.Table 5
Uniform sampling results on real data.
Period
ThinNMF
ThinNTF
MLE
R = 2
R = 3
R = 4
R = 2
R = 3
R = 4
2
0.23
0.24
0.23
0.21
0.25
0.22
0.41
4
0.55
0.52
0.53
0.50
0.48
0.49
0.69
8
0.94
0.93
0.94
0.91
0.90
0.87
0.97
16
1.15
1.11
1.11
1.09
1.05
1.04
1.05
32
1.25
1.24
1.24
1.16
1.13
1.10
1.22
64
1.31
1.29
1.30
1.09
1.06
1.04
1.22Figure 10
Real-world data results with Uniform sampler.
(a)
(b)
(c)
(d)Figure10 also gives the KL values between the true and estimated flow length distributions. While the scale of this metric is different, it gives consistent results with the WMRD. This shows that our model, which minimizes the KL metric, also minimizes the commonly used WMRD metric; hence, the model is suitable for this problem.Another important issue is that for uniform sampling, 3-way factorization is more successful than the 2-way factorization. The periodicity information which is captured by the ThinNTF model helps improve the estimates and makes it a more successful model for this sampling method.In ANLS, all our factorization models gave lower error values than the MLE and unbiased estimator of ANLS as shown in Table6 and Figure 11. Since ANLS is a more powerful sampling method than uniform sampling, the effect framework is slightly less for small sampling parameter u. However, both ThinNMF and ThinNTF gave better result while sampling smaller number of packets (when u is large). Furthermore, since we are trying to recover the same original data in both experiments, we can compare our ThinNMF and ThinNTF models under two sampling methods. We see that, in both methods, as the number of components increases, the models gave lower error rates. However, with the uniform sampling method, 3-dimensional methods give better results, while with ANLS, 2-dimensional models perform slightly better.Table 6
ANLS sampling results on real data.
U
ThinNMF
ThinNTF
MLE
ANLS
R = 2
R = 3
R = 4
R = 2
R = 3
R = 4
0.01
0.03
0.02
0.01
0.05
0.04
0.03
0.05
0.12
0.02
0.04
0.03
0.02
0.06
0.04
0.03
0.08
0.21
0.05
0.04
0.03
0.02
0.07
0.05
0.04
0.13
0.39
0.1
0.06
0.05
0.04
0.08
0.07
0.05
0.17
0.61
0.2
0.08
0.08
0.08
0.10
0.09
0.07
0.21
0.70
0.5
0.13
0.13
0.11
0.16
0.15
0.13
0.33
0.94Figure 11
Real-world data results with ANLS sampler.
(a)
(b)
(c)
(d)


## 6.3. Effect of Clamping

The choice of where to clamp the data can be given by multiple factors. First of all, one can set the clamping valueImax according to a value of special interest. Otherwise, we would like to choose a small Imax so that we deal with a dense tensor and less parameters. On the other hand, we would like to set Imax as high as possible so that the clamped portion of the data is as small as possible.We run the best algorithms found in previous section for uniform and ANLS sampling methods withImax∈25,50,75,100. The WMRD values are given in Figure 12. In both methods, Imax=25 gave relatively poor performance and Imax=100 was generally the best choice. Also the results with Imax≥50 are closer to each other. This is consistent with the graph in Figure 8, where the cumulative flow lengths do not change much after Imax=50. A final remark from this experiment is that as the clamping value increases, estimation becomes harder with small sampling rates. This explains the results in uniform sampling with sampling rate 1/64.Figure 12
Clamping experiments.
(a)
(b)


## 7. Conclusions

In this work, we introduced a novel nonnegative tensor factorization model called ThinNTF, which extends the classic nonnegative tensor factorization with an additional constant factor that can represent a network packet sampling method. We showed that this model can be employed to improve the current reconstruction algorithms in recovering the original flow length distributions.We tested our model with two different types of sampling methods: the uniform packet sampling method and a flow-based packet sampling method, called ANLS. We described how to use these methods by showing how to build their sampling matrices.In order to test our model, we collected high-volume data from a mobile network provider for a long period in order to observe the periodical behavior of the flow length distribution. In experiments on synthetic and real-world data, our models gave promising results by lowering the estimation errors compared to the baselines of each sampling method. We conclude that our model can be used to decrease estimation errors or to decrease the sampling probabilities without increasing the estimation error.An important issue left as future work is the online execution of the ThinNTF model. Theoretically, the ThinNTF model can be used online once sufficient data from the target network is collected, and the flow length distribution components, i.e., theF factor, are inferred. The power of our model is that this inference can be done directly from the sampled observations. Once the F factor is estimated, for each incoming observation, the corresponding entries in other factors can be inferred by keeping F constant during the inference. Moreover, F can be updated periodically, say weekly, in a sliding window fashion and kept up to date with the networks flow length behavior.



---

*Source: 8458016-2019-08-18.xml*