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 nn balls randomly drawn from an urn containing kk 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 k0.5+δ≤Δ≤ckk^{0.5+\delta}\leq\Delta\leq ck in Table 1, we conclude that the sample complexity is sublinear in kk if and only if Δ=k1−o(1)\Delta=k^{1-o(1)}, which also holds for sampling without replacement. To estimate within a constant fraction of balls Δ=ck\Delta=ck for any small constant cc, the sample complexity is Θ(klog⁡k)\Theta(\frac{k}{\log k}), 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 Δ=k\Delta=\sqrt{k}. Although the exact sample complexity near Δ=k\Delta=\sqrt{k} is not completely resolved in the current paper, the lower bound and upper bound in Table 1 differ by a factor of at most log⁡log⁡k\log\log k. In particular, the estimator via interpolation can achieve Δ=k\Delta=\sqrt{k} with n=O(klog⁡log⁡k)n=O(k\log\log k) samples, and achieving a precision of Δ≤k0.5−o(1)\Delta\leq k^{0.5-o(1)} 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 C^\hat{C} denotes an estimator using nn samples with replacements and CC is the number of distinct colors in a kk-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: nn samples are drawn uniformly at random with replacement;

The hypergeometric model: nn 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 α\alpha is shown to be Θ(k/α2)\Theta(k/\alpha^{2}). 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 ckck for a constant c∈(0,12)c\in(0,\frac{1}{2}), the result in [CCMN00] only implies an Ω(1/c)\Omega(1/c) sample complexity lower bound, whereas a much stronger lower bound scales like k1−O(log⁡log⁡klog⁡k)k^{1-O(\sqrt{\frac{\log\log k}{\log k}})} 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 1k\frac{1}{k}, 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 kk-ball urn with replacement can be viewed as i.i.d. samples from a distribution supported on the set {1k,2k,…,kk}\{\frac{1}{k},\frac{2}{k},\dots,\frac{k}{k}\}. 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 O(klog⁡klog⁡kΔ)O(\frac{k}{\log k}\log\frac{k}{\Delta}) samples are sufficient to achieve an additive error of Δ\Delta that satisfies k0.5+O(1)≤Δ≤O(k)k^{0.5+O(1)}\leq\Delta\leq O(k), 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 C=C^seen+UC=\hat{C}_{\rm seen}+U. In turn, our estimator is

where the coefficients uju_{j}’s are to be specified. Since the fingerprints Φ0,Φ1,…\Phi_{0},\Phi_{1},\dots are dependent (for example, they sum up to CC), (3) serves as a linear predictor of U=Φ0U=\Phi_{0} in terms of the observed fingerprints. Equivalently, in terms of histograms, the estimator has the following decomposable form:

and ϕ(a)≜∑j≥1ajuj(n/k)jj!\phi(a)\triangleq\sum_{j\geq 1}a^{j}\frac{u_{j}(n/k)^{j}}{j!} is a (formal) power series with ϕ(0)=0\phi(0)=0. The right-hand side of (5) can be made zero by choosing ϕ\phi to be, e.g., the Lagrange interpolating polynomial that satisfies ϕ(0)=−1\phi(0)=-1 and ϕ(i)=0\phi(i)=0 for i∈[k]i\in[k], namely, ϕ(a)=(−1)k+1k!∏i=1k(a−i)\phi(a)=\frac{(-1)^{k+1}}{k!}\prod_{i=1}^{k}(a-i); however, this strategy results in a high-degree polynomial ϕ\phi 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 LL fingerprints in (3) by setting uj=0u_{j}=0 for all j>Lj>L, where LL is chosen to be proportional to log⁡k\log k. This restricts the polynomial degree in (5) to at most LL 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 Φ0\Phi_{0}. For instance, if red balls are observed for n/2n/2 times, the only information this reveals is that approximately half of the urn are red. In fact, the correlation between Φ0\Phi_{0} and Φj\Phi_{j} decays exponentially (see Appendix B for a proof). Therefore for L=Θ(log⁡k)L=\Theta(\log k), {Φj}j>L\{\Phi_{j}\}_{j>L} offer little predictive power about Φ0\Phi_{0}. Moreover, if a color is observed at most LL times, say, Ni≤LN_{i}\leq L, this implies that, with high probability, ki≤Mk_{i}\leq M, where M=O(kL/n)M=O(kL/n), 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 MM times, i.e., ki∈[M]k_{i}\in[M], for which the bias is at most

