Per-Flow Cardinality Estimation Based On Virtual LogLog Sketching

Flow cardinality estimation is the problem of estimating the number of distinct elements in a data flow, often with a stringent memory constraint. It has wide applications in network traffic measurement and in database systems. The virtual LogLog algorithm proposed recently by Xiao, Chen, Chen and Ling estimates the cardinalities of a large number of flows with a compact memory. The purpose of this thesis is to explore two new perspectives on the estimation process of this algorithm. Firstly, we propose and investigate a family of estimators that generalizes the original vHLL estimator and evaluate the performance of the vHLL estimator compared to other estimators in this family. Secondly, we propose an alternative solution to the estimation problem by deriving a maximum-likelihood estimator. Empirical evidence from both perspectives suggests the near-optimality of the vHLL estimator for per-flow estimation, analogous to the near-optimality of the HLL estimator for single-flow estimation.

today make it infeasible to spend much time on processing each element. According to [1] and [2], on a typical OC-768 2 backbone network link with 40 Gbps traffic speed, the time available to process each packet is at best about tens of nanoseconds, corresponding to no more than a hundred elementary operations.
Secondly, such a high traffic rate necessitates the need to use on-chip cache memory on network processors, in order to achieve real-time processing and maintain high throughput. But most on-chip caches on processors are made of SRAM, typically only a few megabytes [3]. A naïve approach that maintains a look-up dictionary to record all the distinct elements seen thus far is too costly, since it requires O(n) bits. In some applications, n can be on the order of billions.
Therefore, an ideal algorithm should be simple and quick in processing each flow element and be compact in memory usage. The good news is that it is usually not necessary to know the exact value of a flow's cardinality. Many approximation algorithms have been developed to explore the trade-off between estimation accuracy and space/timing efficiency.
Applications Examples: One important direction of application of cardinality estimation is network anomaly detection. Estan et al. [4] suggested three such applications: detecting port scans, detecting denial of service (DoS) attacks and estimating worm spread rate.
Let us take port scan detection as an example to illustrate the idea. A port scan is a probing process that sends service requests from a client to a range of port addresses of a server with the aim of finding an open port, which an attacker can take advantage of to find vulnerabilities on the server. A router can detect such a process by keeping track of packet flows and using cardinality estimation to measure the number of distinct port addresses attempted by each source. Any source that is trying to connect to an abnormally large number of port addresses in a short time interval should be suspected of conducting a port scan.
Other than network anomaly detection, we can apply flow cardinality estimation to many other problems that have a count-distinct nature by adapting the concept of a flow to contain other (maybe abstract) types of data. For example, Google might be interested in estimating the number of distinct users that query certain keywords in a period of time [3]. Here a per-keyword flow can be defined to be all the user IP addresses that query that keyword. Information about the cardinalities of such flows reflects the popularity of those keywords which might be helpful in the optimization of Google's database and search algorithm. Much research on cardinality estimation was motivated by database-related applications (e.g. [5] and [6]).

Motivation and contributions
State-of-the-art algorithms for cardinality estimation comprise two processes. First, the sketching process reads elements from the flow and stores useful information in a compact data structure. Second, in the estimation process, an estimator takes the recorded information as input and outputs the estimated cardinality of the flow.
The HyperLogLog (HLL) algorithm proposed by Flajolet et al. in [7] is near-optimal and has been widely adopted in many industry systems due to its simplicity and excellent performance; according to [3], it is the best existing algorithm for single-flow estimation. Roughly speaking, the HLL algorithm requires hundreds of bytes to make a fair estimation for a single flow with cardinality up to 4 × 10 9 . Details of the HLL algorithm are presented in Section 2.1.1.
However, in some applications the cardinalities of multiple flows need to be estimated at the same time -this is known as per-flow cardinality estimation. In this case, the rate of hundreds of bytes per flow is still too much in some important scenarios where the number of flows can be on the order of tens of millions and memory usage is critical. To this end, Xiao et al. [3] proposed the virtual LogLog algorithm based on memory sharing at register level, which can potentially bring down the memory cost from hundreds of bytes per flow to one bit per flow on average. Details of the sketching process of this algorithm are described in Section 2.2.1. The estimation process of this algorithm uses an estimator called vHLL, summarized in Section 3.2.1.
The main purpose of this thesis is to propose and investigate alternative estimators for the estimation process of the virtual LogLog algorithm in [3]: • The original vHLL estimator of [3] is based on the HLL estimator for single-flow estimation. We show that the HLL estimator can be generalized by a family of estimators parametrized by a value θ (Section 3.1); for HLL, θ = −1. The idea of this generalization can be easily applied to the vHLL estimator for per-flow estimation as well.
Although it is already known that for single-flow estimation, θ = −1 (i.e. HLL) is near-optimal [7], it is not clear whether this near-optimality at θ = −1 extends to the case of per-flow estimation. We provide empirical evidence to show that indeed θ = −1 (i.e. vHLL) is still optimal for per-flow estimation.
• We introduce an alternative approach to the estimation problem -a maximum-likelihood estimator. We find that this estimator does not significantly improve the estimation process, compared to the vHLL estimator.
Empirical evidence from both perspectives suggests the near-optimality of the vHLL estimator for per-flow estimation, on the basis of the virtual LogLog algorithm [3].

