An optimal $(\epsilon,\delta)$-approximation scheme for the mean of random variables with bounded relative variance

TL;DR

This paper introduces an efficient $(\epsilon,\delta)$-approximation scheme for estimating the mean of random variables with bounded relative variance, matching optimal sample complexity without requiring higher moment bounds.

Contribution

It presents a simple, easy-to-implement approximation method that achieves optimal sample complexity for mean estimation under bounded relative variance.

Findings

Achieves optimal sample complexity of approximately $2c^2\epsilon^{-2}\ln(1/\delta)$ samples.

Does not require bounds on third or fourth moments.

Provides a practical scheme for randomized approximation in complex counting problems.

Abstract

Randomized approximation algorithms for many #P-complete problems (such as the partition function of a Gibbs distribution, the volume of a convex body, the permanent of a $\{0,1\}$-matrix, and many others) reduce to creating random variables $X_1,X_2,\ldots$ with finite mean $\mu$ and standard deviation$\sigma$ such that $\mu$ is the solution for the problem input, and the relative standard deviation $|\sigma/\mu| \leq c$ for known $c$. Under these circumstances, it is known that the number of samples from the $\{X_i\}$ needed to form an $(\epsilon,\delta)$-approximation $\hat \mu$ that satisfies $\mathbb{P}(|\hat \mu - \mu| > \epsilon \mu) \leq \delta$ is at least $(2-o(1))\epsilon^{-2} c^2\ln(1/\delta)$. We present here an easy to implement $(\epsilon,\delta)$-approximation $\hat \mu$ that uses $(2+o(1))c^2\epsilon^{-2}\ln(1/\delta)$ samples. This achieves the same optimal running time as other estimators, but without the need for extra conditions such as bounds on third or fourth moments.

Figures

Equations

$$\left|\frac{\sigma }{\mu }\right|\le c.$$

$$\mathbb{P}(|\hat{\mu }-\mu |>ϵ\mu )\le \delta .$$

$$(S+o(1))c^2{ϵ}^{-2}\ln (1/\delta )$$

$$n=\lceil \frac{2c^2{ϵ}^{-2}\ln (4/\delta )}{1-ϵ}\rceil +\lceil 8{ϵ}^{-1}\left(1+c^2\right)\rceil \cdot \left[2\lceil \ln \left(\frac{7}{48\sqrt{\pi }\delta }\right)/\ln \left(\frac{16}{7}\right)\rceil +1\right].$$

$$n\ge 2{ϵ}^{-2}{c}^2\left[\ln \left(\frac{1}{\sqrt{2\pi }\delta }\right)-\ln \left(\frac{2\ln (1/[\sqrt{2\pi }\delta ])+1}{\sqrt{2\ln (1/[\sqrt{2\pi }\delta ])}}\right)\right].$$

$$\Psi (u)=\ln (1+u+u^2/2)\text{mathds}1(u\ge 0)-\ln (1-u+u^2/2)\text{mathds}1(u\le 0).$$

$$X_i={\hat{\mu }}_1+(X_i-{\hat{\mu }}_i)$$

$$W_i={\hat{\mu }}_1+{\alpha }^{-1}\Psi (\alpha \cdot (X_i-{\hat{\mu }}_i)),$$

$$W_i={\hat{\mu }}_1+{\alpha }^{-1}\Psi (\alpha \cdot (X_i-{\hat{\mu }}_1)).$$

$$\mathbb{P}(B\notin [a,b])\le 2\frac{4^r}{\sqrt{\pi r}}\cdot \frac{[(1-(b-a))(b-a){]}^{r+1}}{2(b-a)-1}$$

$${\int }_{x\in [0,a]\cup [b,1]}{f}_B(r)dr\le {\int }_{x\in [0,a]\cup [a,a+1-b]}{f}_B(r)={\int }_{x\in [0,1-(b-a)]}{f}_B(r)dr.$$

$$(\forall x\in [0,t])(x-x^2)\le t(1-t)+(x-t)(1-2t).$$

$${\int }_{x\in [0,t]}[t(1-t)+(x-t)(1-2t){]}^{r+1}dr\le \frac{[t(1-t){]}^r}{(r+1)(1-2t)}.$$