where p(x)≜ϕ(Mx)=∑j=1Lwjxjp(x)\triangleq\phi(Mx)=\sum_{j=1}^{L}w_{j}x^{j}, w=(w1,…,wL)⊤w=(w_{1},\dots,w_{L})^{\top}, 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 O(log⁡k)O(\log k) factor since M=O(klog⁡k/n)M=O(k\log k/n) and n=Ω(k/log⁡k)n=\Omega(k/\log k). In the remainder of this section, we consider two separate cases:

M>LM>L (n≲kn\lesssim k): 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 BB, which is analyzed in Section 2.3 .

M≤LM\leq L (n≳kn\gtrsim k): 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-LL polynomials such that p(0)=0p(0)=0, 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 L<ML<M. 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 nn 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 NiN_{i} and fingerprints Φj\Phi_{j} of nn samples: O(n)O(n); (b) Computing the coefficients ww by solving the least square problem in (6): O(L2(M+L))O(L^{2}(M+L)); (c) Evaluating the linear combination (2): O(n∧k)O(n\wedge k). As shown in Table 1, for an accurate estimation the sample complexity is n=Ω(klog⁡k)n=\Omega(\frac{k}{\log k}), which implies L=O(log⁡k)L=O(\log k) and M=O(log⁡2k)M=O(\log^{2}k). Therefore, the overall time complexity is O(n+log⁡4k)=O(n)O(n+\log^{4}k)=O(n).

Define the following inner product between functions ff and gg:

and the induced norm ∥f∥≜⟨f,f⟩\left\|{f}\right\|\triangleq\sqrt{\langle{f,f}\rangle}. 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 x=0,1,…,M−1x=0,1,\dots,M-1,

and Δm\Delta^{m} denotes the mm-th order forward difference. The polynomials {t0,…,tM−1}\{t_{0},\ldots,t_{M-1}\} are orthogonal with respect to the counting measure over the discrete set {0,1,…,M−1}\left\{0,1,\dots,M-1\right\}; in particular, we have (cf. [Sze75, Sec. 2.8.2, 2.8.3]):

By appropriately shifting and scaling the set of polynomials tmt_{m}, we define an orthonormal basis for the set of polynomials of degree at most L≤M−1L\leq M-1 under the inner product (18) by

Since {ϕm}m=0L\{\phi_{m}\}_{m=0}^{L} constitute a basis for polynomials of degree at most LL, the least square problem (19) can be equivalently formulated as

where ϕ(0)≜(ϕ0(0),…,ϕL(0))\phi(0)\triangleq(\phi_{0}(0),\dots,\phi_{L}(0)), a=(a0,…,aL)a=(a_{0},\ldots,a_{L}), and ⟨⋅,⋅⟩\left\langle\cdot,\cdot\right\rangle denotes vector inner product. Thus, the optimal value is clearly 1∥ϕ(0)∥2\frac{1}{\left\|{\phi(0)}\right\|_{2}}, achieved by a∗=−ϕ(0)∥ϕ(0)∥22a^{*}=-\frac{\phi(0)}{\left\|{\phi(0)}\right\|_{2}^{2}}.

From (21) we have pm(0)=pm(1)=⋯=pm(m−1)=0p_{m}(0)=p_{m}(1)=\dots=p_{m}(m-1)=0. By the formula of tmt_{m} in (20), we obtain

In view of the definition of ϕm\phi_{m} 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 M=L+1M=L+1, then

