---
abstract: 'This work presents new cardinality estimation methods for data sets recorded by HyperLogLog sketches. A simple derivation of the original estimator was found, that also gives insight how to correct its deficiencies. The result is an improved estimator that is unbiased over the full cardinality range, is easy computable, and does not rely on empirically determined data as previous approaches. Based on the maximum likelihood principle a second unbiased estimation method is presented which can also be extended to estimate cardinalities of union, intersection, or relative complements of two sets that are both represented as HyperLogLog sketches. Experimental results show that this approach is more precise than the conventional technique using the inclusion-exclusion principle.'
author:
- Otmar Ertl
bibliography:
- 'bibliography.bib'
title: New Cardinality Estimation Methods for HyperLogLog Sketches
---

Introduction
============

Counting the number of distinct elements in a data stream or large datasets is a common problem in big data processing. In principle, finding the number of distinct elements $\symCardinality$ with a maximum relative error $\symError$ in a data stream requires $\symBigO(\symCardinality)$ space [@Alon1999]. However, probabilistic algorithms that achieve the requested precision only with high probability are able to drastically reduce space requirements. Many different probabilistic algorithms have been developed over the past two decades [@Metwally2008; @Ting2014] until a theoretically optimal algorithm was finally found [@Kane2010]. Although this algorithm achieves the optimal space complexity of $\symBigO(\symError^{-2}+\log \symCardinality)$ [@Alon1999; @Indyk2003], it is not very efficient in practice [@Ting2014].

More practicable and already widely used in many applications is the algorithm [@Flajolet2007] with a near-optimal space complexity $\symBigO(\symError^{-2} \log\log\symCardinality +\log \symCardinality)$. A big advantage of the algorithm is that corresponding sketches can be easily merged, which is a requirement for distributed environments. Unfortunately, the originally proposed estimation method has some problems to guarantee the same accuracy over the full cardinality range. Therefore, a couple of variants have been developed to correct the original estimate by empirical means [@Heule2013; @Rhodes2015; @Sanfilippo2014].

An estimator for sketches, that does not rely on empirical data and that significantly improves the estimation error, is the historic inverse probability estimator [@Ting2014; @Cohen2015]. It trades memory efficiency for mergeability. The estimator needs to be continuously updated while inserting elements and the estimate depends on the insertion order. Moreover, the estimator cannot be further used after merging two sketches, which limits its application to single data streams. If this restriction is acceptable, the self-learning bitmap [@Chen2011], which provides a similar trade-off and also needs less space than the original method, could be used alternatively.

Sometimes not only the number of distinct elements but also a sample of them is needed in order to allow later filtering according to some predicate and estimating the cardinalities of corresponding subsets. In this case the k-minimum values algorithm [@Beyer2007; @Cohen2007] is the method of choice. It needs more space than the algorithm, but also allows set manipulations like construction of intersections, relative complements, or unions [@Dasgupta2015]. The latter operation is the only one that is natively supported by sketches. A sketch that represents the set operation result is not always needed. One approach to estimate the corresponding result cardinality directly is based on the inclusion-exclusion principle, which however can lead to large errors, especially if the result is small compared to the input set sizes [@Dasgupta2015]. Therefore, it was proposed to combine sketches with minwise hashing [@Pascoe2013; @Cohen2016], which improves the estimation error, even though at the expense of significant more space consumption. It was recently pointed out without special focus on sketches, that the application of the method to the joint likelihood function of two probabilistic data structures leads to better cardinality estimates for intersections [@Ting2016].

