Sample complexity of the distinct elements problem
Yihong Wu, Pengkun Yang
Keywords
sampling large population, nonparametric statistics, discrete polynomial approximation, orthogonal polynomials, Vandermonde matrix, minimaxity
AMS 2010 subject classifications
Primary: 62G05; secondary: 62C20, 62D05, 41A05, 41A10
The Distinct Elements problem
The Distinct Elements problem [CCMN00] refers to the following question:
Given balls randomly drawn from an urn containing colored balls, how to estimate the total number of distinct colors in the urn?
Originating from ecology, numismatics, and linguistics, this problem is also known as the species problem in the statistics literature [Lo92, BF93]. Apart from the theoretical interests, it has a wide array of applications in various fields, such as estimating the number of species in a population of animals [FCW43, Goo53], the number of dies used to mint an ancient coinage [Est86], and the vocabulary size of an author [ET76]. In computer science, this problem frequently arises in large-scale databases, network monitoring, and data mining [RRSS09, BYJK+02, CCMN00], where the objective is to estimate the types of database entries or IP addresses from limited observations, since it is typically impossible to have full access to the entire database or keep track of all the network traffic. The key challenge in the Distinct Elements problem is the following: given a small set of samples where most of the colors are not observed, how to accurately extrapolate the number of unseens?
The fundamental limit of the Distinct Elements problem is characterized by the sample complexity, i.e., the smallest sample size needed to estimate the number of distinct colors with a prescribed accuracy and confidence level. A formal definition is the following:
The main results of this paper provide bounds and constant-factor approximations of the sample complexity in various regimes summarized in Table 1, as well as computationally efficient algorithms. Below we highlight a few important conclusions drawn from Table 1:
From the result for in Table 1, we conclude that the sample complexity is sublinear in if and only if , which also holds for sampling without replacement. To estimate within a constant fraction of balls for any small constant , the sample complexity is , which coincides with the general support size estimation problem [VV11a, WY15] (see Section 1.2 for a detailed comparison). However, in other regimes we can achieve better performance by exploiting the discrete nature of the Distinct Elements problem.
The transition from linear to superlinear sample complexity occurs near . Although the exact sample complexity near is not completely resolved in the current paper, the lower bound and upper bound in Table 1 differ by a factor of at most . In particular, the estimator via interpolation can achieve with samples, and achieving a precision of requires strictly superlinear sample size.
To establish the sample complexity, our lower bounds are obtained under zero-one loss and our upper bounds are under the (stronger) quadratic loss. Hence we also obtain the following characterization of the minimax mean squared error (MSE) of the Distinct Elements problem:
where denotes an estimator using samples with replacements and is the number of distinct colors in a -ball urn.
2 Related work
The Distinct Elements problem is equivalent to estimating the number of species (or classes) in a finite population, which has been extensively studied in the statistics (see surveys [BF93, GS04]) and the numismatics literature (see survey [Est86]). Motivated by various practical applications, a number of statistical models have been introduced for this problem, the most popular four being (cf. [BF93, Figure 1]):
The multinomial model: samples are drawn uniformly at random with replacement;
The hypergeometric model: samples are drawn uniformly at random without replacement;
The Bernoulli model: each individual is observed independently with some fixed probability, and thus the total number of samples is a binomial random variable;
The Poisson model: the number of observed samples in each class is independent and Poisson distributed, and thus the total sample size is also a Poisson random variable.
These models are closely related: conditioned on the sample size, the Bernoulli model coincides with the hypergeometric one, and Poisson model coincides with the multinomial one; furthermore, hypergeometric model can simulate multinomial one and is hence more informative. The multinomial model is adopted as the main focus of this paper and the sample complexity in Definition 1 refers to the number of samples with replacement. In the undersampling regime where the sample size is significantly smaller than the population size, all four models are approximately equivalent. See Appendix A for a rigorous justification and detailed comparisons.
Under these models various estimators have been proposed such as unbiased estimators [Goo49], Bayesian estimators [Hil79], variants of Good-Turing estimators [CL92], etc. None of these methodologies, however, have a provable worst-case guarantee. Finally, we mention a closely related problem of estimating the number of connected components in a graph based on sampled induced subgraphs. In the special case where the underlying graph consists of disjoint cliques, the problem is exactly equivalent to the Distinct Elements problem [Fra78].
Computer science literature
The interests in the Distinct Elements problem also arise in the database literature, where various intuitive estimators [HOT88, NS90] have been proposed under simplifying assumptions such as uniformity, and few performance guarantees are available. More recent work in [CCMN00, BYKS01] obtained the optimal sample complexity under the multiplicative error criterion, where the minimum sample size to estimate the number of distinct elements within a factor of is shown to be . For this task, it turns out the least favorable scenario is to distinguish an urn with unitary color from one with almost unitary color, the impossibility of which implies large multiplicative error. However, the optimal estimator performs poorly compared with others on an urn with many distinct colors [CCMN00], the case where most estimators enjoy small multiplicative error. In view of the limitation of multiplicative error, additive error is later considered by [RRSS09, Val11]. To achieve an additive error of for a constant , the result in [CCMN00] only implies an sample complexity lower bound, whereas a much stronger lower bound scales like obtained in [RRSS09], which is almost linear. Determining the optimal sample complexity under additive error is the focus of the present paper.
The Distinct Elements problem can be viewed as a special case of the Support Size problem, where the goal is to estimate the cardinality of the support of an unknown discrete distribution, whose nonzero probabilities are at least , based on independent samples. Improving previous results in [VV11a], the optimal sample complexity has been recently determined in [WY15] to be
Samples drawn from a -ball urn with replacement can be viewed as i.i.d. samples from a distribution supported on the set . From this perspective, any support size estimator, as well as its performance guarantee, is applicable to the Distinct Elements problem.
We briefly describe and compare the strategy to construct estimators in [WY15] and the current paper. Both are based on the idea of polynomial approximation, a powerful tool to circumvent the nonexistence of unbiased estimators [LNS99]. The key is to approximate the function to be estimated by a polynomial, whose degree is chosen to balance the approximation error (bias) and the estimation error (variance). The worst-case performance guarantee for the Support Size problem in [WY15] is governed by the uniform approximation error over an interval where the probabilities may reside. In contrast, in the Distinct Elements problem, samples are generated from a distribution supported on a discrete set of values. Uniform approximation over a discrete subset leads to smaller approximation error and, in turn, improved sample complexity. It turns out that samples are sufficient to achieve an additive error of that satisfies , which strictly improves the sample complexity (1) for the Support Size problem, thanks to the discrete structure of the Distinct Elements problem.
The Distinct Elements problem considered here is not to be confused with the formulation in the streaming literature, where the goal is to approximate the number of distinct elements in the observations with low space complexity, see, e.g., [FFGM07, KNW10]. There, the proposed algorithms aim to optimize the memory consumption, but still require a full pass of every ball in the urn. This is different from the setting in the current paper, where only random samples drawn from the urn are available.
3 Organization
The paper is organized as follows: In Section 2 we describe a unified approach to construct estimators via discrete polynomial approximation, whose bias is analyzed in Section 2.2 and variance is upper bounded in Sections 2.3 and 2.4 separately. In Section 3 we obtain lower bounds on the sample complexity in Table 1 which establish the optimality of the proposed estimators. Section 4 explains how sample complexity bounds summarized in Table 1 follow from various results in Sections 2 and 3. Connections between the four sampling model mentioned in Section 1.2 are detailed in Appendix A. Proofs of auxiliary results are deferred to Appendix B and Appendix C.
4 Notations
Linear estimators via discrete polynomial approximation
The naïve estimator, “what you see is what you get,” is simply the number of observed distinct colors, which can be expressed in terms of fingerprints as
This is typically an underestimator because . In turn, our estimator is
where the coefficients ’s are to be specified. Since the fingerprints are dependent (for example, they sum up to ), (3) serves as a linear predictor of in terms of the observed fingerprints. Equivalently, in terms of histograms, the estimator has the following decomposable form:
and is a (formal) power series with . The right-hand side of (5) can be made zero by choosing to be, e.g., the Lagrange interpolating polynomial that satisfies and for , namely, ; however, this strategy results in a high-degree polynomial with large coefficients, which, in turn, leads to a large variance of the estimator.
To reduce the variance of our estimator, we only use the first fingerprints in (3) by setting for all , where is chosen to be proportional to . This restricts the polynomial degree in (5) to at most and, while possibly incurring bias, reduces the variance. A further reason for only using the first few fingerprints is that higher-order fingerprints are almost uncorrelated with the number of unseens . For instance, if red balls are observed for times, the only information this reveals is that approximately half of the urn are red. In fact, the correlation between and decays exponentially (see Appendix B for a proof). Therefore for , offer little predictive power about . Moreover, if a color is observed at most times, say, , this implies that, with high probability, , where , thanks to the concentration of Poisson random variables. Therefore, effectively we only need to consider those colors that appear in the urn for at most times, i.e., , for which the bias is at most
where , , and
In view of the bias analysis in (6), we have
Proposition 1 suggests that the coefficients of the linear estimator can be chosen by solving the following linear programming (LP):
which is an upper bound of (13), and is in fact within an factor since and . In the remainder of this section, we consider two separate cases:
(): In this case, the linear system in (14) is overdetermined and the minimum is non-zero. Surprisingly, as shown in Section 2.2, the exact optimal value can be found in closed form using discrete orthogonal polynomials. The coefficients of the solution can be bounded using the minimum singular value of the matrix , which is analyzed in Section 2.3 .
(): In this case, the linear system is underdetermined and the minimum in (14) is zero. To bound the variance, it turns out that the coefficients bound obtained from the minimum singular value is not precise enough in this regime. Instead, we express the coefficients in terms of Lagrange interpolating polynomials and use Stirling numbers to obtain sharp variance bounds. This analysis in carried out in Section 2.4.
We finish this subsection with two remarks:
The optimal estimator for the Support Size problem in [WY15] has the same linear form as (2); however, since the probabilities can take any values in an interval, the coefficients are found to be the solution of the continuous polynomial approximation problem
where the infimum is taken over all degree- polynomials such that , achieved by the (appropriately shifted and scaled) Chebyshev polynomial [Tim63]. In contrast, in Section 2.2 we show that the discrete version of (15), which is equivalent to the LP (13), satisfies
provided . The difference between (15) and (16) explains why the sample complexity (1) for the Support Size problem has an extra log factor compared to that of the Distinct Elements problem in Table 1. When the sample size is large enough, interpolation is used in lieu of approximation. See Fig. 1 for an illustration.
The time complexity of the estimator (2) consists of: (a) Computing histograms and fingerprints of samples: ; (b) Computing the coefficients by solving the least square problem in (6): ; (c) Evaluating the linear combination (2): . As shown in Table 1, for an accurate estimation the sample complexity is , which implies and . Therefore, the overall time complexity is .
Define the following inner product between functions and :
and the induced norm . The least square problem (17) can be equivalently formulated as
This can be analyzed using the orthogonal polynomials under the inner product (18), which we describe next.
Recall the discrete Chebyshev polynomial [Sze75, Sec. 2.8]: for ,
and denotes the -th order forward difference. The polynomials are orthogonal with respect to the counting measure over the discrete set ; in particular, we have (cf. [Sze75, Sec. 2.8.2, 2.8.3]):
By appropriately shifting and scaling the set of polynomials , we define an orthonormal basis for the set of polynomials of degree at most under the inner product (18) by
Since constitute a basis for polynomials of degree at most , the least square problem (19) can be equivalently formulated as
where , , and denotes vector inner product. Thus, the optimal value is clearly , achieved by .
From (21) we have . By the formula of in (20), we obtain
In view of the definition of in (22), we have
where the last equality follows from induction since
The second equality in (17) is a direct consequence of Stirling’s approximation. If , then
If , denoting and applying when , we have
where the last step follows from when . In the exponent of (24), the term dominates when . Applying (23) and (24) to the exact solution (17) yields the desired approximation. ∎
3 Minimum singular values of real rectangle Vandermonde matrices
In Proposition 1 the variance of our estimator is bounded by the magnitude of coefficients , which is related to the polynomial coefficients by (7). A classical result from approximation theory is that if a polynomial is bounded over a compact interval, its coefficients are at most exponential in the degree [Tim63, Theorem 2.9.11]: for any degree- polynomial ,
which is tight when is the Chebyshev polynomial. This fact has been applied in statistical contexts to control the variance of estimators obtained from best polynomial approximation [CL11, WY16, WY15, JVHW15]. In contrast, for the Distinct Elements problem, the polynomial is only known to be bounded over the discretized interval. Nevertheless, we show that the bound (25) continues to hold as long as the discretization level exceeds the degree:
provided that (see Remark 3 after Lemma 2). Clearly, (26) implies (25) by sending . If , a coefficient bound like (26) is impossible, because one can add to an arbitrary degree- interpolating polynomial that evaluates to zero at all points.
where denotes the smallest singular value of . Let
which is a well-studied example of ill-conditioned matrices in the numerical analysis literature. In particular, it is known that the condition number of the Hilbert matrix is [Tod54] and the operator norm is , and thus the minimum singular value is exponentially small in the degree. Therefore we expect the discrete moment matrix to behave similarly to the Hilbert matrix when is large enough. Interestingly, we show that this is indeed the case as soon as exceeds (otherwise the minimum singular value is zero).
The inequality (26) follows from Lemma 2 since the coefficient vector satisfies .
The extreme singular values of square Vandermonde matrices have been extensively studied (c.f. [Gau90, Bec00] and the references therein). For rectangular Vandermonde matrices, the focus was mainly with nodes on the unit circle in the complex domain [CGR90, Fer99, Moi15] with applications in signal processing. In contrast, Lemma 2 is on rectangular Vandermonde matrices with real nodes. The result on integers nodes in [EPS01] turns out to be too crude for the purpose of this paper.
For a given , the orthonormal basis in (22) is proportional to the discrete Chebyshev polynomials . The classical asymptotic result for the discrete Chebyshev polynomials shows that [Sze75, (2.8.6)]
where is the Legendre polynomial of degree . This gives the intuition that for real-valued . We have the following non-asymptotic upper bound (proved in Appendix C) for over the complex plane:
Applying (32) on the definition of in (22), for any and any , we have
The right-hand side is increasing with . Therefore,
where in the last inequality we used and . ∎
Recall the connection between and in (7). For , we have . Therefore,
Applying (34) and (35) to Proposition 1, we obtain
Then the desired (33) holds as long as is sufficiently large and is sufficiently small. ∎
4 Lagrange interpolating polynomials and Stirling numbers
When we sample at least a constant faction of the urn, i.e., , we can afford to choose and in (8) so that and is an invertible matrix. We choose the coefficient which is equivalent to applying Lagrange interpolating polynomial and achieves exact zero bias. To control the variance, we can follow the approach in Section 2.3 by using the bound on minimum singular value of the matrix , which implies that the coefficients are and yields a coarse upper bound on the sample complexity. As previously announced in Table 1, this bound can be improved to by a more careful analysis of the Lagrange interpolating polynomial coefficients expressed in terms of the Stirling numbers, which we introduce next.
The Stirling numbers of the first kind are defined as the coefficients of the falling factorial where
Compared to the coefficients expressed by the Lagrange interpolating polynomial:
we obtain a formula for the coefficients in terms of the Stirling numbers:
Consequently, the coefficients of our estimator are given by
The precise asymptotics the Stirling number is rather complicated. In particular, the asymptotic formula of as for fixed is given by [Jor47] and the uniform asymptotics over all is obtained in [MW58] and [Tem93]. The following lemma (proved in Appendix C) is a coarse non-asymptotic version, which suffices for the purpose of constant-factor approximations of the sample complexity.
We construct as in Proposition 1 using the coefficients in (36) to achieve zero bias. The variance upper bound by the coefficients is a direct consequence of the upper bound of Stirling numbers in Lemma 4. Then we obtain the following mean squared error (MSE):
Assume the Poisson sampling model. If for some sufficiently large constant , then
In Proposition 1, fix and so that . Our goal is to show an upper bound of
Here the coefficients are obtained from (36) and, in view of (37), satisfy:
for some universal constant . We consider three cases separately:
In view of (39) and , we have . Then,
as long as and thus . Therefore,
Case II: ηkloglogk≤n≤βklogk𝜂𝑘𝑘𝑛𝛽𝑘𝑘\eta k\log\log k\leq n\leq\sqrt{\beta}k\sqrt{\log k}.
Case III: ηk≤n≤ηkloglogk𝜂𝑘𝑛𝜂𝑘𝑘\eta k\leq n\leq\eta k\log\log k.
We apply the upper bound of expectation by the maximum:
Since , the right-hand side of (39) is decreasing with when , so it suffices to consider . Denoting and , in view of (39), we have , which attains maximum at satisfying . Then,
where the last inequality is because of . Therefore,
Applying the upper bounds in (41), (43) and (44) to Proposition 1 concludes the proof. ∎
It is impossible to bridge the gap near in Table 1 using the technology of interpolating polynomials that aims at zero bias, since its worst-case variance is at least when . To see this, note that the variance term given by (12) is
Optimality of the sample complexity
In this section we develop lower bounds of the sample complexity which certify the optimality of estimators constructed in Section 2. We first give a brief overview of the lower bound in [CCMN00, Theorem 1], which gives the optimal sample complexity under the multiplicative error criterion. The lower bound argument boils down to considering two hypothesis: in the null hypothesis, the urn consists of only one color; in the alternative, the urn contains distinct colors, where balls share the same color as in the null hypothesis, and all other balls have distinct colors. These two scenarios are distinguished if and only if a second color appears in the samples, which typically requires samples. This lower bound is optimal for estimating within a multiplicative factor of , which, however, is too loose for additive error .
In contrast, instead of testing whether the urn is monochromatic, our first lower bound is given by testing whether the urn is maximally colorful, that is, containing distinct colors. The alternative contains colors, and the numbers of balls of two different colors differ by at most one. In other words, the null hypothesis is the uniform distribution on and the alternative is close to uniform distribution with smaller support size. The sample complexity, which is shown in Theorem 3, gives the lower bound in Table 1 for .
Consider the following two hypotheses: The null hypothesis is an urn consisting of distinct colors; The alternative consists of distinct colors, and each color appears either or times. In terms of distributions, is the uniform distribution ; is the closest perturbation from the uniform distribution: randomly pick disjoint sets of indices with cardinality and , where and satisfy
Conditional on , the distribution is given by
Put the uniform prior on the alternative. Denote the marginal distributions of the samples under and by and , respectively. Since the distinct colors in and are separated by , to show that the sample complexity , it suffices to show that no test can distinguish and reliably using samples. A further sufficient condition is a bounded divergence [Tsy09]
The remainder of this proof is devoted to upper bounds of the divergence.
Since and , we have
where is an independent copy of . By the definition of and ,
where , , , and are centered random variables. Applying and Cauchy-Schwarz inequality, we obtain
where the last inequality follows from the fact that is increasing when . Other terms in (49) are bounded analogously and we have
If , the upper bound (51) implies that since the -divergence is finite with samples, using the inequality that for ; if , the lower bound is trivial since .
with for and for . Therefore,
The upper bound (52) yields the sample complexity . ∎
Suppose that we take i.i.d. samples from , which form a -ball urn consisting of distinct colors. By the union bound,
Combining with the sample complexity of the Support Size problem in (1), Lemma 5 leads to the following lower bound for the Distinct Elements problem:
Fix a sufficiently small constant . For any ,
The same lower bound holds for sampling without replacement.
Proof of results in Table 1
Below we explain how the sample complexity bounds summarized in Table 1 are obtained from various results in Section 2 and Section 3:
The upper bounds are obtained from the worst-case MSE in Section 2 and the Markov inequality. In particular, the case of follows from the second and the third upper bounds of Theorem 2; the case of follows from the first upper bound of Theorem 2; the case of follows from Theorem 1. By monotonicity, we have the upper bound when , the upper bound when , and the upper bound when .
The lower bound for follows from Theorem 3; the lower bound for follows from Theorem 4. These further implies the lower bound for by monotonicity.
Appendix A Connections between various sampling models
As mentioned in Section 1.2, four popular sampling models have been introduced in the statistics literature: the multinomial model, the hypergeometric model, the Bernoulli model, and the Poisson model. The connections between those models are explained in details in this section, as well as relations between the respective sample complexities.
The connections between different models are illustrated in Fig. 2. Under the Poisson model, the sample size is a Poisson random variable; conditioned on the sample size, the samples are i.i.d. which is identical to the multinomial model. The same relation holds as the Bernoulli model to the hypergeometric model. Given samples uniformly drawn from a -ball urn without replacement (hypergeometric model), we can simulate drawn with replacement (multinomial model) as follows: for each , let
In view of the connections in Fig. 2, any estimator constructed for one specific model can be adapted to another. The adaptation from multinomial to hypergeometric model is provided by the simulation in (53), and the other direction is given by Lemma 5 (without modifying the estimator). The following result provides a recipe for going between fixed and randomized sample size:
Denote the samples by . Following [RRSS09, Lemma 5.3(a)], define as
;
;
.
;
.
Appendix B Correlation decay between fingerprints
The correlation coefficient between and follows immediately:
where . Note that as . Therefore, for any ,
where is uniform as . Taking derivative, the function on is increasing if and only if , and the maximum is attained at . Therefore, applying ,
Appendix C Proof of auxiliary lemmas
Taking -th derivative of , we obtain
Then the desired (32) follows from (57). ∎
The following uniform asymptotic expansions of the Stirling numbers of the first kind was obtained in [CRT00, Theorem 2]:
where is Euler’s constant, is the unique positive solution to with , , and all terms are uniform in . In the following we consider each range separately and prove the non-asymptotic approximation in (37).
Case I. For , Stirling’s approximation gives
Case III. For , note that , and thus . By [MW58, Lemma 4.1], in this range. Hence,
where is the solution to . Bounding the sum by integrals, we have
If , then , and hence
In view of (58), we have , which is exactly (37) when . If , then , and
Combining (58) yields that , which coincides with (37) since is this range. ∎
Acknowledgment
This research has been supported in part by the National Science Foundation under the grant agreement IIS-14-47879 and CCF-15-27105 and an NSF CAREER award CCF-1651588. The authors thank Greg Valiant for helpful discussion on Lemma 5.14 in his thesis [Val12]. We thank the anonymous referees for constructive comments which have helped to improve the presentation of the paper.