If M≥L+2M\geq L+2, denoting x=L+1Mx=\frac{L+1}{M} and applying n!=2πn(ne)n(1+Θ(1n))n!=\sqrt{2\pi n}(\frac{n}{e})^{n}(1+\Theta(\frac{1}{n})) when n≥1n\geq 1, we have

where the last step follows from (1+x)log⁡(1+x)+(1−x)log⁡(1−x)=Θ(x2)(1+x)\log(1+x)+(1-x)\log(1-x)=\Theta(x^{2}) when 0≤x≤10\leq x\leq 1. In the exponent of (24), the term Θ(Mx2)\Theta(Mx^{2}) dominates when M≥L+2M\geq L+2. 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 uu, which is related to the polynomial coefficients ww 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-LL polynomial p(x)=∑i=0Lwixip(x)=\sum_{i=0}^{L}w_{i}x^{i},

which is tight when pp 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 M≥L+1M\geq L+1 (see Remark 3 after Lemma 2). Clearly, (26) implies (25) by sending M→∞M\to\infty. If M≤LM\leq L, a coefficient bound like (26) is impossible, because one can add to pp an arbitrary degree-LL interpolating polynomial that evaluates to zero at all MM points.

where σmin⁡(B)\sigma_{\min}(B) denotes the smallest singular value of BB. 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 L×LL\times L Hilbert matrix is O((1+2)4LL)O(\frac{(1+\sqrt{2})^{4L}}{\sqrt{L}}) [Tod54] and the operator norm is Θ(1)\Theta(1), and thus the minimum singular value is exponentially small in the degree. Therefore we expect the discrete moment matrix 1MBˉ⊤Bˉ\frac{1}{M}\bar{B}^{\top}\bar{B} to behave similarly to the Hilbert matrix when MM is large enough. Interestingly, we show that this is indeed the case as soon as MM exceeds LL (otherwise the minimum singular value is zero).

The inequality (26) follows from Lemma 2 since the coefficient vector w=(w0,…,wL)w=(w_{0},\ldots,w_{L}) satisfies ∥w∥∞≤∥w∥2≤1σmin⁡(Bˉ)∥Bˉw∥2≤Mσmin⁡(Bˉ)∥Bˉw∥∞\|w\|_{\infty}\leq\|w\|_{2}\leq\frac{1}{\sigma_{\min}(\bar{B})}\|\bar{B}w\|_{2}\leq\frac{\sqrt{M}}{\sigma_{\min}(\bar{B})}\|\bar{B}w\|_{\infty}.

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 MM, the orthonormal basis ϕm(x)\phi_{m}(x) in (22) is proportional to the discrete Chebyshev polynomials tm(Mx−1)t_{m}(Mx-1). The classical asymptotic result for the discrete Chebyshev polynomials shows that [Sze75, (2.8.6)]

where PmP_{m} is the Legendre polynomial of degree mm. This gives the intuition that tm(x)≈Mmt_{m}(x)\approx M^{m} for real-valued x∈[0,M]x\in[0,M]. We have the following non-asymptotic upper bound (proved in Appendix C) for tmt_{m} over the complex plane:

Applying (32) on the definition of ϕm\phi_{m} in (22), for any ∣z∣=1|z|=1 and any M≥L+1M\geq L+1, we have

The right-hand side is increasing with mm. Therefore,

where in the last inequality we used (nk)≥(nk)k\binom{n}{k}\geq(\frac{n}{k})^{k} and n!≥(ne)nn!\geq(\frac{n}{e})^{n}. ∎

Recall the connection between uju_{j} and wjw_{j} in (7). For 1≤j≤L<βlog⁡k1\leq j\leq L<\beta\log k, we have uj=wjj!(βlog⁡k)j≤wjβlog⁡ku_{j}=w_{j}\frac{j!}{(\beta\log k)^{j}}\leq\frac{w_{j}}{\beta\log k}. Therefore,