HyperLogLog Data Structure {#sec:hyperloglog_data_structure}
==========================

$\langle \symBitRepA_1, \ldots, \symBitRepA_\symPrecision,\symBitRepB_1,\ldots,\symBitRepB_\symRegRange\rangle_2 \gets$ $(\symPrecision + \symRegRange)$-bit hash value of $\symDataItem$ $\symIndexI \gets 1+ \langle \symBitRepA_1, \ldots, \symBitRepA_\symPrecision\rangle_2$ $\symRegVal \gets \min(\{\symS\mid \symBitRepB_\symS = 1\}\cup {\{\symRegRange+1\}} )$ $\symRegValVariate_\symIndexI \gets\symRegVal$

The algorithm collects information of incoming elements into a very compact sketching data structure, that finally allows to estimate the number of distinct elements. The data structure consists of $\symNumReg = 2^\symPrecision$ registers. All registers start with zero initial value. The insertion of a data element into a data structure requires the calculation of a $(\symPrecision+\symRegRange)$-bit hash value. The leading $\symPrecision$ bits of the hash value are used to select one of the $2^\symPrecision$ registers. Among the next following $\symRegRange$ bits, the position of the first 1-bit is determined which is a value in the range $[1,\symRegRange+1]$. The value $\symRegRange+1$ is used, if all $\symRegRange$ bits are zeros. If the position of the first 1-bit exceeds the current value of the selected register, the register value is replaced. The complete update procedure is shown in \[alg:insert\].

A sketch can be characterized by the parameter pair $(\symPrecision, \symRegRange)$ where the precision parameter $\symPrecision$ controls the relative estimation error which scales like $1/\sqrt{\symNumReg}$ [@Flajolet2007] while $\symRegRange$ defines the possible range for registers which is $\lbrace 0, 1,\ldots,\symRegRange+1\rbrace$. The case $\symRegRange=0$ corresponds to a bit array and shows that the algorithm can be regarded as generalization of linear counting [@Whang1990]. The number of consumed hash value bits $\symPrecision+\symRegRange$ defines the maximum cardinality that can be tracked. Obviously, if the cardinality reaches values in the order of $2^{\symPrecision+\symRegRange}$, hash collisions will become more apparent and the estimation error will increase.

\[alg:insert\] has some properties which are especially useful for distributed data streams. First, the insertion order of elements has no influence on the final sketch state. Furthermore, any two sketches with same parameters $(\symPrecision, \symRegRange)$ representing two different sets can be easily merged. The sketch that represents the union of both sets can be constructed by taking the register-wise maximum values.

The state of a sketch is described by the vector $\boldsymbol{\symRegValVariate} = (\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg)$. Under the assumption of a uniform hash function the inserted elements are distributed over all $\symNumReg$ registers according to a multinomial distribution with equal probabilities $1/\symNumReg$ [@Flajolet2007]. Therefore, any permutation of $\boldsymbol{\symRegValVariate}$ is equally likely for a given cardinality. Thus, the order of register values $\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg$ contains no information about the cardinality which makes the multiset $\lbrace\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg\rbrace$ a sufficient statistic for $\symCardinality$. Since the values of the multiset are all in the range $[0, \symRegRange+1]$, the multiset can also be written as $\lbrace\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg\rbrace = 0^{\symCountVariate_0}1^{\symCountVariate_1}\cdots\symRegRange^{\symCountVariate_{\symRegRange}}(\symRegRange+1)^{\symCountVariate_{\symRegRange+1}}$ where $\symCountVariate_\symRegVal$ is the multiplicity of value $\symRegVal$. As a consequence, the multiplicity vector $\boldsymbol{\symCountVariate} := (\symCountVariate_0,\ldots,\symCountVariate_{\symRegRange+1})$, which corresponds to the register value histogram, is also a sufficient statistic for the cardinality. By definition we have $\sum_{\symRegVal=0}^{\symRegRange+1}\symCountVariate_{\symRegVal}=\symNumReg$.

Poisson Approximation {#sec:poisson_approximation}
---------------------

The multinomial distribution is the reason that register values are statistically dependent and that further analysis is difficult. For simplification a Poisson model can be used [@Flajolet2007], which assumes that the cardinality itself is distributed according to a Poisson distribution $\symCardinality \sim \text{Poisson}(\symPoissonRate)$. Under the Poisson model the register values are independent and identically distributed according to $$\label{equ:register_value_distribution}
\symProbability(\symRegValVariate \leq \symRegVal\vert\symPoissonRate)
=
\begin{cases}
0 & \symRegVal < 0 \\
e^{-\frac{\symPoissonRate}{\symNumReg 2^\symRegVal}} & 0\leq \symRegVal \leq \symRegRange \\
1 & \symRegVal > \symRegRange.
\end{cases}$$

The Poisson approximation makes it easier to find an estimator $\symPoissonRateEstimate = \symPoissonRateEstimate(\boldsymbol{\symRegValVariate})$ for the Poisson rate $\symPoissonRate$ than for the cardinality $\symCardinality$ under the fixed-size model. Depoissonization finally allows to translate the estimates back to the fixed-size model. Assume we have found an estimator $\symPoissonRateEstimate$ for the Poisson rate that is unbiased $\symExpectation(\symPoissonRateEstimate\vert\symPoissonRate) = \symPoissonRate$ for all $\symPoissonRate\geq 0$. This implies $\symExpectation(\symPoissonRateEstimate\vert\symCardinality) = \symCardinality$ as it is the only solution of $$\symExpectation(\symPoissonRateEstimate\vert\symPoissonRate) 
 = 
\sum_{\symCardinality=0}^\infty \symExpectation(\symPoissonRateEstimate\vert\symCardinality) e^{-\symPoissonRate}\frac{\symPoissonRate^\symCardinality}{\symCardinality!}
= \symPoissonRate
\quad
\text{for all}\ \symPoissonRate\geq 0.$$ Hence, the unbiased estimator $\symPoissonRateEstimate$ conditioned on $\symCardinality$ is also an unbiased estimator for $\symCardinality$, which motivates us to use $\symPoissonRateEstimate$ directly as estimator for the cardinality $\symCardinalityEstimate := \symPoissonRateEstimate$. As our results will show later, this Poisson approximation works well over the full cardinality range, even for estimators that are not exactly unbiased.

Original Estimation Approach {#sec:cardinality_estimation}
============================

The original cardinality estimator [@Flajolet2007] is based on the idea that the number of distinct element insertions a register needs to reach the value $\symRegVal$ is proportional to $\symNumReg  2^{\symRegVal}$. Given that, a rough cardinality estimate can be obtained by averaging the values $\lbrace \symNumReg 2^{\symRegValVariate_1},\ldots,\symNumReg2^{\symRegValVariate_\symNumReg}\rbrace$. The harmonic mean was found to work best as it is less sensitive to outliers. The result is the so-called raw estimator given by $$\label{equ:raw_estimator}
\symCardinalityRawEstimate
=
\alpha_\symNumReg
\frac{\symNumReg}
{\frac{1}{\symNumReg 2^{\symRegValVariate_1}}+\ldots+\frac{1}{\symNumReg2^{\symRegValVariate_\symNumReg}}}
= 
\frac{\symAlpha_\symNumReg \symNumReg^2}{\sum_{\symRegVal=0}^{\symRegRange+1}\symCountVariate_\symRegVal 2^{-\symRegVal}}.$$ Here $\alpha_\symNumReg$ is a bias correction factor [@Flajolet2007] which can be well approximated by $\symAlpha_\infty := \lim_{\symNumReg\rightarrow\infty} \symAlpha_\symNumReg = \frac{1}{2\log 2}$ in practice, because the additional bias is negligible compared to the overall estimation error.

To investigate the estimation error of the raw estimator and other estimation approaches discussed in the following, we filled sketches with up to 50 billion unique elements. To speed up computation we assumed a uniform hash function whose values can be simulated by random numbers. We used the Mersenne Twister random number generator with a state size of bits from the C++ standard library.

\[fig:raw\_estimate\] shows the distribution of the relative error of the raw estimator as function of the true cardinality for $\symPrecision=12$ and $\symRegRange=20$. Corrections for small and large cardinalities have been proposed to reduce the obvious bias. For small cardinalities the sketch can be interpreted as bit array by distinguishing between registers with zero and nonzero values. This allows using the linear counting cardinality estimator [@Whang1990] $$\label{equ:linear_counting_estimator}
\symCardinalityEstimate_\text{small} = \symNumReg \log(\symNumReg/\symCountVariate_0).$$ The corresponding estimation error is small for small cardinalities and is shown in \[fig:small\_range\_estimate\]. It was proposed to use this estimator as long as $\symCardinalityRawEstimate \leq \frac{5}{2}\symNumReg$ where the factor $\frac{5}{2}$ was empirically determined [@Flajolet2007]. For large cardinalities in the order of $2^{\symPrecision+\symRegRange}$, for which a lot of registers are already in a saturated state, meaning that they have reached the maximum possible value $\symRegRange + 1$, the raw estimator underestimates cardinalities. For the 32-bit hash value case $(\symPrecision+\symRegRange=32)$, which was considered in [@Flajolet2007], following correction formula was proposed if $\symCardinalityRawEstimate>2^{32}/30\approx\num{1.43e8}$ to take these saturated registers into account $$\label{equ:large_range_estimate}
\symCardinalityEstimate_\text{large}
=
-2^{32}\log(1-\symCardinalityRawEstimate/2^{32}).$$

The relative estimation error of the original method that includes both corrections is shown in \[fig:original\_estimate\] for the case $\symPrecision=12$ and $\symRegRange=20$. Unfortunately, the ranges where the estimation error is small for $\symCardinalityRawEstimate$ and $\symCardinalityEstimate_\text{small}$ do not overlap, which causes the estimation error to be much larger near the transition region. To reduce the error for cardinalities close to this region, it was proposed to correct the bias of $\symCardinalityRawEstimate$. Empirically collected bias correction data can be either stored as set of interpolation points [@Heule2013], as lookup table [@Rhodes2015], or as best-fitting polynomial [@Sanfilippo2014]. However, all these empirical approaches treat the symptom and not the cause.

![Relative error of the raw estimator for $\symPrecision = 12$ and $\symRegRange=20$.[]{data-label="fig:raw_estimate"}](raw_estimate){width="1\columnwidth"}

![Relative error of the linear counting estimator for bitmap size 4096 which corresponds to $\symPrecision = 12$ and $\symRegRange=0$.[]{data-label="fig:small_range_estimate"}](small_range_estimate){width="1\columnwidth"}

![Relative estimation error of the original method for $\symPrecision = 12$ and $\symRegRange=20$.[]{data-label="fig:original_estimate"}](original_estimate){width="1\columnwidth"}

The large range correction formula is not satisfying either as it does not reduce the estimation error but makes it even worse. Instead of underestimating cardinalities, they are now overestimated. Another indication for the incorrectness of the proposed large range correction is the fact that it is not even defined for all possible states. For instance, consider a $(\symPrecision,\symRegRange)$- sketch with $\symPrecision+\symRegRange=32$ for which all registers are equal to the maximum possible value $\symRegRange+1$. The raw estimate would be $\symCardinalityRawEstimate = \symAlpha_\symNumReg 2^{33}$, which is greater than $2^{32}$ and outside of the domain of the large range correction formula. A simple approach to avoid the need of any large range corrections is to extend the operating range of the raw estimator to larger cardinalities. This can be easily accomplished by increasing $\symPrecision+\symRegRange$, which corresponds to using hash values with more bits. Each additional bit doubles the operating range which scales like $2^{\symPrecision+\symRegRange}$. However, in case $\symRegRange \geq 31$ the number of possible register values, which are $\lbrace 0, 1, \ldots, \symRegRange+1\rbrace$, exceeds 32 and is no longer representable by 5 bits. Therefore, it was proposed to use 6 bits per register in combination with 64-bit hash values [@Heule2013]. Even larger hash values are needless in practice, because it is unrealistic to encounter cardinalities of order $2^{64}$.

Derivation of the Raw Estimator {#sec:derivation_raw_estimator}
-------------------------------

To better understand why the raw estimator fails for small and large cardinalities, we start with a brief and simple derivation without the restriction to large cardinalities ($\symCardinality\rightarrow\infty$) and without using complex analysis as in [@Flajolet2007]. Assume that the register values have following cumulative distribution function $$\label{equ:assumed_register_val_distribution}
\symProbability(\symRegValVariate \leq \symRegVal\vert\symPoissonRate) = e^{-\frac{\symPoissonRate}{\symNumReg 2^{ \symRegVal}}}.$$ For now we ignore that this distribution has infinite support and differs from the register value distribution under the Poisson model , whose support is limited to $[0, \symRegRange+1]$. For a random variable $\symRegValVariate$ obeying the expectation of $2^{-\symRegValVariate}$ is given by $$\label{equ:expectation_power_two}
\symExpectation(2^{-\symRegValVariate})
=
\frac{1}{2}
\sum_{\symRegVal = -\infty}^\infty
2^{-\symRegVal}
e^{-\frac{\symPoissonRate}{\symNumReg 2^{\symRegVal}} }
=
\frac{\symAlpha_\infty\,\symNumReg\,\symPowerSeriesFunc\!\left(\log_2\!\left(\symPoissonRate/\symNumReg\right)\right)}{\symPoissonRate},$$ where the function $$\symPowerSeriesFunc(\symX):= \log(2) \sum_{\symRegVal = -\infty}^\infty
2^{\symRegVal+\symX}
e^{-2^{\symRegVal+\symX}}$$ is a smooth periodic oscillating function with mean 1 and an amplitude that can be bounded by $\symEpsPowerSeriesFunc:=\num{9.885e-6}$ as shown in \[fig:power\_series\_func\]. This limit can also be found using Fourier analysis [@Ertl2017].

![The deviation of $\symPowerSeriesFunc(\symX)$ from 1.[]{data-label="fig:power_series_func"}](power-series-function-minus-1){width="\columnwidth"}

For a large ($\symNumReg\rightarrow \infty$) sample $\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg$, that is distributed according to , we asymptotically have $$\symExpectation\!\left(
\frac{1}{2^
{-\symRegValVariate_1}\!+\ldots+2^{-\symRegValVariate_\symNumReg}}
\right)
\underset{\symNumReg\rightarrow \infty}{=}
\frac{1}{\symExpectation
(2^{-\symRegValVariate_1}+\ldots+2^{-\symRegValVariate_\symNumReg})}
=
\frac{1}{\symNumReg \symExpectation(2^{-\symRegValVariate})}.$$ Together with we obtain $$\symPoissonRate
=
\symExpectation\left(
\frac{\symAlpha_\infty\,\symNumReg^2\,\symPowerSeriesFunc(\log_2( \symPoissonRate/\symNumReg))}{2^{-\symRegValVariate_1}+\ldots+2^{-\symRegValVariate_\symNumReg}}
\right)
\quad
\text{for $\symNumReg\rightarrow \infty$}.$$ Therefore, the asymptotic relative bias of $$\symPoissonRateEstimate 
= 
\frac{\symAlpha_\infty\,\symNumReg^2}{2^{-\symRegValVariate_1}+\ldots+2^{-\symRegValVariate_\symNumReg}}$$ is bounded by $\symEpsPowerSeriesFunc$, which makes this statistic a good estimator for the Poisson parameter. It also corresponds to the raw estimator , if the Poisson parameter estimate is used as cardinality estimate as discussed in \[sec:poisson\_approximation\].

Limitations of the Raw Estimator
--------------------------------

The raw estimator is based on two prerequisites. In practice, only the first requiring $\symNumReg$ to be sufficiently large is satisfied. However, the second assuming that the distribution of register values can be approximated by is not always true. A random variable $\symRegValVariate'$ with cumulative distribution function can be transformed into a random variable $\symRegValVariate$ with cumulative distribution function using $$\label{equ:dist_transformation}
\symRegValVariate = \min\!\left(\max\!\left(\symRegValVariate',0\right), \symRegRange+1\right).$$ Therefore, register values $\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg$ can be seen as the result after applying this transformation to a sample $\symRegValVariate'_1,\ldots,\symRegValVariate'_\symNumReg$ from . If all registers values are in the range $[1,\symRegRange]$, they must be identical to the values $\symRegValVariate'_1,\ldots,\symRegValVariate'_\symNumReg$. In other words, the observed register values are also a plausible sample of distribution . Hence, as long as all or at least most register values are in the range $[1,\symRegRange]$, which is the case if $2^\symPrecision \ll \symPoissonRate \ll 2^{\symPrecision
+\symRegRange}$, the approximation of by is valid. This explains why the raw estimator works best for intermediate cardinalities. However, for small and large cardinalities many register values are equal to 0 or $\symRegRange+1$, respectively, which contradicts and ends up in the observed bias.

Improved Estimator
==================

If we knew the values $\symRegValVariate'_1,\ldots,\symRegValVariate'_\symNumReg$ for which transformation led to the observed register values $\symRegValVariate_1,\ldots,\symRegValVariate_\symNumReg$, we would be able to use the raw estimator $$\label{equ:raw_estimate_with_multiplicities}
\symPoissonRateEstimate 
= 
\frac{\symAlpha_\infty\,\symNumReg^2}{\sum_{\symRegVal=-\infty}^\infty \symCountVariate'_\symRegVal 2^{-\symRegVal}}$$ where $\symCountVariate'_\symRegVal := \vert\lbrace \symIndexI\vert \symRegVal = \symRegValVariate'_\symIndexI\rbrace\vert$ are the multiplicities of value $\symRegVal$ in $\lbrace\symRegValVariate'_1,\ldots,\symRegValVariate'_\symNumReg\rbrace$. Due to , the multiplicities $\symCountVariate'_\symRegVal$ and the multiplicities $\symCountVariate_\symRegVal$ for the observed register values have following relationships $$\label{equ:multiplicity_transformation}
\symCountVariate_0 = \textstyle\sum_{\symRegVal = -\infty}^{0} \symCountVariate'_\symRegVal,\
\symCountVariate_\symRegVal =  \symCountVariate'_\symRegVal\ (1\leq\symRegVal\leq\symRegRange),\
\symCountVariate_{\symRegRange+1} = \textstyle\sum_{\symRegVal = \symRegRange + 1}^{\infty} \symCountVariate'_\symRegVal.$$ The idea is now to find estimates $\hat{\symCount}'_\symRegVal$ for all $\symRegVal\in\mathbb{Z}$ and use them as replacements for $\symCountVariate'_\symRegVal$ in . For $\symRegVal\in[1,\symRegRange]$ we can use the trivial estimators $\hat{\symCount}'_\symRegVal := \symCountVariate_\symRegVal$. Estimators for $\symRegVal\leq 0$ and $\symRegVal\geq\symRegRange+1$ can be found by considering the expectation of $\symCountVariate'_\symRegVal$ $$\symExpectation(\symCountVariate'_\symRegVal)
=
\symNumReg\,
\symProbability(\symRegValVariate' = \symRegVal\vert\symPoissonRate)
=
\symNumReg\,e^{-\frac{\symPoissonRate}{\symNumReg 2^{\symRegVal}}}\left(1-e^{-\frac{\symPoissonRate}{\symNumReg 2^{\symRegVal}}}\right).$$ According to we have $\symExpectation(\symCountVariate_0/\symNumReg)=
e^{-\frac{\symPoissonRate}{\symNumReg}}$ and $\symExpectation(1-\symCountVariate_{\symRegRange+1}/\symNumReg)=
e^{-\frac{\symPoissonRate}{\symNumReg 2^\symRegRange}}$ which gives $$\begin{aligned}
\symExpectation(\symCountVariate'_\symRegVal)
&=
\symNumReg\,
(\symExpectation(\symCountVariate_0/\symNumReg))^{2^{-\symRegVal}}
\left(1-\left(\symExpectation(\symCountVariate_0/\symNumReg)\right)^{2^{-\symRegVal}}\right),
\\
&=
\symNumReg\,
(\symExpectation(1-\symCountVariate_{\symRegRange+1}/\symNumReg))^{2^{\symRegRange-\symRegVal}}
\left(1-(\symExpectation(1-\symCountVariate_{\symRegRange+1}/\symNumReg))^{2^{\symRegRange-\symRegVal}}\right).\end{aligned}$$ These two expressions for the expectation suggest to use $$\begin{aligned}
\hat{\symCount}'_\symRegVal
&=
\symNumReg\,
(\symCountVariate_0/\symNumReg)^{2^{-\symRegVal}}
\left(1-(\symCountVariate_0/\symNumReg)^{2^{-\symRegVal}}\right)\quad\text{and}
\\
\hat{\symCount}'_\symRegVal
&=
\symNumReg\,
(1-\symCountVariate_{\symRegRange+1}/\symNumReg)^{2^{\symRegRange-\symRegVal}}
\left(1-(1-\symCountVariate_{\symRegRange+1}/\symNumReg)^{2^{\symRegRange-\symRegVal}}\right)\end{aligned}$$ as estimators for $\symRegVal \leq 0$ and $\symRegVal \geq \symRegRange+1$, respectively. Both estimators conserve the mass of zero-valued and saturated registers, because is satisfied, if $\symCountVariate'_\symRegVal$ is replaced by $\hat{\symCount}'_\symRegVal$. Plugging all these estimators into as replacements for $\symCountVariate'_\symRegVal$ finally gives $$\label{equ:correctedestimator}
\symPoissonRateEstimate 
= 
\frac{\symAlpha_\infty\symNumReg^2}
{
\symNumReg\,\symSmallCorrectionFunc(\symCountVariate_0/\symNumReg) + \sum_{\symRegVal=1}^\symRegRange \symCountVariate_\symRegVal 2^{-\symRegVal} + \symNumReg\, \symLargeCorrectionFunc(1-\symCountVariate_{\symRegRange+1}/\symNumReg) 2^{-\symRegRange}
}$$ which we call the improved estimator. Here $\symNumReg\,\symSmallCorrectionFunc(\symCountVariate_0/\symNumReg)$ and $2\symNumReg\,\symLargeCorrectionFunc(1-\symCountVariate_{\symRegRange+1}/\symNumReg)$ are replacements for $\symCountVariate_0$ and $\symCountVariate_{\symRegRange+1}$ in the raw estimator , respectively. The functions $\symSmallCorrectionFunc$ and $\symLargeCorrectionFunc$ are defined as $$\begin{aligned}
\label{equ:sigma}
\symSmallCorrectionFunc(\symX) &:= 
\symX
+
\sum_{\symRegVal=1}^\infty
\symX^{2^\symRegVal} 2^{\symRegVal-1},
\\
\label{equ:tau}
\symLargeCorrectionFunc(\symX)
&:=
\frac{1}{3}
\left(
1-\symX
-
\sum_{\symRegVal=1}^\infty
\left(
1-
\symX^{2^{-\symRegVal}}
\right)^{\!2}
2^{-\symRegVal}
\right).\end{aligned}$$ We can cross-check the new estimator for the linear counting case with $\symRegRange=0$. Using the identity $$\label{equ:sigma_tau_relationship}
\symSmallCorrectionFunc(\symX) + \symLargeCorrectionFunc(\symX) = \symAlpha_\infty\symPowerSeriesFunc(\log_2(\log(1/\symX)))/\log(1/\symX),$$ we get $$\symPoissonRateEstimate 
= 
\frac{\symAlpha_\infty\symNumReg}{
\symSmallCorrectionFunc(\symCountVariate_0/\symNumReg)
+
\symLargeCorrectionFunc(\symCountVariate_0/\symNumReg)
}
=
\frac{\symNumReg\log\!\left(\symNumReg/\symCountVariate_0\right)}{\symPowerSeriesFunc\!\left(\log_2\!\left(\log\!\left(\symNumReg/\symCountVariate_0\right)\right)\right)}$$ which is as expected almost identical to the linear counting estimator , because $\symPowerSeriesFunc(\symX)\approx 1$.

The new estimator can be directly translated into an estimation algorithm that does not depend on magic numbers or special cases as previous approaches. Since $\symCountVariate_\symRegVal\in\lbrace0,\ldots,\symNumReg\rbrace$, the values for $\symNumReg\,\symSmallCorrectionFunc(\symCountVariate_0/\symNumReg)$ and $\symNumReg\, \symLargeCorrectionFunc(1-\symCountVariate_{\symRegRange+1}/\symNumReg)$ can be precalculated and kept in lookup tables of size $\symNumReg+1$. In this way a complete branch-free cardinality estimation can be realized. However, on-demand calculation of $\symSmallCorrectionFunc$ and $\symLargeCorrectionFunc$ is very fast as well. The series converges quadratically for all $\symX\in[0,1)$ and its terms can be calculated recursively using elementary operations. The case $\symX=1$ needs special handling, because the series diverges and causes an infinite denominator in and therefore a vanishing cardinality estimate. As this case only occurs if all register values are in initial state ($\symCountVariate_0=\symNumReg$), this is exactly what is expected. Apart from the trivial roots at $\symX=0$ and $\symX=1$, the calculation of $\symLargeCorrectionFunc$ is slightly more expensive, because it involves square root evaluations and its series converges only linearly. Since $1-\symX^{2^{-\symRegVal}}\leq -\log(\symX) 2^{-\symRegVal}$ for $\symX\in(0,1]$ its convergence speed is comparable to a geometric series with ratio $1/8$. It is also thinkable to calculate $\symLargeCorrectionFunc$ using the approximation $\symLargeCorrectionFunc(\symX)\approx \symAlpha_\infty/\log(1/\symX)-\symSmallCorrectionFunc(\symX)$ which can be obtained from and $\symPowerSeriesFunc(\symX)\approx 1$. The advantage is that the calculation of $\symSmallCorrectionFunc$ and the additional logarithm is slightly faster than the direct calculation of $\symLargeCorrectionFunc$ using . However, for arguments close to 1 both terms may become very large compared to their difference, which requires special care to avoid numerical cancellation. The calculation of $\symLargeCorrectionFunc$ can be omitted at all, if the parameters are chosen such that $2^{\symPrecision+\symRegRange}$ is much larger than the expected cardinality. The number of saturated registers is negligible in this case ($\symCountVariate_{\symRegRange+1}\approx 0$) and therefore $\symLargeCorrectionFunc(1 - \symCountVariate_{\symRegRange+1}/\symNumReg)\approx\symLargeCorrectionFunc(1) = 0$.

Results
-------

![Relative error of the improved estimator for $\symPrecision = 12$ and $\symRegRange=20$.[]{data-label="fig:raw_corrected_estimation_error_12_20"}](raw_corrected_estimate_12_20){width="1\columnwidth"}

![Standard deviations of the relative error of different cardinality estimators for $\symPrecision = 12$ and $\symRegRange=20$.[]{data-label="fig:stdev_comparison"}](stdev_comparison){width="1\columnwidth"}

\[sec:corrected\_raw\_estimation\_error\] \[fig:raw\_corrected\_estimation\_error\_12\_20\] shows the relative estimation error of the new improved estimator, again based on randomly generated sketches with parameters $\symPrecision=12$ and $\symRegRange=20$. The experimental results show that new estimator is unbiased up to cardinalities of 10 billions which is a clear improvement over the raw estimator (compare \[fig:raw\_estimate\]). The improved estimator beats the precision of existing methods that apply bias correction on the raw estimator [@Heule2013; @Rhodes2015; @Sanfilippo2014]. Based on the simulated data we have empirically determined the bias correction function $\symBiasCorrectionFunc$ that satisfies $\symCardinality = \symExpectation(\symBiasCorrectionFunc(\symCardinalityRawEstimate)\vert\symCardinality)$ for all cardinalities. By definition, the estimator $\symCardinalityRawEstimate':=\symBiasCorrectionFunc(\symCardinalityRawEstimate)$ is unbiased and a function of the raw estimator. Its standard deviation is compared with that of the improved estimator in \[fig:stdev\_comparison\]. For cardinalities smaller than the empirical bias correction approach is not very precise. This is the reason why all previous approaches had to switch over to the linear counting estimator at some point. The standard deviation of the linear counting estimator is also shown in \[fig:stdev\_comparison\]. Obviously, the previous approaches cannot do better than given by the minimum of both curves for linear counting and raw estimator. In practice, the standard deviation of the combined approach is even larger, because the choice between both estimators must be made based on an estimate and not on the true cardinality, for which the intersection point of both curves, which is approximately in \[fig:stdev\_comparison\], would be the ideal transition point. In contrast, the improved estimator performs well over the entire cardinality range.

We also investigated the estimation error of the improved estimator for the case $\symPrecision=22$, $\symRegRange=10$ as shown in \[fig:raw\_corrected\_estimation\_error\_22\_10\]. The standard deviation is by a factor 32 smaller according to the $1/\sqrt{\symNumReg}$ error scaling law. Again, cardinalities up to order $2^{\symPrecision+\symRegRange}\approx\num{4e9}$ can be well estimated, because $\symPrecision+\symRegRange=32$ was kept the same for both investigated configurations. For $\symPrecision=22$ a small oscillating bias becomes apparent, which is caused by approximating the periodic function $\symPowerSeriesFunc$ by a constant (see \[sec:derivation\_raw\_estimator\]).

![Relative error of the improved estimator for $\symPrecision = 22$ and $\symRegRange=10$.[]{data-label="fig:raw_corrected_estimation_error_22_10"}](raw_corrected_estimate_22_10){width="1\columnwidth"}

Maximum Likelihood Estimation {#sec:max_likelihood_estimation}
=============================

We know from \[sec:poisson\_approximation\] that any unbiased estimator for the Poisson parameter is also an unbiased estimator for the cardinality. Moreover, we know that under suitable regularity conditions of the probability mass function the estimator is asymptotically efficient [@Casella2002]. This means, if the number of registers $\symNumReg$ is sufficiently large, the estimator should be unbiased.

The log-likelihood function for given register values $\boldsymbol{\symRegValVariate}$, which are assumed to be distributed according to under the Poisson model, is $$\label{equ:log_likelihood_single}
\log \mathcal{\symLikelihood}(\symPoissonRate\vert\boldsymbol{\symRegValVariate}) = 
\sum_{\symRegVal=1}^{\symRegRange+1}
\log\left(1-e^{-\frac{\symPoissonRate}{\symNumReg 2^{\min(\symRegVal,\symRegRange)}}}\right)
\symCountVariate_\symRegVal
-
\frac{\symPoissonRate}{\symNumReg}\sum_{\symRegVal=0}^{\symRegRange}\frac{\symCountVariate_\symRegVal}{2^\symRegVal}
.$$ After differentiation and multiplying by $\symPoissonRate$ we find that the estimate $\symPoissonRateEstimate$ is the unique root of the monotone decreasing function $$\symFunc(\symPoissonRate)
=
\sum_{\symRegVal=1}^{\symRegRange+1} \frac{\frac{\symPoissonRate}{\symNumReg 2^{\min(\symRegVal,\symRegRange)}}}{e^{\frac{\symPoissonRate}{\symNumReg 2^{\min(\symRegVal,\symRegRange)}}}-1}
\symCountVariate_\symRegVal
-
\frac{\symPoissonRate}{\symNumReg}\sum_{\symRegVal=0}^\symRegRange \frac{\symCountVariate_\symRegVal}{2^\symRegVal}.$$ Since $\symFunc(0)=\symNumReg-\symCountVariate_0\geq 0$ and $\symFunc$ is at least linear decreasing as long as $\symCountVariate_{\symRegRange+1} < \symNumReg$, there exists a unique root. The special case $\symCountVariate_{\symRegRange+1} = \symNumReg$, for which all registers have reached the maximum possible value, the estimate would be positive infinite.

Using $\symFunc(\symPoissonRateEstimate)=0$ and $1-\frac{\symX}{2}\leq \frac{\symX}{e^\symX-1}\leq 1$ we obtain following bounds for $\symPoissonRateEstimate$ $$\label{equ:inequality}
\frac{\symNumReg(\symNumReg-\symCountVariate_0)}
{\symCountVariate_0+\frac{3}{2}\sum_{\symRegVal=1}^{\symRegRange}\frac{\symCountVariate_\symRegVal}{2^\symRegVal} + \frac{\symCountVariate_{\symRegRange+1}}{2^{\symRegRange+1}}}
\leq
\symPoissonRateEstimate
\leq
\frac{\symNumReg (\symNumReg-\symCountVariate_0)}
{\sum_{\symRegVal=0}^{\symRegRange}
\frac{\symCountVariate_\symRegVal}{2^\symRegVal}}.$$ If the cardinality is in the intermediate range, where $\symCountVariate_0=\symCountVariate_{\symRegRange+1}=0$ the lower and the upper bound differ only by a constant factor and both are proportional to the harmonic mean of $\lbrace \symNumReg 2^{\symRegValVariate_1},\ldots,\symNumReg 2^{\symRegValVariate_\symNumReg} \rbrace$. Hence, consequent application of the method would have directly suggested to use a cardinality estimator that is proportional to the harmonic mean without knowing the raw estimator in advance. The history of the algorithm shows that the raw estimator was first found after several attempts using the geometric mean [@Flajolet2007; @Durand2003].

For $\symRegRange=0$, which corresponds to the linear counting case, the estimator can be found by analytical means. The result is exactly the linear counting estimator which makes us optimistic that the method in combination with the Poisson model also works well for the more general case $\symRegRange>0$. Since $f$ is convex, Newton-Raphson iteration and the secant method [@Press2007] will both converge to the root, provided that the function is positive for the chosen starting points. Even though the secant method has the disadvantage of slower convergence, a single iteration is simpler to calculate as it does not require the evaluation of the first derivative. Possible starting points are zero and the lower bound in . The iteration can be stopped, if the relative increment is below a certain limit $\symStopDelta$. Since the expected estimation error scales according to $1/\sqrt{\symNumReg}$, it makes sense to choose $\symStopDelta = \symStopEpsilon/\sqrt{\symNumReg}$ with some constant $\symStopEpsilon$. For our results presented below we used $\symStopEpsilon = 10^{-2}$. In practice, only a handful of iterations are necessary to satisfy the stop criterion.

Results {#sec:maximum_likelihood_estimation_error}
-------

![Relative error of the estimator for $\symPrecision = 12$ and $\symRegRange=20$.[]{data-label="fig:max_likelihood_estimation_error_12_20"}](max_likelihood_estimate_12_20){width="1\columnwidth"}

![Relative error of the estimator for $\symPrecision = 22$ and $\symRegRange=10$.[]{data-label="fig:max_likelihood_estimation_error_22_10"}](max_likelihood_estimate_22_10){width="1\columnwidth"}

We have investigated the estimation error of the estimator for both configurations as for the improved estimator. \[fig:max\_likelihood\_estimation\_error\_12\_20,fig:max\_likelihood\_estimation\_error\_22\_10\] look very similar to \[fig:raw\_corrected\_estimation\_error\_12\_20,fig:raw\_corrected\_estimation\_error\_22\_10\], respectively. For $\symPrecision=12$, $\symRegRange=20$ the estimator has a somewhat smaller median bias for cardinalities around . In addition, the standard deviation of the relative error is marginally better than that of the improved estimator for cardinalities between and as shown in \[fig:stdev\_comparison\]. For $\symPrecision=22$, $\symPrecision=10$ we see that the the estimator does not have the small oscillating bias as the improved estimator. Since all observed improvements are not very relevant in practice, the improved estimator should be preferred over the method, because it leads to a simpler and faster algorithm.

Joint Estimation {#sec:cardinality_estimation_set_intersections}
================

While the union of two sets that are represented by sketches can be straightforwardly obtained by taking the register-wise maximums, the computation of cardinalities for other set operations like intersections and relative complements is more challenging. The conventional approach uses the inclusion-exclusion principle $$\label{equ:conventional_approach}
\begin{aligned}
&\left\vert\symSetS_1 \setminus \symSetS_2\right\vert
=
\left\vert\symSetS_1 \cup \symSetS_2\right\vert
-
\left\vert\symSetS_2\right\vert,
\quad
\left\vert\symSetS_2 \setminus \symSetS_1\right\vert
=
\left\vert\symSetS_1 \cup \symSetS_2\right\vert
-
\left\vert\symSetS_1\right\vert,
\\
&\left\vert\symSetS_1 \cap \symSetS_2\right\vert
=  
\left\vert\symSetS_1\right\vert  
+
\left\vert\symSetS_2\right\vert
-
\left\vert\symSetS_1 \cup \symSetS_2\right\vert.
\end{aligned}$$ It allows to express set operation cardinalities in terms of the union cardinality. However, this approach can lead to very large estimation errors, especially if the result is small compared to the operand cardinalities [@Dasgupta2015]. In the worst case, the estimate could be negative without artificial restriction to nonnegative values.

Motivated by the good results we have obtained for a single sketch using the method in combination with the Poisson approximation, we applied the same approach also for the estimation of set operation result sizes. Assume two given sketches with register values $\boldsymbol{\symRegValVariate}_1=(\symRegValVariate_{11},\ldots,\symRegValVariate_{1\symNumReg})$ and $\boldsymbol{\symRegValVariate}_2=(\symRegValVariate_{21},\ldots,\symRegValVariate_{2\symNumReg})$ representing sets $\symSetS_1$ and $\symSetS_2$, respectively. The goal is to find estimates for the cardinalities of the pairwise disjoint sets $\symSetA = \symSetS_1\setminus\symSetS_2$, $\symSetB = \symSetS_2\setminus\symSetS_1$, and $\symSetX = \symSetS_1\cap\symSetS_2$. The Poisson approximation allows us to assume that pairwise distinct elements are inserted into the sketches representing $\symSetS_1$ and $\symSetS_2$ at rates $\symPoissonRate_\symSetASuffix$ and $\symPoissonRate_\symSetBSuffix$, respectively. Furthermore, we assume that further unique elements are inserted into both sketches simultaneously at rate $\symPoissonRate_\symSetXSuffix$. We expect that good estimates $\symPoissonRateEstimate_\symSetASuffix$, $\symPoissonRateEstimate_\symSetBSuffix$, and $\symPoissonRateEstimate_\symSetXSuffix$ for the rates are also good estimates for the cardinalities $\vert\symSetA\vert$, $\vert\symSetB\vert$, and $\vert\symSetX\vert$.

Joint Log-Likelihood Function
-----------------------------

In order to apply the method, we need to find the joint probability distribution of both sketches. Under the Poisson model the individual registers are independent and identically distributed. Therefore, we first derive the joint probability distribution for a single register that has value $\symRegValVariate_1$ in the first sketch that represents $\symSetS_1$ and value $\symRegValVariate_2$ in the second that represents $\symSetS_2$. The first one can be thought to be constructed by merging two sketches representing $\symSetA$ and $\symSetX$, respectively. Analogously, the sketch for $\symSetS_2$ could have been obtained from sketches for $\symSetB$ and $\symSetX$. Let $\symRegValVariate_\symSetASuffix$, $\symRegValVariate_\symSetBSuffix$, and $\symRegValVariate_\symSetXSuffix$ be the values of the considered register in the sketches for $\symSetA$, $\symSetB$, and $\symSetX$, respectively. The corresponding values in sketches for $\symSetS_1$ and $\symSetS_2$ are given by $\symRegValVariate_1 = \max(\symRegValVariate_\symSetASuffix, \symRegValVariate_\symSetXSuffix)$ and $\symRegValVariate_2 = \max(\symRegValVariate_\symSetBSuffix, \symRegValVariate_\symSetXSuffix)$. Since $\symSetA$, $\symSetB$, and $\symSetX$ are disjoint and therefore $\symRegValVariate_\symSetASuffix$, $\symRegValVariate_\symSetBSuffix$, and $\symRegValVariate_\symSetXSuffix$ are independent, the joint cumulative probability function of $\symRegValVariate_1$ and $\symRegValVariate_2$ is $$\begin{gathered}
\label{equ:joint_cdf}
\symProbability(
\symRegValVariate_1 \leq \symRegVal_1
\wedge
\symRegValVariate_2 \leq \symRegVal_2
)
=
\\
\symProbability(
\symRegValVariate_\symSetASuffix \leq \symRegVal_1)
\,
\symProbability(
\symRegValVariate_\symSetBSuffix \leq \symRegVal_2)
\,
\symProbability(
\symRegValVariate_\symSetXSuffix \leq \min(\symRegVal_1, \symRegVal_2)
)
=
\\
\begin{cases}
0 & \symRegVal_1 < 0 \vee \symRegVal_2 < 0
\\
e^{
-
\frac{\symPoissonRate_{\symSetASuffix}}{\symNumReg 2^{\symRegVal_1}}
-
\frac{\symPoissonRate_{\symSetBSuffix}}{\symNumReg 2^{\symRegVal_2}}
-
\frac{\symPoissonRate_{\symSetXSuffix}}{\symNumReg 2^{\min(\symRegVal_1, \symRegVal_2)}}
}
& 0\leq\symRegVal_1 \leq \symRegRange \wedge 0\leq\symRegVal_2\leq\symRegRange
\\
e^{
-
\frac{\symPoissonRate_{\symSetBSuffix} + \symPoissonRate_{\symSetXSuffix}}{\symNumReg 2^{\symRegVal_2}}
}
& 0\leq\symRegVal_2 \leq \symRegRange < \symRegVal_1
\\
e^{
-
\frac{\symPoissonRate_{\symSetASuffix} + \symPoissonRate_{\symSetXSuffix}}{\symNumReg 2^{\symRegVal_1}}
}
&  0\leq\symRegVal_1 \leq \symRegRange < \symRegVal_2
\\
1
&
\symRegRange < \symRegVal_1 \wedge \symRegRange < \symRegVal_2.
\end{cases}\end{gathered}$$ Here we used that $\symRegValVariate_\symSetASuffix$, $\symRegValVariate_\symSetBSuffix$, and $\symRegValVariate_\symSetXSuffix$ follow under the Poisson model. The corresponding probability mass function is $$\begin{gathered}
\symProbabilityMass(\symRegVal_1,\symRegVal_2)
=
\symProbability(
\symRegValVariate_1 \leq \symRegVal_1
\wedge
\symRegValVariate_2 \leq \symRegVal_2
)
-
\symProbability(
\symRegValVariate_1 \leq \symRegVal_1-1
\wedge
\symRegValVariate_2 \leq \symRegVal_2
)
\\
-\symProbability(
\symRegValVariate_1 \leq \symRegVal_1
\wedge
\symRegValVariate_2 \leq \symRegVal_2-1
)
+\symProbability(
\symRegValVariate_1 \leq \symRegVal_1-1
\wedge
\symRegValVariate_2 \leq \symRegVal_2-1
).\end{gathered}$$ Since the values for different registers are independent under the Poisson model, the joint probability mass function for all registers in both sketches is $\symProbabilityMass(\boldsymbol{\symRegVal}_1,\boldsymbol{\symRegVal}_2)
=
\prod_{\symIndexI = 1}^{\symNumReg}
\symProbabilityMass(\symRegVal_{1\symIndexI},\symRegVal_{2\symIndexI})$. The estimates $\symPoissonRateEstimate_\symSetASuffix$, $\symPoissonRateEstimate_\symSetBSuffix$, and $\symPoissonRateEstimate_\symSetXSuffix$ can be finally obtained by maximizing the log-likelihood function given by $$\begin{gathered}
\label{equ:log_likelihood_pair}
\log \symLikelihood(
\symPoissonRate_\symSetASuffix,
\symPoissonRate_\symSetBSuffix,
\symPoissonRate_\symSetXSuffix
\vert
\boldsymbol{\symRegValVariate}_1,
\boldsymbol{\symRegValVariate}_2
)
=
\sum_{\symIndexI = 1}^{\symNumReg}
\log(\symProbabilityMass(\symRegValVariate_{1\symIndexI},\symRegValVariate_{2\symIndexI}))
=
\\
\sum_{\symRegVal=1}^{\symRegRange}
\log\left(1-e^{-\frac{\symPoissonRate_\symSetASuffix+\symPoissonRate_\symSetXSuffix}{\symNumReg 2^{\symRegVal}}}\right)
\symCountVariate^{<}_{1\symRegVal}
+
\log\left(1-e^{-\frac{\symPoissonRate_\symSetBSuffix+\symPoissonRate_\symSetXSuffix}{\symNumReg 2^{\symRegVal}}}\right)
\symCountVariate^{<}_{2\symRegVal}
\\
+
\sum_{\symRegVal=1}^{\symRegRange+1}
\log\left(1-e^{-\frac{\symPoissonRate_\symSetASuffix}{\symNumReg 2^{\min\left(\symRegVal,\symRegRange\right)}}}\right)
\symCountVariate^{>}_{1\symRegVal}
+
\log\left(1-e^{-\frac{\symPoissonRate_\symSetBSuffix}{\symNumReg 2^{\min\left(\symRegVal,\symRegRange\right)}}}\right)
\symCountVariate^{>}_{2\symRegVal}
\\
+
\sum_{\symRegVal=1}^{\symRegRange+1}
\log\left(
1
-e^{-\frac{\symPoissonRate_\symSetASuffix+\symPoissonRate_\symSetXSuffix}{\symNumReg 2^{\min\left(\symRegVal,\symRegRange\right)}}}
-
e^{-\frac{\symPoissonRate_\symSetBSuffix+\symPoissonRate_\symSetXSuffix}{\symNumReg 2^{\min\left(\symRegVal,\symRegRange\right)}}}
+
e^{-\frac{\symPoissonRate_\symSetASuffix+\symPoissonRate_\symSetBSuffix+\symPoissonRate_\symSetXSuffix}{\symNumReg 2^{\min\left(\symRegVal,\symRegRange\right)}}}
\right)
\symCountVariate^{=}_{\symRegVal}
\\
-
\frac{\symPoissonRate_\symSetASuffix}{\symNumReg}
\sum_{\symRegVal=0}^{\symRegRange}
\frac{
  \symCountVariate^{<}_{1\symRegVal}+
  \symCountVariate^{=}_{\symRegVal}+
  \symCountVariate^{>}_{1\symRegVal}
}{2^{\symRegVal}}
-
\frac{\symPoissonRate_\symSetBSuffix}{\symNumReg}
\sum_{\symRegVal=0}^{\symRegRange}
\frac{
  \symCountVariate^{<}_{2\symRegVal}+
  \symCountVariate^{=}_{\symRegVal}+
  \symCountVariate^{>}_{2\symRegVal}
}{2^{\symRegVal}}
\\
-
\frac{\symPoissonRate_\symSetXSuffix}{\symNumReg}
\sum_{\symRegVal=0}^{\symRegRange}
\frac{
  \symCountVariate^{<}_{1\symRegVal}+
  \symCountVariate^{=}_{\symRegVal}+
  \symCountVariate^{<}_{2\symRegVal}
}{2^{\symRegVal}}\end{gathered}$$ where $\symCountVariate^{<}_{1\symRegVal}$, $\symCountVariate^{>}_{1\symRegVal}$, $\symCountVariate^{<}_{2\symRegVal}$, $\symCountVariate^{>}_{2\symRegVal}$, and $\symCountVariate^{=}_{\symRegVal}$ are defined as $$\label{equ:sufficient_joint_statistic}
\begin{aligned}
\symCountVariate^{<}_{1\symRegVal}&:=\left|\lbrace\symIndexI\vert\symRegVal=\symRegValVariate_{1\symIndexI}<
\symRegValVariate_{2\symIndexI}\rbrace\right|,&
\symCountVariate^{>}_{1\symRegVal}&:=\left|\lbrace\symIndexI\vert\symRegVal=\symRegValVariate_{1\symIndexI}>
\symRegValVariate_{2\symIndexI}\rbrace\right|,
\\
\symCountVariate^{<}_{2\symRegVal}&:=\left|\lbrace\symIndexI\vert\symRegVal=\symRegValVariate_{2\symIndexI}<
\symRegValVariate_{1\symIndexI}\rbrace\right|,&
\symCountVariate^{>}_{2\symRegVal}&:=\left|\lbrace\symIndexI\vert\symRegVal=\symRegValVariate_{2\symIndexI}>
\symRegValVariate_{1\symIndexI}\rbrace\right|,
\\
\symCountVariate^{=}_{\symRegVal}&:=\left|\lbrace\symIndexI\vert\symRegVal=\symRegValVariate_{1\symIndexI}=
\symRegValVariate_{2\symIndexI}\rbrace\right|
\end{aligned}$$ for $0\leq \symRegVal \leq \symRegRange+1$. These $5(\symRegRange+2)$ values represent a sufficient statistic for estimating $\symPoissonRate_\symSetASuffix$, $\symPoissonRate_\symSetBSuffix$, and $\symPoissonRate_\symSetXSuffix$ and greatly reduce the number of terms and also the evaluation costs of the log-likelihood function. The derived formula is the generalization of to two sketches and therefore has a similar structure.

Numerical Optimization
----------------------

The estimates $\symPoissonRateEstimate_\symSetASuffix$, $\symPoissonRateEstimate_\symSetBSuffix$, and $\symPoissonRateEstimate_\symSetXSuffix$ are obtained by maximizing . Since the three parameters are all nonnegative, this is a constrained optimization problem. The transformation $\symPoissonRate=e^{\symPhi}$ helps to get rid of these constraints and also translates relative accuracy limits into absolute ones, because $\Delta\symPhi = \Delta\symPoissonRate/\symPoissonRate$. Many optimizer implementations allow the definition of absolute limits rather than relative ones. Quasi-Newton methods are commonly used to find the maximum of such multi-dimensional functions. They all require the computation of the gradient which can be straightforwardly derived for . Among these methods the algorithm [@Press2007] is very popular. We used the implementation provided by the Dlib C++ library [@King2009] for our experiments. Good initial guess values are important to ensure fast convergence for any optimization algorithm. An obvious choice are the cardinality estimates obtained by application of the inclusion-exclusion principle . However, in order to ensure that their logarithms are all defined, we require that the initial values are not smaller than 1. The optimization loop is continued until the relative changes of $\symPoissonRate_\symSetASuffix$, $\symPoissonRate_\symSetBSuffix$, and $\symPoissonRate_\symSetXSuffix$ are all smaller than a predefined threshold. For the results presented below we used again $\symStopDelta=\symStopEpsilon/\sqrt{\symNumReg}$ with $\symStopEpsilon = 10^{-2}$. A few tens of iterations are typically necessary to satisfy this stop criterion for starting points determined by the inclusion-exclusion principle.

Results {#results-1}
-------

To evaluate the estimation error of the new joint cardinality estimation approach we have randomly generated independent pairs of sketches for which the relative complement cardinalities $|\symSetA|$ and $|\symSetB|$ and the intersection cardinality $|\symSetX|$ are known. Each pair is constructed by randomly generating three sketches filled with $|\symSetA|$, $|\symSetB|$, and $|\symSetX|$ distinct elements, respectively. Then we merged the first with the third and the second with the third to get sketches for $\symSetS_1 = \symSetA\cup\symSetX$ and $\symSetS_2=\symSetB\cup\symSetX$, respectively.

\[tbl:joint\_estimation\_results\] shows the results for sketches with parameters $\symPrecision=16$ and $\symRegRange=16$. For different cardinality configurations $|\symSetA|$, $|\symSetB|$, and $|\symSetX|$ we have compared the conventional approach using the inclusion-exclusion principle with the new joint approach. Among the considered cases there are also cardinalities that are small compared to the number of registers in order to prove that the new approach also covers the small cardinality range where many registers are still in initial state. We have determined the relative from the estimates for $|\symSetA| = |\symSetS_1\setminus\symSetS_2|$, $|\symSetB| = |\symSetS_2\setminus\symSetS_1|$, and $|\symSetX| = |\symSetS_1\cap\symSetS_2|$ for all generated examples. In addition, we investigated the estimation error of $\symPoissonRateEstimate_\symSetASuffix+\symPoissonRateEstimate_\symSetBSuffix
+
\symPoissonRateEstimate_\symSetXSuffix$ for the union cardinality $|\symSetU| = |\symSetS_1\cup\symSetS_2|$. We also calculated the improvement factor which we defined as the ratio between both approaches. Since we only observed values greater than one, the new estimation approach improves the precision for all investigated cases. For some cases the improvement factor is even clearly greater than two. Due to the square root scaling of the error, this means that we would need four times more registers to get the same error when using the conventional approach. As the results suggest the new method works well over the full cardinality range without the need of special handling of small or large cardinalities. Obviously, the joint estimation algorithm is also able to reduce the estimation error for unions by a significant amount.

A reason why the method performs better than the inclusion-exclusion method is that the latter only uses a fraction of the available information given by the sufficient statistic , because the corresponding estimator can be expressed as a function of just the three vectors $(\boldsymbol\symCountVariate^{<}_{1} + \boldsymbol\symCountVariate^{=} + \boldsymbol\symCountVariate^{>}_{1})$, $(\boldsymbol\symCountVariate^{<}_{2} + \boldsymbol\symCountVariate^{=} + \boldsymbol\symCountVariate^{>}_{2})$, and $(\boldsymbol\symCountVariate^{>}_{1} + \boldsymbol\symCountVariate^{=} + \boldsymbol\symCountVariate^{>}_{2})$. In contrast, the method uses all the information as it incorporates each individual value of the sufficient statistic.

Outlook {#sec:future_work}
=======

As we have shown, the method is able to improve the cardinality estimates for results of set operations between two sketches. Unfortunately, joint cardinality estimation is much more expensive than for a single sketch, because it requires maximization of a multi-dimensional function. Since we have found the improved estimator which is almost as precise as the estimator for the single sketch case, we could imagine that there also exists a faster algorithm for joint cardinality estimation of two sketches. It is expected that such a new algorithm makes use of all the information given by the sufficient statistic .

The method can also be used to estimate distance measures such as the Jaccard distance of two sets that are represented as sketches. This directly leads to the question whether the algorithm could be used for locality-sensitive hashing [@Leskovec2014; @Wang2014]. The algorithm itself can be regarded as hashing algorithm as it maps sets to register values. For sufficiently large cardinalities we can use the Poisson approximation and assume that the number of zero-valued registers can be ignored. Furthermore, if $\symPrecision+\symRegRange$ is chosen large enough, the number of saturated registers can be ignored as well. As a consequence, we can simplify and assume that $\symRegValVariate_1$ and $\symRegValVariate_2$ are distributed according to $$\symProbability(
\symRegValVariate_1 \leq \symRegVal_1
\wedge
\symRegValVariate_2 \leq \symRegVal_2
)
=
e^{
-
\frac{\symPoissonRate_{\symSetASuffix}}{\symNumReg 2^{\symRegVal_1}}
-
\frac{\symPoissonRate_{\symSetBSuffix}}{\symNumReg 2^{\symRegVal_2}}
-
\frac{\symPoissonRate_{\symSetXSuffix}}{\symNumReg 2^{\min(\symRegVal_1, \symRegVal_2)}}
}$$ which yields $$\begin{gathered}
\symProbability(\symRegValVariate_1=\symRegValVariate_2)
=
\sum_{\symRegVal=-\infty}^\infty \symProbability(\symRegValVariate_1 = \symRegVal \wedge
\symRegValVariate_2 = \symRegVal)
=
\\
\sum_{\symRegVal = -\infty}^{\infty}
e^{-\frac{\symPoissonRate_\symSetASuffix + \symPoissonRate_\symSetBSuffix + \symPoissonRate_\symSetXSuffix}
{\symNumReg 2^\symRegVal}
}
\left(
1
-
e^{-\frac{\symPoissonRate_\symSetASuffix +  \symPoissonRate_\symSetXSuffix}
{\symNumReg 2^\symRegVal}
}
-
e^{-\frac{\symPoissonRate_\symSetBSuffix + \symPoissonRate_\symSetXSuffix}
{\symNumReg 2^\symRegVal}
}
+
e^{-\frac{\symPoissonRate_\symSetASuffix + \symPoissonRate_\symSetBSuffix + \symPoissonRate_\symSetXSuffix}
{\symNumReg 2^\symRegVal}
}
\right).\end{gathered}$$ for the probability that a register has the same value in both sketches. Using the approximation $\sum_{\symRegVal=-\infty}^\infty
e^{-\frac{\symX}{2^\symRegVal}} - e^{-\frac{\symY}{2^\symRegVal}}
\approx
2\symAlpha_\infty (\log(\symY)-\log(\symX))$, which is obtained by integrating $\symPowerSeriesFunc(\symX)\approx 1$ on both sides, we get $$\symProbability(\symRegValVariate_1=\symRegValVariate_2)
\approx
1
+
2\symAlpha_\infty
\log\left(
1
-
\textstyle\frac{1}{2}
\symDistanceMeasure
+
\textstyle\frac{1}{4}
\symDistanceMeasure^2
\textstyle\frac{
\symPoissonRate_\symSetASuffix
\symPoissonRate_\symSetBSuffix
}{\left(
\symPoissonRate_\symSetASuffix+
\symPoissonRate_\symSetBSuffix
\right)^2}
\right)$$ where $\symDistanceMeasure = \frac{\symPoissonRate_\symSetASuffix
+\symPoissonRate_\symSetBSuffix}
{\symPoissonRate_\symSetASuffix+\symPoissonRate_\symSetBSuffix+
\symPoissonRate_\symSetXSuffix}$ corresponds to the Jaccard distance. Since $
\frac{
\symPoissonRate_\symSetASuffix
\symPoissonRate_\symSetBSuffix
}{\left(
\symPoissonRate_\symSetASuffix+
\symPoissonRate_\symSetBSuffix
\right)^2}
$ is always in the range $[0,\frac{1}{4}]$, the probability for equal register values can be bounded by $
1
+
2\symAlpha_\infty
\log(
1
-
\textstyle\frac{1}{2}
\symDistanceMeasure
)
\lesssim
\symProbability(\symRegValVariate_1=\symRegValVariate_2)
\lesssim
1
+
2\symAlpha_\infty
\log(
1
-
\textstyle\frac{1}{2}
\symDistanceMeasure
+
\textstyle\frac{1}{16}
\symDistanceMeasure^2
)
$ as shown in \[fig:equal\_register\_probability\].

![The approximate probability range of equal register values as a function of the Jaccard distance.[]{data-label="fig:equal_register_probability"}](equal-register-probability){width="\columnwidth"}

This dependency on the Jaccard distance makes the algorithm an interesting candidate for locality-sensitive hashing. Furthermore, the method described in \[sec:cardinality\_estimation\_set\_intersections\] would allow a more detailed estimation of the Jaccard distance by using the estimates for intersection and union sizes. This could be used for additional more precise filtering when searching for similar items. Since sketches can be efficiently constructed, because only a single hash function evaluation is needed for each item, the preprocessing step would be very fast. In contrast, preprocessing is very costly for minwise hashing, because of the many required permutations [@Li2011].

Conclusion {#sec:conclusion}
==========

Based on the Poisson approximation we have presented two new methods to estimate the number of distinct elements from sketches. Unlike previous approaches that use patchworks of different estimators or rely on empirically determined data, the presented ones are inherently unbiased over the full cardinality range. The first extends the original estimator by theoretically derived correction terms in order to obtain an improved estimator that can be straightforwardly translated into a fast algorithm. Due to its simplicity we believe that it has the potential to become the standard textbook cardinality estimator for sketches. The second estimation approach is based on the method and can also be applied to two sketches which significantly improves cardinality estimates for corresponding intersections, relative complements, and unions compared to the conventional approach using the inclusion-exclusion principle.