Thesis overview
The rest of this thesis is organized as follows. Chapter 2 provides the background knowledge and preliminaries to the estimators discussed in this thesis. Chapter 3 and Chapter 4 investigate alternative estimators of the virtual LogLog algorithm from the aforementioned two perspectives. Chapter 3 introduces the vLL θ estimator, a generalization of the vHLL estimator, and evaluates its performances for different values of θ. Chapter 4 derives the maximum-likelihood estimator and compares its performance with that of the vLL θ estimator. Chapter 5 summarizes the conclusions of the thesis and points to possible future works.
This chapter briefly overviews cardinality estimation algorithms and points to relevant references. The focus is on summarizing previous works that are essential to the understanding of the main body of this thesis: the LogLog [8] and HyperLogLog [7] algorithms for single-flow estimation and the virtual LogLog algorithm [3] for per-flow estimation.

Single-flow cardinality estimation
Given a single flow of elements x 1 , x 2 , x 3 , · · · , the goal of cardinality estimation is to estimate the number of distinct elements in the flow. There are two general approaches to this problem: sampling and streaming. A sampling-based algorithm collects a small sample of elements from the flow while bypassing other elements and infers the cardinality of the flow from sample information. A streaming-based algorithm, in contrast, processes every element in the flow. Most state-of-the-art algorithms are based on streaming. For a review of these two kinds of algorithms and a comparison between them, see [9].
A streaming-based algorithm for cardinality estimation has three integral components: • A hash function H, mapping each element x i of the flow into a hashed value H(x i ) of a chosen type. The purpose of H is generally twofold: filtering out duplicate elements (relying on the fact that multiple appearances of the same element have the same hash) and providing randomization.
• A compact data structure D, capturing certain low-dimensional statistical information of the flow's cardinality by recording features of the hashed values. In some literature, the data structure is known as a sketch, meaning it is a concise summary of the flow. The process of updating the sketches is called sketching.
• An estimator E that takes the content of D as input and outputs the estimated cardinality, thus n = E(D).
Usually we want the estimator to be unbiased (i.e. E[ n|n] = n). One commonly used metric for the evaluation of the performance of an estimator is the relative standard error √ Var( n) n . For two unbiased estimators using the same amount of memory, the one with a smaller relative standard error is better. For a detailed classification and comparison of existing single-flow cardinality estimation algorithms, see [10,11]. We take a closer look at the LogLog and HyperLogLog algorithms in Section 2.1.1.

LogLog and HyperLogLog
The HyperLogLog (HLL) [7] algorithm is the best-known algorithm in practice, due to its simplicity and good performance. HLL is a successor of LogLog [8]: the two algorithms use the same hash function and data structure (therefore have the same sketching process) and only differ in their choice of estimators.
The data structure of the LogLog/HLL algorithm is a single array of k registers (counters), denoted as S, where S[i] refers to the i th register in the array (or the value stored in it, depending on the context). Suppose we have a hash function H that maps each element of the flow to a sufficiently long bit string. Define a function ρ that takes a bit string and outputs the position of the leftmost 1-bit of the string, e.g. ρ(000010 · · · ) = 5. The sketching process of LogLog/HLL is summarized as follows.
LogLog/HLL sketching process: 1. Initialize all the registers in S with value 0. Let l = log 2 k.
2. For each element x in flow: . Hash x into a sufficiently long bit string.
(ii) Divide the bit string into two parts: Use the prefix l bits to select a register in S.
Use the remaining bits to update the selected register. Here are a few points to note about the above sketching process. First, k is chosen to be some integer power of 2 so that l is an integer and j ∈ {0, 1, · · · , k − 1}. Assume that for a random element x, H(x) is a random bit string with independent and uniform bits. Then j can be regarded as a uniform random variable in {0, · · · , k − 1}, which means that the k registers in S are equally likely to be selected for updating for a randomly given element. Second, the values of the registers in S after the sketching process should not be affected by duplicate elements, nor by the order of the elements appearing in the flow. This is guaranteed by hashing and the max operation in step (iii).
Third, how large can a register's value get? With the previous assumption on the independence and uniformity of the bits in the hash output string, for a random element x we have ρ(q) = i with probability 1 2 i (i ≥ 1). That is, ρ(q) can be regarded as a geometric random variable with parameter 1/2. By the same assumption, we expect about n k distinct elements distributed to any particular register S[j]. Suppose that this number n k is exact, then at the end of the sketching process, S[j] is just the maximum of n k independent Geo(1/2) random variables. According to [8], previous study has shown that the expectation of S[j] is close to log n k with a small additive bias. Therefore, on average each register needs about log log n k = O(log log n) bits to store its value -this is why such a register is also known as a LogLog sketch. Based on this discussion, a rough estimator for n can be n = k2 S[j] . Following this idea, the LogLog estimator [8] replaces S[j] with the average of the k register values: n = LL(S) ηk2 where η is a suitable bias correction factor. The HLL algorithm [7] uses a different estimator: with a different bias correction factor γ. Note that 2 is commonly known as the geometric mean of the values 2 We will introduce two other estimators for single-flow estimation in Section 3.1 (LL θ estimator) and Section 4.4 (maximum-likelihood estimator), respectively.
It has been shown in [8] and [7] that the LogLog and HLL estimators are asymptotically approximately unbiased in the sense that, as n → ∞, E[ n] n is very close to 1 with a practically negligible fluctuation. It has been shown that as n → ∞, the bias correction constants are independent of n. They do depend on k, however; but as k gets large (e.g. k ≥ 64), they can be safely replaced by constants for all practical purposes. Practical values of η and γ are 0.39701 and 0.7213 respectively.
The relative standard error √ Var( n) n is approximately 1.30 √ k for the LogLog estimator [8] and approximately 1.04 √ k for the HLL estimator [7], as n → ∞. Therefore, using the same number of registers (same amount of memory), HLL's estimate is more accurate than that of LogLog. With the HLL algorithm, we can achieve an estimate accuracy (in terms of relative standard error) of about 5% by choosing k = 512. If we use 5 bits for each register (for measuring cardinalities up to 2 2 5 ≈ 4 × 10 9 ), then in total we need 5k = 2560 bits = 320 bytes for one estimation.
According to Section 4 of [7], the HLL algorithm is near-optimal, in the sense that its relative standard error 1.04 √ k is quite close to the lower bound 1 √ k for a wider class of algorithms based on order statistics. For a more detailed discussion of this lower bound, see [2].
One problem of the HLL algorithm is that it is highly biased and inaccurate for the estimation of small cardinalities. Methods to remedy this flaw are considered in Section 4 of [7] and further in [6].