$$\mathbb{P}(\mathrm{med}(\{A_i\})\notin [\mu (1-ϵ),\mu (1+ϵ)])<\frac{({\nu }^2)(1-{\nu }^2)}{\sqrt{\pi r}(1-2{\nu }^2)}\exp (r\ln (4({\nu }^2)(1-{\nu }^2))).$$

$$\mathbb{P}(A_i\notin [\mu (1-ϵ),\mu (1+ϵ)])\le \frac{{\nu }^2{ϵ}^2{\mu }^2}{(ϵ\mu {)}^2}={\nu }^2.$$

$$\mathbb{P}(\mathrm{med}(\{U_i\})\in [x_1,x_2])=\mathbb{P}(\mathrm{med}(\{A_i\})\in [a,b]).$$

$$\mathbb{P}(\mathrm{med}(\{U_i\})\notin [x_1,x_2])\le \frac{({\nu }^2)(1-{\nu }^2)[4({\nu }^2)(1-{\nu }^2){]}^r}{\sqrt{\pi r}(1-2{\nu }^2)}$$

$$S=(X_1+\cdots +X_k)/k.$$

$$\mathbb{P}(|\mathrm{med}(\{S_i\})-\mu |>ϵ\mu )\le \delta .$$

$${\Psi }_U(x)=\ln (1+x+x^2/2),{\Psi }_L(x)=-\ln (1-x+x^2/2).$$

$${\Psi }_L(x)\le \Psi (x)\le {\Psi }_U(x).$$

$${\Psi }_L(x)\le \Psi (x)$$

$$\iff 1\le 1+x^4/4,$$

$$W_{L,i}$$

$$W_{U,i}$$

$$\mathbb{P}(\bar{W}>\mu (1+ϵ))\le \exp [-(n{ϵ}^2/(2c^2))\cdot (1-{\xi }^2(1+1/c^2)))].$$