Applying (34) and (35) to Proposition 1, we obtain

Then the desired (33) holds as long as β\beta is sufficiently large and α\alpha is sufficiently small. ∎

4 Lagrange interpolating polynomials and Stirling numbers

When we sample at least a constant faction of the urn, i.e., n=Ω(k)n=\Omega(k), we can afford to choose α\alpha and β\beta in (8) so that L=ML=M and BB is an invertible matrix. We choose the coefficient w=B−11w=B^{-1}\mathbf{1} 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 BB, which implies that the coefficients are exp⁡(O(L))\exp(O(L)) and yields a coarse upper bound O(klog⁡k1∨log⁡Δ2k)O(k\frac{\log k}{1\vee\log\frac{\Delta^{2}}{k}}) on the sample complexity. As previously announced in Table 1, this bound can be improved to O(klog⁡log⁡k1∨log⁡Δ2k)O(k\log\frac{\log k}{1\vee\log\frac{\Delta^{2}}{k}}) 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 (x)n(x)_{n} where

Compared to the coefficients ww expressed by the Lagrange interpolating polynomial:

we obtain a formula for the coefficients ww in terms of the Stirling numbers:

Consequently, the coefficients of our estimator uju_{j} are given by

The precise asymptotics the Stirling number is rather complicated. In particular, the asymptotic formula of s(n,m)s(n,m) as n→∞n\rightarrow\infty for fixed mm is given by [Jor47] and the uniform asymptotics over all mm 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 C^\hat{C} as in Proposition 1 using the coefficients uju_{j} in (36) to achieve zero bias. The variance upper bound by the coefficients uu 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 n>ηkn>\eta k for some sufficiently large constant η\eta, then

In Proposition 1, fix β=3.5\beta=3.5 and α=βkn\alpha=\frac{\beta k}{n} so that L=ML=M. Our goal is to show an upper bound of

Here the coefficients uju_{j} are obtained from (36) and, in view of (37), satisfy:

for some universal constant η\eta. We consider three cases separately:

In view of (39) and j≥1j\geq 1, we have ∣uj∣≤(Θ(k/n)log⁡M)j|u_{j}|\leq(\Theta(k/n)\log M)^{j}. Then,

as long as n≳klog⁡log⁡kn\gtrsim k\log\log k and thus klog⁡Mn≲1\frac{k\log M}{n}\lesssim 1. Therefore,

Case II: η​k​log⁡log⁡k≤n≤β​k​log⁡k𝜂𝑘𝑘𝑛𝛽𝑘𝑘\eta k\log\log k\leq n\leq\sqrt{\beta}k\sqrt{\log k}.

Case III: η​k≤n≤η​k​log⁡log⁡k𝜂𝑘𝑛𝜂𝑘𝑘\eta k\leq n\leq\eta k\log\log k.

We apply the upper bound of expectation by the maximum:

Since ηkn≤1\frac{\eta k}{n}\leq 1, the right-hand side of (39) is decreasing with jj when j≥M/ej\geq M/e, so it suffices to consider j≤M/ej\leq M/e. Denoting x=log⁡Mjx=\log\frac{M}{j} and τ=Θ(kn)\tau=\Theta(\frac{k}{n}), in view of (39), we have ∣uj∣≤exp⁡(Me−xlog⁡(τx))|u_{j}|\leq\exp(Me^{-x}\log(\tau x)), which attains maximum at x∗x^{*} satisfying e1/x∗x∗=τ\frac{e^{1/x^{*}}}{x^{*}}=\tau. Then,

where the last inequality is because of τ>1x∗\tau>\frac{1}{x^{*}}. Therefore,

Applying the upper bounds in (41), (43) and (44) to Proposition 1 concludes the proof. ∎