Per-flow cardinality estimation with memory sharing
In many real-world applications we need to estimate the cardinalities of multiple flows at the same time. For example, consider a stream of data packets from many flows observed by a network monitor device (such as a router). Let each packet be abstracted as a 2-tuple (f, x), where f is the IP source address of the packet (the flow ID) and x is the destination address (i.e. an element of the flow). The goal is, at the end of the measuring period, to estimate the number of distinct destination addresses from each given source address, i.e. the cardinality of each flow. For example, if the stream of packets is (A, 2), (C, 9) An immediate idea to solve this problem is to allocate a separate block of memory for each flow and use any of the existing single-flow cardinality estimation algorithms for each flow. Since we do not know the cardinalities of the flows beforehand, it is inevitable to allocate the maximum amount of memory for each flow; with the HLL algorithm, we still need hundreds of bytes per flow. However, in many applications most of the flows have small cardinalities.  In view of this waste of memory, some algorithms have been proposed that allow memory to be shared among flows. Such an algorithm usually combines the following components/ideas: • A memory pool, which is the actual source of memory for all estimations.
• A virtual data structure for each flow. This data structure does not physically exist but is logically constructed by (often randomly) pointing to memory units in the pool. Consequently, difference flows may share some common parts of the memory.
• In the sketching process, for each incoming element, the algorithm applies the sketching process of an existing single-flow estimation algorithm to update the virtual data structure allocated for the corresponding flow.
• After the sketching process, the algorithm uses an estimator to estimate the cardinality of any given flow (with its ID). The estimator needs to take into account the interference among flows due to memory sharing. This estimation process can be done off-line (in comparison to the online sketching process).
Such algorithms may differ from each other in any of the above aspects. A review and comparison of some per-flow cardinality estimation algorithms based on memory sharing can be found in [3]. All these algorithms, however, share memory at bit level (meaning that the basic unit of the memory pool and of the virtual data structure is a bit), which [3] claims to be inherently too noisy.
Another algorithm that shares memory at bit level, which [3] did not mention, is the virtual FM sketches algorithm proposed in [13]. The algorithm constructs multiple virtual bit arrays from a bit pool for each flow and relies on the sketching process of the FM sketches 2 algorithm for single-flow estimation. The algorithm adopts a maximum-likelihood estimator in the estimation process based on bit patterns in the sketches. We will also consider a maximum-likelihood estimator solution in Chapter 4, but for a different sketching scheme. The performance reported in [13] shows that with memory cost of 1 bit per flow, the algorithm can estimate flows with cardinalities up to 3000 and average flow cardinality about 2.5; it can achieve relative standard error of about 20% for flows with cardinalities less than 500 and 10% for flows with cardinalities in the range 500 to 3000. For relatively small flows with cardinalities in the range 100 to 500, this algorithm may outperform the virtual LogLog algorithm with the vHLL estimator [3]; but it does not appear to be as competitive for larger ranges of flow cardinality.
In view of the inherent drawback of bit-level sharing, [3] proposed the virtual LogLog algorithm based on memory sharing at register level and showed through experiments that such an algorithm outperforms previous ones. The data structure and sketching process of this algorithm are presented in the following subsection.