$$\mathbb{P}({\bar{W}}_U>\mu (1+ϵ)$$

$$=\mathbb{P}(\exp (\alpha (W_{U,1}+\cdots +W_{U,n}))>\exp (\alpha n\mu (1+ϵ)))$$

$$\le \mathbb{E}[\exp (\alpha (W_{U,1}+\cdots +W_{U,n}))]/\exp (\alpha n\mu (1+ϵ))$$

$$={\left[\frac{\mathbb{E}[\exp (\alpha W_{U,1})]}{\exp (\alpha \mu (1+ϵ)}\right]}^n$$

$$\exp (\alpha W_{U,1})$$

$$=\exp (\alpha {\hat{\mu }}_1)[1+\alpha (X-{\hat{\mu }}_1)+({\alpha }^2/2)(X-{\hat{\mu }}_1{)}^2]$$

$$=\exp (\alpha {\hat{\mu }}_1)[1+\alpha (X-\mu +\gamma )+({\alpha }^2/2)(X-\mu +\gamma {)}^2]$$

$$=\exp (\alpha {\hat{\mu }}_1)[1+\alpha (X-\mu )+\alpha \gamma +({\alpha }^2/2)((X-\mu {)}^2+2(X-\mu )\gamma +{\gamma }^2)]$$

Full text

1 Introduction

Suppose $X_{1},X_{2},\ldots\sim X$ are iid with mean $\mathbb{E}[X]=\mu$ and variance $\mathbb{V}(X)=\sigma^{2}$. The relative standard deviation is $\sigma/|\mu|$ and the relative variance is $\sigma^{2}/\mu^{2}$. Say the relative standard deviation is bounded by $c$ if

[TABLE]

Suppose $\mu$ and $\sigma$ are unknown, but $c$ is known. Then the goal is to use as few $X_{i}$ as possible to find an estimate $\hat{\mu}$ for $\mu$ that is an $(\epsilon,\delta)$-randomized approximation, that is

[TABLE]

Suppose that an $(\epsilon,\delta)$-randomized approximation requires

[TABLE]

samples for a constant $S$ and where the little-o notation refers to $\epsilon\rightarrow 0$. Then call $S$ the scale factor of the algorithm.

In this work we present a simple algorithm that is both easy to implement and which achieves the optimal scale factor $S=2$ without any additional assumptions about the random variables such as bounded higher moments.

This basic problem arises often in randomized algorithms. For instance, problems for approximating the partition function of the Ising model [14], the permanent of a $\{0,1\}$-matrix [17], the volume of a convex body [7, 18], the number of solutions to a DNF logical expression [16], the number of linear extensions of a poset [20, 11], and many more all have this problem as a subproblem. Any improvement in the ability to deal with this basic problem directly translates into better approximation algorithms for all of these problems.

This problem has a long history, stretching back to Nemirovsky and Yudin [23] who used the median-of-means estimator in the context of stochastic optimization. Jerrum, Valiant, and Vazirani [16] developed a similar estimator for the purposes of creating randomized approximation schemes for #P complete problems. By 1999 [1], this method was in wide use for online algorithms. Hsu and Sabato [9, 10] analyzed the basic median-of-means estimator and proved that it had a scale factor of 121.5 for small enough $\epsilon$.

Catoni [4] greatly advanced the area by presenting an approximation that used an $M$-estimator. This was not an $(\epsilon,\delta)$-randomized approximation algorithm, rather it gave a confidence bound based on the samples and specific values of the parameters used in the estimate. While it could bound the confidence interval for the estimate based on the parameters and the unknown $\mu$ and $\sigma^{2}$ for $X_{i}$, there seems to be no way of setting the parameters ahead of time for given $\epsilon$ and $\delta$ without some additional information. Certainly no such method was given in [4].

Devroye et. al. [6] showed that if the kurtosis of the random variables is bounded above, then the optimal scale factor $S=2$ could be attained with a simpler estimator. Unfortunately, in order to run their algorithm, the user needed this upper bound on the kurtosis. Bounding the kurtosis can be much more challenging mathematically than bounding the variance. Minsker and Strawn [21] returned to the original median-of-means estimator. When the random variables have bounded third moment, the Berry-Esseen Theorem can be used to show quick convergence to normality, and they showed that this gave the simple median-of-means algorithm a scale factor of 4.5. As with the Devroye et. al. method, this requires that the bound on the third moment be given explicitly before the algorithm can be used.

The approach here takes the Catoni $M$-estimator in a new direction. There is no unique approach to getting the extra information required to turn the Catoni $M$-estimator into an $(\epsilon,\delta)$-approximation. One approach is to use a two-step process that works as follows. Before running the $M$-estimator, first generate a weaker estimate $\hat{\mu}_{1}$ that is an $(\sqrt{\epsilon},\delta/2)$ randomized approximation to $\mu$. Then use this estimate to set the parameters of the Catoni $M$-estimator to give an output that is provably an $(\epsilon,\delta)$-approximation.

While this two-step process works, it (like all $M$-estimators) requires finding the root of a nonlinear equation. Analysis of the number of steps needed to get an close approximation to the root was not done in [4], and would need to be accomplished before the running time of the method is known.

Previous algorithms either had too large a scale factor, required rootfinding, or required knowledge of higher moments. The new method presented here solves all these difficulties.

• •

It achieves the optimal scale factor $S=2$.

• •

No rootfinding step is required. Instead, first a function is randomly chosen by some initial samples, and then the final estimator is a sample average this random function applied to new data.

• •

No bound on higher moments is necessary. In fact, even if the second moment is the highest moment that exists for the random variables, the new method is still an $(\epsilon,\delta)$-randomized approximation.

Our main result concerning this new method is as follows.

Theorem 1**.**

Let $X_{1},X_{2},\ldots\sim X$ with $\mathbb{E}[X]=\mu$ and $\mathbb{V}[X]=\sigma^{2}$ satisfying $\sigma^{2}/\mu^{2}\leq c^{2}.$ For $\epsilon<1$, there exists an $(\epsilon,\delta)$-randomized approximation algorithm that uses $n$ samples where

[TABLE]

This constant in the leading order term is the best possible.

Theorem 2**.**

Given $\epsilon$ and $\delta$ positive, let $\hat{\mu}:^{n}\rightarrow$ be an $(\epsilon,\delta)$-randomized algorithm for all distributions $X$ with $\mathbb{E}[X]=\mu$ and $\mathbb{V}(X)=\sigma^{2}$ satisfying $\sigma^{2}/\mu^{2}=c^{2}$. Then

[TABLE]

The remainder of the paper is organized as follows. The next section describes the two-step algorithm and proves correctness for each step. Section 3 then shows the lower bound on the number of samples needed, and Section 4 considers several of the applications mentioned in the introduction in more detail.

2 The Algorithm

Define the $\Psi$ function as follows.

[TABLE]

This function was used in [4] as part of the $M$-estimator.

For $u\in[-1,1]$, the value of $\Psi(u)$ is approximately $u$ (see Figure 1.) For $u$ greater than 1 in magnitude, the value of $\Psi(u)/u$ becomes close to 0. For a constant $\alpha>0$, $\alpha^{-1}\Psi(\alpha u)$ is a scaled version of $\Psi$ that is close to $u$ for $u\in[-1/\alpha,1/\alpha]$.

Suppose that $\hat{\mu}_{1}$ is an initial estimate for $\mu$. Then

[TABLE]

has mean $\mu$ but is also susceptible to outliers in the $X_{i}$ distribution. By replacing this with

[TABLE]

the value of $W_{i}$ will be close to $X_{i}$ when $|X_{i}-\hat{\mu}_{i}|\leq\alpha^{-1}$, but always has a light-tailed distribution because of the logarithm function.

The algorithm proceeds as follows. The first step uses a median-of-means approach to find $\hat{\mu}_{1}$ that is a $(\sqrt{\epsilon},\delta/2)$ approximation of $\mu$. Given that the first step did not fail, the next step then uses the sample average of the $W_{i}$ variables to create the final estimate $\hat{\mu}$ that is an $(\epsilon,\delta/2)$ approximation. The chance that either step fails is at most $\delta$.

The first step is to construct a median-of-means estimator for $\mu$ [16]. Let $\epsilon_{1}=\sqrt{\epsilon(c^{2}/(1+c^{2}))},$ $k=\lceil 8c^{2}\epsilon_{1}^{-2}\rceil,$ and $m=2\lceil\ln(7/[48\sqrt{\pi}\delta)/\ln(16/7)\rceil+1$. Let $S$ have the distribution of the sample average of $k$ independent draws from $X$. Draw $S_{1},\ldots,S_{m}\sim S$ independently, and let $\hat{\mu}_{1}=\text{median}(\{S_{i}\})/(1-\epsilon_{1}^{2})$. 2. 2.

Let $n=\lceil 2c^{2}\epsilon^{-2}\ln(4/\delta)/(1-\epsilon)\rceil$, and $\alpha=\epsilon/[c^{2}\hat{\mu}_{1}]$. Draw $X_{1},\ldots,X_{n}$ independently. For all $i$, let

[TABLE]

Set $\hat{\mu}=(W_{1}+\cdots+W_{n})/n$.

2.1 The first step of the algorithm

The first step of the algorithm is the powering method of Jerrum, Valiant, and Vazirani [16] applied to the sample averages. This technique was also used in [1], and later referred to as the median-of-means method [9]. These authors did not attempt to optimize the constants in their arguments, and so we repeat the proof here so we can see exactly how the choice of constant enters into the failure bound.

Suppose that we have random variables whose relative standard deviation is at most $\nu\epsilon$. What is the chance that the median of $m=2r+1$ draws from the random variable falls into $[\mu(1-\epsilon),\mu(1+\epsilon)]$?

To answer this question, first consider the probability that a beta distributed random variable with both parameters equal to an integer $r$ falls into a subinterval of $[0,1]$.

Lemma 1**.**

Let $B\sim\textsf{Beta}(r+1,r+1)$ denote a random variable with density $f_{B}(x)=[(2r+1)!/(r!r!)]x^{r}(1-x)^{r}\mathds{1}(x\in[0,1])$. For any $0\leq a\leq 1/2\leq b\leq 1$ with $1-(b-a)\leq 1/2$,

[TABLE]

Proof.

Density $f_{B}$ is symmetric about its unique local maximum at $1/2$, so

[TABLE]

Note $x^{r}(1-x)^{r}=(x-x^{2})^{r}$. Let $t=1-(b-a).$ Then $t\leq 1/2$, $[x-x^{2}]^{\prime}=1-2x>0$ and $[x-x^{2}]^{\prime\prime}=-2x\leq 0$, so the function lies below its tangent line at $t$. That is,

[TABLE]

Then

[TABLE]

Using Stirling’s formula to give $(2r+1)(2r)!/(r!r!)\leq(2r+1)4^{r}/\sqrt{\pi r}$ completes the proof. ∎

Lemma 2**.**

Let $A_{1},\ldots,A_{2r+1}$ be iid with mean $\mu$ and variance at most $\nu^{2}\epsilon^{2}\mu^{2}$ where $\nu^{2}\leq 1/2$. Then

[TABLE]

Proof.

By Chebyshev’s inequality,

[TABLE]

Let $x_{1}=\mathbb{P}(A_{i}<a)$ and $x_{2}=\mathbb{P}(A_{i}\leq b)$. Construct a uniform random variable over $[0,1]$ as follows. If $A_{i}<a$, then let $U_{i}\sim\textsf{Unif}([0,x_{1}))$. If $A_{i}\in[a,b]$, let $U_{i}\sim\textsf{Unif}([x_{1},x_{2}])$. Finally, if $A_{i}>b$, then let $U_{i}\sim\textsf{Unif}((x_{2},1])$. Note that

[TABLE]

The median of $2r+1$ iid uniform $[0,1]$ random variables is well known to have a beta distribution: $\operatorname{med}(\{U_{i}\})\sim\textsf{Beta}(r+1,r+1)$. So the previous lemma can be used to state

[TABLE]

∎

Hence the failure probability is going down exponentially at rate $\ln(4(\nu^{2})(1-\nu^{2})).$ Now for an integer $k$, consider $X_{1},\ldots,X_{k}\sim X$ iid, and

[TABLE]

Then $\mathbb{E}[S]=\mathbb{E}[X]=\mu$ and $\mathbb{V}[S]=k\mathbb{V}[X]/k^{2}=\sigma^{2}/k.$

In particular, for $\sigma^{2}/\mu^{2}\leq c^{2}$, and $k=\lceil\epsilon^{-2}\nu^{-2}c^{2}\rceil$, $\mathbb{V}[S]\leq(\nu\epsilon\mu)^{2}.$ To take the median of $2r+1$ draws of the sample average of $k$ draws from the $\{X_{i}\}$ takes $\Theta(kr)=\Theta(-1/(\nu^{2}\ln(4\nu^{2}(1-\nu^{2}))))$ samples.

Lemma 3** (Median-of-means).**

Suppose $X_{1},X_{2},\ldots$ are as in (1). For $k=\lceil 8c^{2}\epsilon^{-2}\rceil$ let $S$ be distributed as $(X_{1}+\cdots+X_{k})/k$. Let $m=2\lceil\ln(7/[48\sqrt{\pi}\delta])/\ln(16/7)\rceil+1.$ Let $S_{1},\ldots,S_{m}\sim S$ then

[TABLE]

Proof.

Just set $\nu^{2}=1/8$ in the previous lemma to get a rough minimum for the bound. ∎

Now suppose that instead of bounding $\hat{\mu}_{1}/\mu-1$, we wish to bound $\xi=\mu/\hat{\mu}_{1}-1$. The next lemma shows how to build a biased estimate where $|\xi|\leq\epsilon$ from an estimate with relative error at most $\epsilon$.

Lemma 4**.**

Suppose $\hat{\mu}_{1}$ is an estimate for $\mu$ with $|\hat{\mu}_{1}/\mu-1|\leq\epsilon$. Then $\hat{\mu}_{2}=\hat{\mu}_{1}/(1-\epsilon^{2})$ has $|\mu/\hat{\mu}_{2}-1|\leq\epsilon$.

Proof.

The proof follows from simplifying the appropriate inequalities. ∎

This is why in the first step of the algorithm, the estimate found from median-of-means is divided by $1-\epsilon_{1}^{2}$ before moving to the next step.

2.2 The second step of the algorithm

To analyze this step, it helps to have two new functions that upper and lower bound $\Psi$.

[TABLE]

Lemma 5**.**

For all $x\in$,

[TABLE]

Proof.

First consider $\Psi_{L}(x)\leq\Psi(x)$. These are equal when $x\leq 0$, so suppose $x\geq 0$. Exponentiating gives

[TABLE]

therefore the inequality holds. The other inequality is shown similarly. ∎

Now set

[TABLE]

By the previous lemma, $W_{L,i}\leq W_{i}\leq W_{U,i}$ for all $i$.

Lemma 6**.**

Denote $\bar{W}_{U}=(W_{U,1}+\cdots+W_{U,n})/n$. Then

[TABLE]

Proof.

Take a Chernoff bound [5] style approach. Since $\alpha>0$ and $\exp$ is a strictly increasing function,

[TABLE]

First consider the expression inside the mean in the numerator. Setting $\gamma=\mu-\hat{\mu}_{1}$ gives

[TABLE]

Since $\mathbb{E}[X-\mu]=0$ and $\mathbb{E}[(X-\mu)^{2}]=\sigma^{2}$, we have

[TABLE]

Note

[TABLE]

Next use $1+x\leq\exp(x)$ to state

[TABLE]

Since $\alpha=\epsilon/[c^{2}\hat{\mu}_{1}^{2}]$, $\gamma=\mu-\hat{\mu}_{1}$ and $\xi=\mu/\hat{\mu}_{1}-1$,

[TABLE]

Similarly, using $\sigma^{2}/\mu^{2}\leq c^{2}$,

[TABLE]

Note $2(1+\xi)-(1+\xi)^{2}=1-\xi^{2}$.

Putting this together with the $\alpha^{2}\gamma^{2}/2$ term gives

[TABLE]

∎

Note that at the end of Step 1 of the algorithm, $\xi^{2}\leq\epsilon(c^{2}/(1+c^{2}))$, which means $\xi^{2}(1+1/c^{2})\leq\epsilon$, and $\mathbb{P}(\bar{W}_{U}>\mu(1+\epsilon))\leq\exp(-n(1-\epsilon)\epsilon^{2}/[2c^{2}]).$

Lemma 7**.**

Denote $\bar{W}_{L}=(W_{L,1}+\cdots+W_{L,n})/n$. Then

[TABLE]

Proof.

The proof is similar to the previous lemma: first multiply by $-\alpha$ and exponentiate to get

[TABLE]

and the rest of the proof is the same as the previous lemma. ∎

Putting these results together gives the following.

Lemma 8**.**

For $n\geq 2c^{2}\epsilon^{-2}\ln(2/\delta)(1-\epsilon)^{-1}$ and $|\hat{\mu}_{1}-\mu|\leq\sqrt{\epsilon c^{2}/(1+c^{2})}\hat{\mu}_{1}$,

[TABLE]

Proof.

Apply the previous lemma using $\xi^{2}\leq\epsilon c^{2}/(1+c^{2}).$ ∎

Theorem 1 immediately follows.

3 Lower bound on the number of samples

Begin with a rephrasing of Proposition 6.1 from [4].

Lemma 9**.**

Let $\hat{\mu}:^{n}\rightarrow$ be any estimator of the mean of $n$ iid random variables. Let $Y_{1},\ldots,Y_{n}\sim\textsf{N}(\mu,\sigma^{2})$ and $\bar{Y}=(Y_{1}+\cdots+Y_{n})/n$. Then it holds that either $\mathbb{P}(\hat{\mu}\geq\mu(1+\epsilon))\geq\mathbb{P}(\bar{Y}\geq\mu(1+\epsilon))$ or $\mathbb{P}(\hat{\mu}\leq\mu(1-\epsilon))\geq\mathbb{P}(\bar{Y}\leq\mu(1-\epsilon))$ for $\bar{Y}=(Y_{1}+\cdots Y_{n})/n$.

In other words, for any estimator of the mean for normal random variables, there is either a higher chance that the estimate is in the upper tail than for the sample average, or there is a higher chance that the estimate falls in the lower tail than the sample average does. Note that for $Y_{i}\sim\textsf{N}(\mu,c^{2}\mu^{2})$, then $\bar{Y}\sim\textsf{N}(\mu,c^{2}\mu^{2}/n)$. Let $Z\sim\textsf{N}(0,1)$. From the scaling properties of normal random variables,

[TABLE]

Since $\mathbb{P}(Z\leq-a)=\mathbb{P}(Z\geq a)$ for all $a$, we need only bound one tail of the normal.

Lemma 10**.**

Let $Z\sim\textsf{N}(0,1)$ and $a_{\delta}$ satisfy $\mathbb{P}(Z\geq a_{\delta})=\delta$ where $\delta\leq 1/\sqrt{2\pi}$. Then

[TABLE]

Proof.

Gordon [8] showed that for $a\geq 0$,

[TABLE]

Without the $a/(a^{2}+1)$ factor, the right hand side equals $\delta/2$ when $a_{1}=\sqrt{2\ln(1/(\sqrt{2\pi}\delta))}$. Since $a/(a^{2}+1)\leq 1$ we have $a_{\delta}\leq a_{1}$. Also, $a/(a^{2}+1)$ is a decreasing function, so

[TABLE]

Solving gives

[TABLE]

as desired. ∎

Putting $a_{\delta}=\epsilon\sqrt{n}/c$ then gives Theorem 2.

4 Applications

Jerrum, Valiant, and Vazirani [16] showed that for a large class of self-reducible problems, the ability to sample from a density in polynomial time leads to an $(\epsilon,\delta)$-randomized approximation scheme for the normalizing constant of the unnormalized density. Since finding that normalizing constant is often a #P-complete problem, this has been used in many settings. Each of these leads to a problem such as that considered here where a random variable has mean $\mu$ equal to the target with bounded relative standard deviation. This method was expanded to more examples later by Jerrum and Sinclair [15].

The idea is as follows. Suppose that the goal is to find $\#A_{0}$ which is the size of a set (either number of elements for a finite set or the Lebesgue measure for $A_{0}\subset^{n}$.) Suppose that we can find a sequence of decreasing sets $A_{0}\supseteq A_{1}\supseteq A_{2}\supseteq\cdots\supseteq A_{k}$ where $\#A_{k}$ is known. If each of the sets $A_{i}$ represents an instance of the original problem (perhaps with a different input), the problem is self-reducible. If there is an efficient method for generating samples uniformly from the $A_{i}$, then for each $i\in\{0,1\ldots,k-1\}$, let $X_{i,1},\ldots,X_{i,m}\sim\textsf{Unif}(A_{i})$, and let $R_{i}=m^{-1}\sum_{j}\mathds{1}(X_{i,j}\in A_{i+1})$ be the percentage of values that fall into $A_{i+1}$. Then

[TABLE]

so let $\hat{r}=R_{0}\cdots R_{k-1}$ be the unbiased product estimator for $\#A_{k}/\#A_{0}$.

Then

[TABLE]

Let $r_{i}=\#A_{i}/\#A_{i+1}$. Then $X_{i,1}$ has a Bernoulli distribution with mean $r_{i}$ and variance $r_{i}(1-r_{i})$. As the sample average of $m$ iid draws from $X_{i,1}$, $R_{i}$ has mean $r_{i}$ and variance $r_{i}(1-r_{i})/m$. Then

[TABLE]

so if $r_{i}\geq 1/M$ for all $i$, using $1+x\leq\exp(x)$ gives

[TABLE]

There are $k$ different $R_{i}$ each requiring $m$ samples, therefore $km$ are needed to generate one value of $\hat{r}$. From the above the variance is $\exp(k(M-1)/m)-1\approx k(M-1)/m$ for large $m$. Hence for large $m$ (such as $k(M-1)\epsilon^{-2}$) using the algorithm presented here has the total number of samples needed for an $(\epsilon,\delta)$-approximation is (to leading order) $2k(M-1)\epsilon^{-2}\ln(4/\delta)$, with the 2 being the optimal value of the constant.

4.1 Linear extensions of a poset

For a direct application of this process, consider the problem of counting the number of linear extensions of a partially ordered set (poset). A poset on $n$ objects $\{1,\ldots,n\}$ is an ordering $\preceq$ with three properties. Let $i,j,k\in\{1,\ldots,n\}$. First, $i\preceq i$. Second, if $i\preceq j$ and $j\preceq i$, then $i=j$. Third, if $i\preceq j$ and $j\preceq k$, then $i\preceq k$. A linear extension of the poset is a permutation $\tau$ such that $\tau(i)\preceq\tau(j)\Rightarrow i\leq j$.

Brightwell and Winkler [2] showed that counting the number of linear extensions of an arbitrary poset is a #P-complete problem. Finding the number of linear extensions has applications in nonparametric statistics [22].

A sequence of results [19, 20, 3, 11] culminated in an $O(n^{3}\ln(n))$ method for generating samples uniformly from the set of linear extensions. To convert this method into a method for approximately counting the number of linear extensions, use self-reducibility.

Let $n_{\ell}$ be any element of $\{1,\ldots,n\}$ which is not preceded by another element in the set. Then an easy Markov chain argument gives that the probability that a uniformly chosen linear extension has $\tau(n_{\ell})=n$ is at least $1/n$. Fixing $\tau(n_{\ell})=n$ in the permutation leaves a linear extension problem of size $n-1$. So the methods of this section can be applied with $k=n$ and $M=n$. Hence (to first order) $2n^{2}\epsilon^{-2}\ln(4/\delta)$ samples are needed to give an $(\epsilon,\delta)$-approximation to the number of linear extensions.

4.2 Permanent of a $\{0,1\}$-matrix

Let ${\cal S}_{n}$ be the set of permutations on $\{1,\ldots,n\}$. Then the permanent of a matrix $A$ with entries $a_{ij}$ is

[TABLE]

Calculating the permanent exactly was shown by Valiant [25] to be a #P-complete problem.

Note that if $a_{ij}\in\{0,1\}$, then the only permutations $\tau$ that contribute to the sum have $a_{i,\tau(i)}=1$ for all $i$. So the permanent is the normalizing constant of the distribution over ${\cal S}_{n}$ with unnormalized density $f(\tau)=\prod_{i=1}^{n}a_{i,\tau(i)}.$

Jerrum, Sinclair, and Vigoda [17] developed a polynomial time algorithm for approximately sampling from the density $f(\tau)$. As with the previous problem of linear extensions, for such a problem on permutations there exists a value $i$ such that $\mathbb{P}(\tau(n)=i)\geq 1/n$. This can then be used with the basic self-reducibility process to get an $(\epsilon,\delta)$-approximation for the permanent. Without going into details (as the method of [17] for approximation was more complex than the basic approach) the result is the same as for linear extensions: use of the methods of this paper immediately reduces the constant in the leading term down to the optimal value.

4.3 Gibbs distributions

These distributions arise in statistical physics and other applications.

Definition 1**.**

$\{\pi_{\beta}\}_{\beta\in}$ is a Gibbs distribution with parameter $\beta$ over finite state space $\Omega$ if there exists a Hamiltonian function $H(x):\Omega\rightarrow$ such that for $X\sim\pi_{\beta}$,

[TABLE]

where $Z(\beta)=\sum_{x\in\Omega}\exp(-\beta H(x))$ is called the partition function of the distribution.

A famous example of a Gibbs distribution is the Ising model [13], where the state space consists of labellings of the nodes of a graph $G=(V,E)$ by either 0 or 1, and $H(x)=\sum_{\{v,w\}\in E}-(x(v)-x(w))^{2}.$ In [14] finding the partition function of the ferromagnetic Ising model (where $\beta>0$) was shown to be a #P-complete problem for general graphs, but that same work showed how to generate (approximately) samples from the distribution in time polynomial in the size of the graph.

Typically it is easy to find $Z(0)$ for these problems. For the Ising model, $Z(0)=2^{\#V}$. In [24] it was shown how to build an estimate for $Z_{\beta}/Z(0)$ using samples from $\pi$ where the ratio $\sigma^{2}/\mu^{2}$ was bounded. In [12], it was shown how to build two random variables $W$ and $V$ such that $\mathbb{E}[W]/\mathbb{E}[V]=Z(\beta)/Z(0)$ and each had relative variance bounded above by $2e$.

Let $\epsilon^{\prime}=[-1+\sqrt{1+\epsilon^{2}}]/\epsilon\leq\epsilon/2-\epsilon^{3}(1.5-\sqrt{2})$ for $\epsilon\in[0,1]$. If

[TABLE]

then

[TABLE]

Then it is straightforward to show that $\hat{\mu}=[\hat{\mu}_{W}/\hat{\mu}_{V}]\sqrt{1+\epsilon^{2}}$ satisfies

[TABLE]

thereby giving an $(\epsilon,\delta)$-approximation that (to leading order) requires $2(4\epsilon^{-2})(2e)^{2}\ln(4/\delta)$ samples to estimate the partition function value.