It is impossible to bridge the gap near Δ=k\Delta=\sqrt{k} in Table 1 using the technology of interpolating polynomials that aims at zero bias, since its worst-case variance is at least k1+Ω(1)k^{1+\Omega(1)} when n=O(k)n=O(k). 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 2Δ+12\Delta+1 distinct colors, where k−2Δk-2\Delta 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 Ω(k/Δ)\Omega(k/\Delta) samples. This lower bound is optimal for estimating within a multiplicative factor of Δ\sqrt{\Delta}, which, however, is too loose for additive error Δ\Delta.

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 kk distinct colors. The alternative contains k−2Δk-2\Delta 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 [k][k] 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 Δ≤k\Delta\leq\sqrt{k}.

Consider the following two hypotheses: The null hypothesis H0H_{0} is an urn consisting of kk distinct colors; The alternative H1H_{1} consists of k−2Δk-2\Delta distinct colors, and each color appears either b1≜⌊kk−2Δ⌋b_{1}\triangleq\lfloor{\frac{k}{k-2\Delta}}\rfloor or b2≜⌈kk−2Δ⌉b_{2}\triangleq{\lceil{\frac{k}{k-2\Delta}}\rceil} times. In terms of distributions, H0H_{0} is the uniform distribution Q=(1k,…,1k)Q=(\frac{1}{k},\dots,\frac{1}{k}); H1H_{1} is the closest perturbation from the uniform distribution: randomly pick disjoint sets of indices I,J⊆[k]I,J\subseteq[k] with cardinality ∣I∣=c1|I|=c_{1} and ∣J∣=c2|J|=c_{2}, where c1c_{1} and c2c_{2} satisfy

Conditional on θ≜(I,J)\theta\triangleq(I,J), the distribution Pθ=(pθ,1,…,pθ,k)P_{\theta}=(p_{\theta,1},\dots,p_{\theta,k}) is given by

Put the uniform prior on the alternative. Denote the marginal distributions of the nn samples X=(X1,…,Xn)X=(X_{1},\dots,X_{n}) under H0H_{0} and H1H_{1} by QXQ_{X} and PXP_{X}, respectively. Since the distinct colors in H0H_{0} and H1H_{1} are separated by 2Δ2\Delta, to show that the sample complexity n∗(k,Δ)≥nn^{*}(k,\Delta)\geq n, it suffices to show that no test can distinguish H0H_{0} and H1H_{1} reliably using nn samples. A further sufficient condition is a bounded χ2\chi^{2} divergence [Tsy09]

The remainder of this proof is devoted to upper bounds of the χ2\chi^{2} divergence.

Since PX∣θ=Pθ⊗nP_{X|\theta}=P_{\theta}^{\otimes n} and QX=Q⊗nQ_{X}=Q^{\otimes n}, we have

where θ′\theta^{\prime} is an independent copy of θ\theta. By the definition of PθP_{\theta} and QQ,

where A1≜b12k(∣I∩I′∣−c12k)A_{1}\triangleq\frac{b_{1}^{2}}{k}(|I\cap I^{\prime}|-\frac{c_{1}^{2}}{k}), A2≜b22k(∣J∩J′∣−c22k)A_{2}\triangleq\frac{b_{2}^{2}}{k}(|J\cap J^{\prime}|-\frac{c_{2}^{2}}{k}), A3=b1b2k(∣I∩J′∣−c1c2k)A_{3}=\frac{b_{1}b_{2}}{k}(|I\cap J^{\prime}|-\frac{c_{1}c_{2}}{k}), and A4=b1b2k(∣J∩I′∣−c1c2k)A_{4}=\frac{b_{1}b_{2}}{k}(|J\cap I^{\prime}|-\frac{c_{1}c_{2}}{k}) are centered random variables. Applying 1+x≤ex1+x\leq e^{x} and Cauchy-Schwarz inequality, we obtain

where the last inequality follows from the fact that x↦ex−1−xx\mapsto e^{x}-1-x is increasing when x>0x>0. Other terms in (49) are bounded analogously and we have