Virtual LogLog register sharing and sketching process
Virtual Data Structure: The virtual LogLog algorithm [3] keeps a memory pool in the form of a register array, denote as R. Denote the size of the array by m, which is usually a very large number. R[j] refers to the j th register in the array (or the value stored in the register, depending on the context). For each flow f , we form a virtual data structure (virtual register array) denoted as R f , which is a logically constructed array of k registers with the i th register denoted as R f [i]. The registers of R f are randomly selected from R by using k independent hash functions G 0 , G 1 , · · · , G k−1 , each mapping the flow ID uniformly to an integer in {0, · · · , m − 1}, i.e.
The k hash functions can be implemented by a single master hash function G as follows: where "|" is the concatenation operator. It should be emphasized that R f does not need to be physically constructed (thus it is "virtual"). A simple example is shown in Figure 2.2 to illustrate the concepts. In the example, say we want to update the third register of flow f 2 , i.e. R f 2 [3], what actually will be updated is R[5] -we do not even need to know where the other registers in R f 2 are. There are three flows and three corresponding virtual register arrays: A caveat should be pointed out here, which the original paper [3] omitted: it could happen that two logically distinct registers of a virtual register array are mapped from the same physical register in the pool, i.e. G i (f ) = G j (f ) for some i = j. Certainly we wish to avoid this situation, but it is tolerable if the number of such "collisions" is very small.
To see how likely it is for such a "collision" to occur, we can consider a bins-and-balls analogous problem: suppose we have k balls and n bins and we throw the balls sequentially, independently and uniformly at random into the bins. In the end how many bins will contain more than one balls? If k is reasonably large, then the number of balls distributed at each bin can be approximated by an independent Poisson random variable X with mean k m (see Chapter 5.4 of [14]). So we have where the last approximation assumes k m. Therefore out of the m bins, we expect to have approximately k 2 m bins that hold more than one ball. Ideally we want k 2 m to be small, which means with high probability the registers in a virtual register array are all mapped from distinct physical registers in the pool. This factor should be included in the design consideration.
Sketching Process: The sketching process of virtual LogLog is almost identical to that of LogLog/HLL (recall from Section 2.1.1): for each packet (f, x) in the stream, we process x to obtain j (the register selector) and ρ(q) (to be compared with the selected register's value); except that the last step (iii) becomes That is, we are treating the virtual register array as the actual register array. By combining (2.3) and (2.4), the above expression can be re-written as Again, it shows that updates are actually made in the physical registers in R. It should now be clear why the algorithm is called virtual LogLog: it is based on virtual register arrays, where each register is a LogLog sketch. After the sketching process, we obtain m register values. The remaining problem is to infer flow cardinalities from these register values: this is the estimation process. We will consider two kinds of estimators in Chapters 3 and 4, respectively.

A per-flow estimator performance metric -Weighted square error
Before we discuss specific estimators for per-flow cardinality estimation, we consider a metric based on weighted square error to evaluate the performance of any given estimator.
Suppose that the incoming data stream contains M flows, with cardinalities n 1 , n 2 , · · · , n M . Given any estimator E, suppose its corresponding estimates for the flows' cardinalities are n 1 , n 2 , · · · , n M . The weighted square error of estimator E is defined as where w(·) is a weight function mapping the cardinality of a flow to a positive number. For two estimators using the same amount of memory, the one with a lower WSE is better.

Weight function
The choice of the weight function depends on the specific application. For example, if each flow is considered equally important regardless of its cardinality, then we can simply let w(n) = 1 for all n. In many other situations, large flows are considered to be more important than small flows, then we want w(n) be an increasing function in n. The weight function used in this thesis is presented and explained as follows. First, assume that the cardinality of a randomly chosen flow can be modeled as a random variable N with pmf p N (n) Pr(N = n). We use the following weight function: . (2.7) The motivation for using this weight function is explained here. The integral is an approximation of the total weight of flows whose cardinalities are in the range [m 1 , m 2 ). With the weight function in (2.7), we have w(n)p N (n) = 1 n ; it is easy to verify that in this case the integral for > 0 is independent of n * . If we plot the curve w(n)p N (n) as a function of n on log scale of n, then the area under the curve should be approximately the same in each decade interval: [1,10), [10,100), [100, 1000), [1000, 10000), etc. In another word, by choosing such a weight function, we put approximately the same total weight to the aggregate of flow cardinalities in each of these intervals.

Zipf model
The weight function in (2.7) can be applied to any given flow cardinality distribution. The particular distribution used for simulation in this thesis is the Zipf distribution. A random variable N following the Zipf(π, n max ) distribution has the following pmf: where n max is an upper bound of the flow cardinality, π is a parameter that controls the shape of the Zipf distribution (π > 0) and C is a constant such that nmax n=1 p N (n) = 1. The adoption of the Zipf model is motivated by the fact that Zipf's law underlies many Internet applications (see [15]). In Figure 2  With N ∼ Zipf(π, n max ), by (2.7) and (2.8), we have where the constant C can be omitted (set to 1) because when we compare the weighted square errors of two estimators, we only care about their relative values, which are unaffected by a constant factor. We remark that the specific estimators we will discuss in Chapters 3 and 4 do not reply on any particular distribution model. The Zipf model here is for two purposes: to complete the definition of the weighted square error for the performance evaluation of estimators and to generate random flow cardinalities for simulation; both purposes are independent of the estimation process.
Finally, we make a note on experimental evaluation of per-flow estimators. Given the close match between the real trace data and simulated trace data with a Zipf distribution, we will evaluate the estimators using simulated trace data. We generated 100 simulated trace files. Each trace file contains 10 6 flows, with the cardinality of each flow randomly and independently generated according to a Zipf(2.25, 10 5 ) distribution. A flow is presented in the file as a collection of packets with the same source address and distinct destination addresses; different flows have different source addresses (flow IDs). For one experiment, we process one such trace file and estimate the cardinalities of all the 10 6 flows in it. For a given estimator, we perform 100 independent experiments and evaluate its performance based on the statistics averaged over the 100 experiments. All experiments are performed with m = 200, 000 and k = 512. If each register uses 5 bits, this setting uses 10 6 bits to measure the cardinalities of 10 6 flows, leading to one bit per flow on average.

Chapter 3 VIRTUAL LOGLOG ESTIMATOR WITH PARAMETER θ
The sketching process of the virtual LogLog algorithm has been described in Section 2.2.1.
After the sketching process, we can offload the register values from the network measuring device for off-line query. The following per-flow estimation problem is to be solved: In this chapter, we first introduce LL θ , a family of generalized LogLog estimators parameterized by θ, for single-flow estimation. Then, with the similar idea from LL θ , we propose vLL θ , a family of generalized virtual LogLog estimators parameterized by θ, for per-flow estimation.

The LL θ estimator for single-flow estimation
The LL θ estimator is a class of estimators parameterized by θ, which unifies and generalizes the LogLog and HLL estimators described in Section 2.1.1.

Generalized mean
We start with the concept of generalized mean. Given k positive numbers x 1 , x 2 , · · · , x k , their generalized mean parameterized by θ is defined for nonzero θ by In the case of θ = 0, we let which is in fact the limit of A θ (x 1 , · · · , x k ) as θ → 0. Note that

Unification and generalization of LogLog and HLL estimators
Recall from Section 2.1.1, the LogLog and HLL estimators use the geometric mean and harmonic mean, respectively. Based on the generalized mean notation, we attempt to unify these two estimators by proposing the LL θ estimator: where ξ is a suitable coefficient. It is desirable if there exists a value of ξ to make the estimator approximately unbiased and for which we can identify such value. Specifically, ξ should not depend on n. We already know that for θ = 0 (i.e. LogLog [8]) and θ = −1 (i.e. HLL [7]) such values of ξ exist (by letting ξ = η and ξ = γ respectively). But we do not know if it is true for other values of θ. The analysis of this estimator for general θ is difficult. 1 We resort to simulations to explore empirical evidence of the existence of this coefficient ξ.
, · · · , 2 S[k−1] for short. In   From the plots we see: • For the same value of θ and k, n kA θ is almost a constant when n gets large (say n > 10000).
• For the same value of θ and n, n kA θ is almost a constant for large k (256, 512 or 1024). We conclude that, at least for θ ∈ [−3.0, 1.0], there exists a value for the aforementioned coefficient ξ that approximately depends only on θ for large n and k. To reflect this dependence of ξ on θ, we denote this coefficient as ξ θ instead.
To investigate the relationship between ξ θ and θ,  We observe an approximately linear relationship between ξ θ and θ by which is plotted as the red line in Figure 3.2. The LL θ estimator is now completely defined by replacing ξ in (3.4) with ξ θ specified in (3.5). The claim that this family of estimators unifies the LogLog and HLL estimators can be verified by checking that ξ 0 ≈ η (for LogLog) and ξ −1 ≈ γ for (HLL). A natural question one may ask next is if there exists an optimal value of θ for the LL θ estimator. To answer this question, we empirically calculate the relative standard error of the LL θ estimator for different values of θ (and for selected k's). The results are plotted in  According to the experiment results, for each fixed value of k, the optimal value of θ for the LL θ estimator is either −1 or −0.9; in the cases where LL −0.9 is the optimal, the difference between LL −1 and LL −0.9 , in terms of relative standard error, is negligible. Recall the LL −1 estimator is the same as the HLL estimator. This verifies the claim that the HLL is near-optimal [7], which we briefly described in Section 2.1.1.

The vLL θ estimator for per-flow estimation
With the LL θ estimator defined, we now introduce the vLL θ estimator for per-flow estimation. We start with a high-level idea of the estimator and the motivation for introducing the parameter θ here.
Just like for single-flow estimation where we infer the cardinality of a flow from its register array, for per-flow estimation we can infer the cardinality of a given flow f from its virtual register array R f . We can directly apply the LL θ estimator on R f to give an estimate of the total number of distinct elements distributed to R f . But this estimate involves the noise brought by other flows that have registers shared with flow f ; so we probably need a different bias correction coefficient to calibrate this rough estimate. While for the single-flow case θ = −1 is near-optimal, it is not immediately clear whether θ = −1 with a suitable bias correction coefficient is still near-optimal in the per-flow case.
Recall that A θ with a lesser value of θ is more robust to abnormally high values in obtaining the mean. For per-flow estimation, due to register sharing, large flows can cause much noise to some small flows by causing abnormally high register values. Therefore, we speculate that an estimator with a lesser value of θ might work better in the per-flow case due to its robustness against bursty noise from large flows. This motivates our investigation of the vLL θ estimator, introduced in the rest of this section.

Review of vHLL estimator and the generalization to vLL θ estimator
We start with a review of the virtual HyperLogLog (vHLL) estimator introduced in [3]. Let n T be the total aggregate cardinality of all flows, i.e. the sum of the cardinalities of all flows in the packet stream. Let n f + be the total number of distinct elements distributed to the register array R f , which include the distinct elements from flow f and those from other flows (we call noise elements). Suggested in [3] is the following approximate relationship between n f , n f + and n T : which can be interpreted as: the number of noise elements received by R f is the total number of elements from flows other than f scaled by k m -the ratio of the number of registers in R f to the number of registers in R. The assumption here is that noise elements are roughly uniformly distributed to all the m registers in the pool, which is a fair approximation when the number of flows and the number of registers for each flow are both sufficiently large [3]. Rearranging (3.6) gives us Now the values of n T and n f + are not directly available so we need their estimates n T and n f + . It has been shown in Section 6.1 of [3] that if n f + and n T are close to n f + k m (n T − n f ) (which approximately equals to n f + by (3.6)) and n T respectively, then n f is an approximately unbiased estimator of n f .
There are two possible methods to estimate n T . The first method is to treat n T as the cardinality of a grand flow -the flow containing all the distinct packets in the stream. In another word, if a packet is abstracted as a (f, x) pair, to estimate n T is to estimate the total number of distinct f − x pairs (e.g. source-destination pair) in the stream. This is a singleflow cardinality estimation problem and, as discussed in Section 2.1, can be solved using the HLL algorithm with an additional few hundreds of bytes, which is negligible compared to the memory for main register pool R. The second method is to use n T = HLL(R) as a rough estimator, based on the assumption that all the elements of the grand flow are distributed approximately uniformly over R. Either method works fine in practice. So n T is relatively easy to obtain. Now, for n f + , the vHLL estimator applies the HLL single-flow estimator on R f to obtain an estimate, i.e. (3.8) Combining the ideas above, the vHLL estimator for flow f is summarized as: We consider a generalization of the vHLL estimator by replacing the above harmonic mean A −1 with a generalized mean A θ . More specifically, given any value θ, let be a rough estimate of n f . Then we calibrate this rough estimate by suitable additive and multiplicative constants α θ and β θ to obtain a better estimate, i.e.
We discuss how to set the values of α θ and β θ in Section 3.2.2.

How to set α θ and β θ
For a given θ, one way to set the values of α θ and β θ is through empirical error minimization, explained as follows. Recall from the end of Section 2.3.2, we have 100 simulated trace files, each can be considered as a training sample. Suppose that in one training sample we have M flows with cardinalities n 1 , n 2 , · · · , n M , and by the vLL θ estimator specified in (3.11) we obtain corresponding estimates n 1 (θ), n 2 (θ), · · · , n M (θ). Then we can find values of α θ and β θ that minimize the weighted square error defined in Section 2.3 on this training sample. That is (n i − (α + β n i (θ))) 2 w(n i ), (3.12) where we let w(n i ) = n π−1 i and π = 2.25. The solution to the above minimization problem is standard: .
Since we have 100 training samples, we can obtain the values of α θ and β θ for each of the samples using the above method and use their average values for the estimator. The values of α θ and β θ obtained in this way are optimal in the sense that they minimize the weighted square error of the training samples (i.e. the empirical error ). The drawback of this approach is obvious. It only works well for the training samples generated according to a specific Zipf distribution; there is no performance guarantee for other data sets. Also, it relies on the specific definition of the weighted square error (including the weight function and the cardinality distribution model), which may not be universal for all applications. We suggest the following alternative and practical choice of α θ and β θ for the vLL θ estimator, which does not reply on any assumption on weight function or cardinality distribution model: where ξ θ can be calculated using (3.5) and n T can be calculated using either of the two methods mentioned in Section 3.2.1 (e.g. n T = HLL(R)). Note that this is a direct generalization of the vHLL estimator in (3.9). Also note that in this case α θ does not depend on θ.
In Figure 3.4, we plot the empirical values of α θ and β θ obtained by these two different methods through experiments on the 100 training samples. The plots show that our suggested values of α θ and β θ are good as they are close to the optimal values of α θ and β θ for these training samples.

Summary of the vLL θ estimator
We summarize the vLL θ estimator here. For a given θ ∈ [−3.0, 1.0], the vLL θ estimator estimates the cardinality of a given flow f by and n T is an estimate of the total aggregate cardinality of all the flows and can be obtained by calculating HLL(R) in practice. Note that it could happen that Equation (3.15) produces an estimate n f < 1, which is impossible in practice. In such case, we can simply reset n f = 1.

Experimental performance evaluation
We ran experiments on the 100 simulated trace files to evaluate the performance of the vLL θ estimator for different values of θ, in terms of the weighted square error defined in Section 2.3. The results are plotted in Figure 3.5. Results for θ ≥ 0.5 are not plotted in the graph because the values are comparatively too large (which means the corresponding vLL θ estimators are bad). According to the experiment results, the vLL −1 estimator has the best performance, which refutes our speculation at the beginning of Section 3.2 that a value of θ lesser than −1 might be optimal for the per-flow estimation case.
A possible explanation of this phenomenon is as follows. With the weight function defined in Section 2.3 and Zipf(2.25, 10 5 ) used for flow cardinality distribution, we have that w(n) ∝ n 1.15 . This means we put much larger weights on large flows compared to small flows. Therefore the weighted square error of an estimator is largely determined by how well the estimator estimates large flows. But the influence of bursty noise on large flows is much weaker than that on small flows; in another word, the signal-to-noise ratio for the cardinality estimation of a large flow is larger than that for small flows. Therefore, as we have discussed for the single-flow case where θ = −1 is near-optimal, for per-flow case θ = −1 is also likely to give good estimates for large flows and hence results in a smaller weighted square error.

Chapter 4 MAXIMUM LIKELIHOOD ESTIMATOR
In this chapter we derive an alternative approach to the same per-flow estimation problem stated at the beginning of Chapter 3: a maximum-likelihood estimator.

Formulation
We model the values of the k registers in R f after the sketching process as a random vector where Z f,i ∈ {0, · · · , r max }. The r max is the maximum value that can be stored in a register. For example, if the register has 5 bits, then r max = 2 5 − 1 = 31. Register values R f = (R f [0], · · · , R f [k − 1]) is an instance of the random vector Z f . Let p n (z 0 , · · · , z k−1 ) be the pmf of Z f (i.e. the joint pmf of the k random variables) given that n f = n, i.e. it is the likelihood function of n: p n (z 0 , · · · , z k−1 ) Pr Z f = (z 0 , · · · , z k−1 ) n f = n . (4.1) The maximum likelihood estimator finds the value of n that maximizes the likelihood function based on the instance of register values in R f . An equivalent formulation is to maximize the natural log of the likelihood function The maximum-likelihood estimator of flow f 's cardinality is then given by n f,ML = arg max n L f (n), n ∈ {1, 2, · · · , n max }, where n max is an upper bound of the cardinality of a flow.

Derivation of the likelihood function
In this section we derive an analytical expression for the likelihood function p n (z 0 , · · · , z k−1 ).

Outline of derivation
Suppose we have two copies of the register pool R: R (1) and R (2) , both initialized with 0's. Consider the following process: 1. Repeat the virtual LogLog sketching process described in Section 2.2.1 for all elements, except those from flow f , with register pool R (1) .
2. Repeat the sketching process only for elements from flow f with register pool R (2) .
3. Merge the register values in R (1) and R (2) by a register-wise max operation.
The register values in the merged pool as described above are exactly the same as those in R after the original sketching process. This is because the register values at the end of the sketching process are not affected by the order of arrival of the elements in the stream. Let R (1) f ) denote the virtual register array for flow f constructed by selecting registers from R (1) (R (2) ), whose indices in R (1) (R (2) ) are the same as those of R f in R. Then define the following random variables: Z the value of an arbitrary register in R (1) .
W n the value of an arbitrary register in R (2) f , conditioned on n f = n.
Z n the value of an arbitrary register in R f , conditioned on n f = n.
Z represents the background noise in our estimation caused by elements from flows other than f . W n represents the impact of elements from flow f itself. We model the values of the k corresponding registers in R (1) f and R (2) f also as random vectors, respectively, i.e.
Then we have f,1 , · · · , Z (1) Here is an outline of the derivation of an analytical expression for p n (z 0 , · · · , z k−1 ): 1. Find the distribution of Z and argue that Z (1) f,1 , · · · , Z (1) f,k−1 are independent. We will see that it is difficult to obtain the exact distribution of Z , but there exists a convenient and close approximation to it.
2. Find the distribution of W n and argue that Z f,k−1 are independent under the assumption that n is large.
3. With the distribution of Z and W n , find the distribution of Z n by (4.4).

Distribution of Z and independence of Z
(1) f,0 , · · · , Z (1) f,k−1 One possible approach to obtain an approximate expression of the distribution of Z is as follows. First assume that the total number of flows (excluding flow f ) is known and a flow's cardinality can be modeled as a random variable (e.g. Zipf). Then we can model the number of distinct elements distributed to an arbitrary register by a random variable Y and calculate the distribution of Y by adding up the influences from all the individual independent flows on that register. This calculation involves a high-dimensional convolution, which is computationally costly. We can perhaps avoid the convolution by appealing to the central limit theory for approximation. However this approximation gives us a continuous distribution for Y , from which we need to derive the CDF of Z that only takes discrete values (and really concentrates on only a few values as we shall see soon). This calculation may be highly inaccurate and is somewhat complicated. Moreover, in practice we usually do not know the total number of flows beforehand and it is undesirable that our estimator relies on any particular cardinality distribution model.
Another approach to approximate the distribution of Z , which circumvents the above difficulties and complications, is the following. Define a random variable Z: Z the value of an arbitrary register in R.
We claim that, for any practical purpose, the distribution of Z can be directly approximated by the distribution of Z, i.e.
where F Z and F Z are the CDFs of Z and Z, respectively. This approximation is based on the following two considerations: • Register values in R (1) differ from those in R at no more than k registers, because flow f can only affect k registers. If k m, the difference in the CDF of Z compared to that of Z is small.
• Assume there are many flows in the stream and n f is small compared to the total cardinalities of all other flows. As packets are being distributed to a register, it becomes harder and harder for the register's value to further increase, because it requires an much less likely (with geometrically decaying probability) hashed value to occur.
The analytical distribution of Z is also difficult to find. However, we can obtain an empirical distribution of Z directly from the register values in the pool at the end of the sketching process: In Figure 4.1 we show a sample distribution of Z generated using real trace data (the same data used for plotting Figure 2.1). From the figure we can see that the values of most of the registers are centered around 2 to 7.
Since Z is the distribution of an arbitrary register value in R (1) , including the k registers in R (2) f . We have Z (1) f,k−1 ∼ Z . Since the k registers are randomly selected from the pool with replacement, their independence is guaranteed by a good choice of the hash function G.

Distribution of W n and independence of Z
(2) f,0 , · · · , Z (2) f,k−1 To find the distribution of W n , we first define a random variable X n : X n the number of distinct elements distributed to an arbitrary register in R (2) f , conditioned on n f = n.
It is easy to see X n ∼ Binom n, 1 k . The number of distinct elements distributed to each of the k registers in R (2) f has the same distribution as X n , but they are not independent of each other (because they sum to n). However, if n is large, we can approximate n by a Poisson random variable with mean n. Under this Poissonization approximation, X n ∼ Poisson n k and the number of distinct elements distributed to each of the k registers in R (2) f will then be independent (see Chapter 5.4 of [14] for this Poissonization trick).
As discussed in Section 2.1.1, W n can be regarded as the maximum of X n i.i.d. geometric random variables with parameter 1/2. Given that X n ∼ Poisson n k , the CDF of W n can be calculated as follows. For 0 ≤ i ≤ r max , Under the Poissonization approximation, Z f,0 , Z f,1 , · · · , Z (2) f,k−1 are independent and have the same distribution with W n . We emphasize that this independence follows from the assumption that n is large.

Distribution of Z n and independence of
Recall that, by construction, Z n = max{Z , W n }. Therefore the CDF of Z n is: where F Z (i) and F Wn (i) are the CDF of Z and W n respectively. We have used that Z and W n are independent of each other and the distribution of Z can be approximated by that of Z. Hence, the pmf of Z n can be approximated by where with Poissonization of n. (4.9) We have argued the independence of Z (1) f,k−1 and the independence of Z (2) f,i by construction. So Z f,0 , · · · , Z f,k−1 are also independent of each other under the same conditions. The consequence of this independence is that we can factor the k-variable likelihood function p n (z 0 , · · · , z k−1 ) into the product of k single-variable pmfs, based on the distribution of Z n . That is, under this approximation, the log likelihood function in (4.2) can be re-written as (4.10)

Implementation
It would be nice if we could obtain a closed-form analytical expression of the value of n that maximizes the log-likelihood function L f (n). Unfortunately it turns out to be difficult. In order to search for this value, we explore properties of L f (n) that might be helpful.

Concavity of the log-likelihood function
In this subsection we show that L f (n) has the decreasing increment property with respect to n. That is, L f (n) is the restriction of a concave function to integer values. A by-product of this property is that L f (n) must have a global maximum over the possible values of n. We will call this property "concave" or "concavity" for convenience henceforth. Suppose for now that the possible values of n is a continuous range. First, since R f [j] ∈ {0, · · · , r max }, we have c i ln p Zn (i), (4.11) for some non-negative integer constants c 0 , c 1 , · · · , c rmax such that rmax i=1 c i = k (i.e. the constants represent the empirical pmf of (R f [0], · · · , R f [k − 1])). Therefore a sufficient condition for L f (n) to be concave in n is that ln p Zn (i) is concave in n for each fixed value i ∈ {0, · · · , r max }. Recall the pmf of Z n from (4.8) and (4.9). Let us denote a i F Z (i) and When i = 0, ln p n (i) = ln a i b n i = ln a i + n ln b i , which is linear in n and thus concave in n. When 1 ≤ i ≤ r max , we prove that ln p Zn (i) is concave in n by showing that its second derivative with respect to n is non-positive: We conclude that L f (n) is concave in n.
In Figure 4.2, we show the plot of L f (n) as a function of n for four selected flows whose cardinalities are 150, 611, 1000 and 31536, respectively. The plots are generated based on the same trace data used for We can also observe that in each plot, the value of n at which the curve peaks is close to the actual flow cardinality n f (though not exact the same due to randomness and noise).

Bi-section search implementation
We have shown that L f (n) is concave in n. Since n can only take integer values in {1, 2, · · · , n max }, we can do a bi-section search to find the optimal n.
The pseudo-code above should be self-explanatory. In actual implementation, it is found that the algorithm reaches the optimal n at a faster rate if we replace mid 1 = lb+ub 2 with mid 1 = √ lb · ub, which is equivalent to log mid 1 = 1 2 (log lb + log ub), i.e. a bi-section search on log scale. This faster rate is because, as we already mentioned, most of the flows have small cardinalities. If n max = 10 6 , then at the first while loop iteration, with mid 1 = lb+ub 2 , we have mid 1 = 5 × 10 5 ; but with mid 1 = √ lb · ub we have mid 1 = 1000 -the latter search is quicker if the optimal n is actually small.
With bi-section search, the number of searches required for the estimation of one flow is upper bounded by log n max . This is not ideal compared to the vLL θ estimator which only needs one search/estimation. However, since the estimation process is off-line, the worst-case time complexity of log n max for searching is still acceptable.

Summary of the maximum-likelihood estimator
We wrap up the ideas in the previous sections of this chapter and name this estimator vLL-MLE, meaning that it uses maximum likelihood estimation based on virtual LogLog sketching.
The estimator is summarized as follows. To estimate the cardinality of flow f : 1. Find the values of the k registers in R f : R f [0], · · · , R f [k − 1].
2. From the register values in R, find the CDF of Z by (4.6).
3. Do a bi-section search on integer value n in the range {1, 2, . . . , n max } to maximize the log likelihood function L f (n), output this value of n as the estimate.

Aside: MLE for single-flow estimation
As an aside, we can also use the vLL-MLE estimator for single-flow estimation, with some minor changes. Recall that in single-flow estimation, the goal is to estimate the flow's cardinality n from its k register values. This is similar to the per-flow estimation case, where the goal is to estimate a given flow f 's cardinality n f from its k register values R f [0], · · · , R f [k − 1], except that in the single-flow case, there is no background noise in the register values. That is, the vLL-MLE estimator for per-flow estimation can be used for single-flow estimation by simply letting Z = 0. In this case, we do not need the assumptions related to Z anymore. But the Poissonization approximation of n (i.e. assuming n is large) is still necessary.
To evaluate the performance of the MLE on single-flow estimation, we run experiments to obtain empirical values of its relative standard error for selected values of k and compare it with the HLL estimator. The experiments are performed with n = 10 7 . The results are summarized in Table 4.1. The experiment results show that the MLE's performance is very close to that of HLL, but not any better. It reinforces the claim that the HLL estimator is near-optimal for single-flow estimation.

Experimental performance evaluation
In this section we evaluate the performance of the vLL-MLE estimator and compare it with that of the vLL θ estimator introduced in Chapter 3. Since we have shown that vLL −1 (i.e. vHLL) is the best within the family of the vLL θ estimators, we compare vLL-MLE with vLL −1 in particular. Results for each estimator are generated by running experiments on the same 100 simulated trace files. Figure 4.3 shows the estimation results of both estimators by directly plotting the estimated cardinalities vs. the corresponding actual cardinalities. Each point in the graphs represents one flow, with its x-coordinate value being the actual cardinality of the flow and its y-coordinate value being the estimated cardinality. The more clustered the points are to the equality line y = x, the more accurate the estimator is. We can see that the two estimators have comparable performances.  ) vs. the actual flow cardinality for both estimators. Since there are not many flows for some cardinalities (especial the large cardinalities), we divide the horizontal axis into bins of width 1000 for cardinalities ≤ 10000 and of width 5000 for cardinalities > 10000. In each bin, we calculate the empirical relative bias and relative standard error of the data and interpolate the values for each bin on the graph to form the plots. Again we see that the two estimators have very similar performances.  In both figures, the curves are plotted for cardinalities larger than 1000 only. This is for better graph layouts: we found from the experiment results that estimates for flows with cardinalities less than 1000 are very inaccurate. Accurate estimation for small flows is difficult because of the noise caused by large flows in their register values. Figure 4.4 shows that both estimators are approximately unbiased for large flows (with cardinalities > 1000). Figure 4.5 shows that in general, for both estimators, the estimation accuracy improves as the cardinality gets larger. We do observe a slight curving up at the high end of the graph though, meaning that the accuracy stops improving once the actual flow cardinality reaches a certain level (probably > 40000). Finally, in Figure 4.6 we compare the weighted square error of the vLL-MLE estimator with that of the vLL θ estimators. Experiment results show that the vLL-MLE estimator outperforms the vLL θ estimator for all values of θ; compared to the vLL −1 estimator, the vLL-MLE estimator has a slight improvement of about 3.5%. We conclude that the two estimators, vLL-MLE and vLL −1 (i.e. vHLL), have comparable performances.
identify large flows and estimate their cardinalities directly and more efficiently are useful in this context.