If k−2Δ≥kk-2\Delta\geq\sqrt{k}, the upper bound (51) implies that n∗(k,Δ)≥Ω(k−2Δk)n^{*}(k,\Delta)\geq\Omega(\frac{k-2\Delta}{\sqrt{k}}) since the χ2\chi^{2}-divergence is finite with O(k−2Δk)O(\frac{k-2\Delta}{\sqrt{k}}) samples, using the inequality that ex−1−x≤x22e^{x}-1-x\leq\frac{x^{2}}{2} for x≥0x\geq 0; if k−2Δ≤kk-2\Delta\leq\sqrt{k}, the lower bound is trivial since k−2Δk≤1\frac{k-2\Delta}{\sqrt{k}}\leq 1.

with t=4nkt=\frac{4n}{k} for i=1,2i=1,2 and t=−4nkt=-\frac{4n}{k} for i=3,4i=3,4. Therefore,

The upper bound (52) yields the sample complexity n∗(k,Δ)≥Ω(karccosh⁡(1+k4Δ2))n^{*}(k,\Delta)\geq\Omega(k\operatorname{arccosh}(1+\frac{k}{4\Delta^{2}})). ∎

Suppose that we take kk i.i.d. samples from P=(p1,p2,… )P=(p_{1},p_{2},\dots), which form a kk-ball urn consisting of CC 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 cc. For any 1≤Δ≤ck1\leq\Delta\leq ck,

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 Δ≤k(log⁡k)−δ\Delta\leq\sqrt{k}(\log k)^{-\delta} follows from the second and the third upper bounds of Theorem 2; the case of k≤Δ≤k0.5+δ\sqrt{k}\leq\Delta\leq k^{0.5+\delta} follows from the first upper bound of Theorem 2; the case of k1−δ≤Δ≤ckk^{1-\delta}\leq\Delta\leq ck follows from Theorem 1. By monotonicity, we have the O(klog⁡log⁡k)O(k\log\log k) upper bound when k(log⁡k)−δ≤Δ≤k\sqrt{k}(\log k)^{-\delta}\leq\Delta\leq\sqrt{k}, the O(klog⁡k)O(\frac{k}{\log k}) upper bound when Δ≥ck\Delta\geq ck, and the O(k)O(k) upper bound when k0.5+δ≤Δ≤k1−δk^{0.5+\delta}\leq\Delta\leq k^{1-\delta}.

The lower bound for Δ≤k\Delta\leq\sqrt{k} follows from Theorem 3; the lower bound for k0.5+δ≤Δ≤ckk^{0.5+\delta}\leq\Delta\leq ck follows from Theorem 4. These further implies the Ω(k)\Omega(k) lower bound for k≤Δ≤k0.5+δ\sqrt{k}\leq\Delta\leq k^{0.5+\delta} 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 (Y1,…,Yn)(Y_{1},\dots,Y_{n}) uniformly drawn from a kk-ball urn without replacement (hypergeometric model), we can simulate (X1,…,Xn)(X_{1},\dots,X_{n}) drawn with replacement (multinomial model) as follows: for each i=1,…,ni=1,\dots,n, 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 X1,…,XNX_{1},\dots,X_{N}. Following [RRSS09, Lemma 5.3(a)], define C^′\hat{C}^{\prime} as

nH∗(k,Δ,δ)≤nM∗(k,Δ,δ)n^{*}_{H}(k,\Delta,\delta)\leq n^{*}_{M}(k,\Delta,\delta);

nP∗(k,Δ,δ)≤n⇒nM∗(k,Δ,δ+(e/4)n)≤2nn_{P}^{*}(k,\Delta,\delta)\leq n\Rightarrow n_{M}^{*}(k,\Delta,\delta+(e/4)^{n})\leq 2n;

nM∗(k,Δ,δ)≤n⇒nP∗(k,Δ,δ+(2/e)n)≤2nn_{M}^{*}(k,\Delta,\delta)\leq n\Rightarrow n_{P}^{*}(k,\Delta,\delta+(2/e)^{n})\leq 2n.

nB∗(k,Δ,δ)≤n⇒nH∗(k,Δ,δ+(e/4)n)≤2nn_{B}^{*}(k,\Delta,\delta)\leq n\Rightarrow n_{H}^{*}(k,\Delta,\delta+(e/4)^{n})\leq 2n;

nH∗(k,Δ,δ)≤n⇒nB∗(k,Δ,δ+(2/e)n)≤2nn_{H}^{*}(k,\Delta,\delta)\leq n\Rightarrow n_{B}^{*}(k,\Delta,\delta+(2/e)^{n})\leq 2n.

Appendix B Correlation decay between fingerprints

The correlation coefficient between Φ0\Phi_{0} and Φj\Phi_{j} follows immediately:

where λi=npi\lambda_{i}=np_{i}. Note that max⁡x>0e−xxjj!=e−jjjj!→0\max_{x>0}\frac{e^{-x}x^{j}}{j!}=\frac{e^{-j}j^{j}}{j!}\rightarrow 0 as j→∞j\to\infty. Therefore, for any x>0x>0,

where oj(1)o_{j}(1) is uniform as j→∞j\to\infty. Taking derivative, the function x↦e−2xxj1−e−xx\mapsto\frac{e^{-2x}x^{j}}{1-e^{-x}} on x>0x>0 is increasing if and only if x+ex(j−2x)−j>0x+e^{x}(j-2x)-j>0, and the maximum is attained at x=j/2+oj(1)x=j/2+o_{j}(1). Therefore, applying j!>(j/e)jj!>(j/e)^{j},

Appendix C Proof of auxiliary lemmas

Taking mm-th derivative of pmp_{m}, 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 γ\gamma is Euler’s constant, RR is the unique positive solution to h′(x)=0h^{\prime}(x)=0 with h(x)≜log⁡Γ(x+n+1)Γ(x+1)xmh(x)\triangleq\log\frac{\Gamma(x+n+1)}{\Gamma(x+1)x^{m}}, H=R2h′′(R)H=R^{2}h^{\prime\prime}(R), and all o(1)o(1) terms are uniform in mm. In the following we consider each range separately and prove the non-asymptotic approximation in (37).

Case I. For 1≤m≤log⁡n1\leq m\leq\sqrt{\log n}, Stirling’s approximation gives

Case III. For log⁡n≤m≤n−n1/3\sqrt{\log n}\leq m\leq n-n^{1/3}, note that h(x)=∑i=1nlog⁡(x+i)−mlog⁡xh(x)=\sum_{i=1}^{n}\log(x+i)-m\log x, and thus H=R2h′′(R)=m−∑i=1nR2(R+i)2≤mH=R^{2}h^{\prime\prime}(R)=m-\sum_{i=1}^{n}\frac{R^{2}}{(R+i)^{2}}\leq m. By [MW58, Lemma 4.1], H=ω(1)H=\omega(1) in this range. Hence,

where RR is the solution to x(1x+1+⋯+1x+n)=mx(\frac{1}{x+1}+\dots+\frac{1}{x+n})=m. Bounding the sum by integrals, we have

If log⁡n≤m≤ne\sqrt{\log n}\leq m\leq\frac{n}{e}, then R≍mlog⁡(n/m)R\asymp\frac{m}{\log(n/m)}, and hence

In view of (58), we have ∣s(n+1,m+1)∣=n!(Θ(R))m|s(n+1,m+1)|=\frac{n!}{(\Theta(R))^{m}}, which is exactly (37) when m≤n/em\leq n/e. If n/e≤m≤n−n1/3n/e\leq m\leq n-n^{1/3}, then R≍n2n−mR\asymp\frac{n^{2}}{n-m}, and

Combining (58) yields that ∣s(n+1,m+1)∣=n!(Θ(1n))m|s(n+1,m+1)|=n!(\Theta(\frac{1}{n}))^{m}, which coincides with (37) since n≍mn\asymp m 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.

References