Properly Learning Poisson Binomial Distributions in Almost Polynomial Time

Introduction

The Poisson binomial distribution (PBD) is the discrete probability distribution of a sum of mutually independent Bernoulli random variables. PBDs comprise one of the most fundamental nonparametric families of discrete distributions. They have been extensively studied in probability and statistics [Poi37, Che52, Hoe63, DP09b], and are ubiquitous in various applications (see, e.g., [CL97] and references therein). Recent years have witnessed a flurry of research activity on PBDs and generalizations from several perspectives of theoretical computer science, including learning [DDS12, DDO+13, DKS15b, DKT15, DKS15a], pseudorandomness and derandomization [GMRZ11, BDS12, De15, GKM15], property testing [AD15, CDGR15], and computational game theory [DP07, DP09a, DP14a, DP14b, GT14].

Despite their seeming simplicity, PBDs have surprisingly rich structure, and basic questions about them can be unexpectedly challenging to answer. We cannot do justice to the probability literature studying the following question: Under what conditions can we approximate PBDs by simpler distributions? See Section 1.2 of [DDS15] for a summary. In recent years, a number of works in theoretical computer science [DP07, DP09a, DDS12, DP14a, DKS15b] have studied, and essentially resolved, the following questions: Is there a small set of distributions that approximately cover the set of all PBDs? What is the number of samples required to learn an unknown PBD?

We study the following natural computational question: Given independent samples from an unknown PBD $\mathbf{P}$, can we efficiently find a hypothesis PBD $\mathbf{Q}$ that is close to $\mathbf{P}$, in total variation distance? That is, we are interested in properly learning PBDs, a problem that has resisted recent efforts [DDS12, DKS15b] at designing efficient algorithms. In this work, we propose a new approach to this problem that leads to a significantly faster algorithm than was previously known. At a high-level, we establish an interesting connection of this problem to algebraic geometry and polynomial optimization. By building on this connection, we provide a new structural characterization of the space of PBDs, on which our algorithm relies, that we believe is of independent interest. In the following, we motivate and describe our results in detail, and elaborate on our ideas and techniques.

Distribution Learning. We recall the standard definition of learning an unknown probability distribution from samples [KMR+94, DL01]: Given access to independent samples drawn from an unknown distribution $\mathbf{P}$ in a given family ${\cal C}$, and an error parameter $\epsilon>0$, a learning algorithm for ${\cal C}$ must output a hypothesis $\mathbf{H}$ such that, with probability at least $9/10$, the total variation distance between $\mathbf{H}$ and $\mathbf{P}$ is at most $\epsilon$. The performance of a learning algorithm is measured by its sample complexity (the number of samples drawn from $\mathbf{P}$) and its computational complexity.

In non-proper learning (density estimation), the goal is to output an approximation to the target distribution without any constraints on its representation. In proper learning, we require in addition that the hypothesis $\mathbf{H}$ is a member of the family ${\cal C}$. Note that these two notions of learning are essentially equivalent in terms of sample complexity (given any accurate hypothesis, we can do a brute-force search to find its closest distribution in ${\cal C}$), but not necessarily equivalent in terms of computational complexity. A typically more demanding notion of learning is that of parameter estimation. The goal here is to identify the parameters of the unknown model, e.g., the means of the individual Bernoulli components for the case of PBDs, up to a desired accuracy $\epsilon$.

Discussion. In many learning situations, it is desirable to compute a proper hypothesis, i.e., one that belongs to the underlying distribution family ${\cal C}$. A proper hypothesis is typically preferable due to its interpretability. In the context of distribution learning, a practitioner may not want to use a density estimate, unless it is proper. For example, one may want the estimate to have the properties of the underlying family, either because this reflects some physical understanding of the inference problem, or because one might only be using the density estimate as the first stage of a more involved procedure. While parameter estimation may arguably provide a more desirable guarantee than proper learning in some cases, its sample complexity is typically prohibitively large.

For the class of PBDs, we show (Proposition 14, Appendix A) that parameter estimation requires $2^{\Omega(1/\epsilon)}$ samples, for PBDs with $n=\Omega(1/\epsilon)$ Bernoulli components, where $\epsilon>0$ is the accuracy parameter. In contrast, the sample complexity of (non-)proper learning is known to be $\widetilde{O}(1/\epsilon^{2})$ [DDS12]. Hence, proper learning serves as an attractive middle ground between non-proper learning and parameter estimation. Ideally, one could obtain a proper learner for a given family whose running time matches that of the best non-proper algorithm.

This work is part of a broader agenda of systematically investigating the computational complexity of proper distribution learning. We believe that this is a fundamental goal that warrants study for its own sake. The complexity of proper learning has been extensively investigated in the supervised setting of PAC learning Boolean functions [KV94, Fel15], with several algorithmic and computational intractability results obtained in the past couple of decades. In sharp contrast, very little is known about the complexity of proper learning in the unsupervised setting of learning probability distributions.

2 Our Results and Comparison to Prior Work.

We are ready to formally describe the main contributions of this paper. As our main algorithmic result, we obtain a near-sample optimal and almost polynomial-time algorithm for properly learning PBDs:

In this work, we circumvent the $(1/\epsilon)^{\Omega(\log(1/\epsilon))}$ cover size lower bound by establishing a new structural characterization of the space of PBDs. Very roughly speaking, our structural result allows us to reduce the proper learning problem to the case that the underlying PBD has $O(\log(1/\epsilon))$ distinct parameters. Indeed, as a simple corollary of our main structural result (Theorem 4 in Section 2), we obtain the following:

We note that in subsequent work [DKS15a] the authors generalize the above theorem to Poisson multinomial distributions.

Remark. We remark that Theorem 2 is quantitatively tight, i.e., $O(\log(1/\epsilon))$ distinct parameters are in general necessary to $\epsilon$-approximate PBDs. This follows directly from the explicit cover lower bound construction of [DKS15b].

We view Theorem 2 as a natural structural result for PBDs. Alas, its statement does not quite suffice for our algorithmic application. While Theorem 2 guarantees that $O(\log(1/\epsilon))$ distinct parameters are enough to consider for an $\epsilon$-approximation, it gives no information on the multiplicities these parameters may have. In particular, the upper bound on the number of different combinations of multiplicities one can derive from it is $(1/\epsilon)^{O(\log(1/\epsilon))}$, which is not strong enough for our purposes. The following stronger structural result (see Theorem 4 and Lemma 5 for detailed statements) is critical for our improved proper algorithm:

Now suppose we would like to properly learn an unknown PBD with $O(\log(1/\epsilon))$ distinct parameters and known multiplicities for each parameter. Even for this very restricted subset of PBDs, the construction of [DKS15b] implies a cover lower bound of $(1/\epsilon)^{\Omega(\log(1/\epsilon))}$. To handle such PBDs, we combine ingredients from Fourier analysis and algebraic geometry with careful Taylor series approximations, to construct an appropriate system of low-degree polynomial inequalities whose solution approximately recovers the unknown distinct parameters.

In the following subsection, we provide a detailed intuitive explanation of our techniques.

3 Techniques.

The starting point of this work lies in the non-proper learning algorithm from our recent work [DKS15b]. Roughly speaking, our new proper algorithm can be viewed as a two-step process: We first compute an accurate non-proper hypothesis $\mathbf{H}$ using the algorithm in [DKS15b], and we then post-process $\mathbf{H}$ to find a PBD $\mathbf{Q}$ that is close to $\mathbf{H}$. We note that the non-proper hypothesis $\mathbf{H}$ output by [DKS15b] is represented succinctly via its Discrete Fourier Transform; this property is crucial for the computational complexity of our proper algorithm. (We note that the description of our proper algorithm and its analysis, presented in Section 3, are entirely self-contained. The above description is for the sake of the intuition.)

We now proceed to explain the connection in detail. The crucial fact, established in [DKS15b] for a more general setting, is that the Fourier transform of a PBD has small effective support (and in particular the effective support of the Fourier transform has size roughly inverse to the effective support of the PBD itself). Hence, in order to learn an unknown PBD $\mathbf{P}$, it suffices to find another PBD, $\mathbf{Q}$, with similar mean and standard deviation to $\mathbf{P}$, so that the Fourier transform of $\mathbf{Q}$ approximates the Fourier transform of $\mathbf{P}$ on this small region. (The non-proper algorithm of [DKS15b] for PBDs essentially outputs the empirical DFT of $\mathbf{P}$ over its effective support.)

Note that the Fourier transform of a PBD is the product of the Fourier transforms of its individual component variables. By Taylor expanding the logarithm of the Fourier transform, we can write the log Fourier transform of a PBD as a Taylor series whose coefficients are related to the moments of the parameters of $\mathbf{P}$ (see Equation (2)). We show that for our purposes it suffices to find a PBD $\mathbf{Q}$ so that the first $O(\log(1/\epsilon))$ moments of its parameters approximate the corresponding moments for $\mathbf{P}$. Unfortunately, we do not actually know the moments for $\mathbf{P}$, but since we can easily approximate the Fourier transform of $\mathbf{P}$ from samples, we can derive conditions that are sufficient for the moments of $\mathbf{Q}$ to satisfy. This step essentially gives us a system of polynomial inequalities in the moments of the parameters of $\mathbf{Q}$ that we need to satisfy.

Unfortunately, the above structural result is not strong enough, as in order to set up an appropriate system of polynomial inequalities for the parameters of $\mathbf{Q}$, we must first guess the multiplicities to which the distinct parameters appear. A simple counting argument shows that there are roughly $k^{\log(1/\epsilon)}$ ways to choose these multiplicities. To overcome this second obstacle, we need the following refinement of our structural result on distinct parameters: We divide the parameters of $\mathbf{Q}$ into categories based on how close they are to or $1$. We show that there is a tradeoff between the number of parameters in a given category and the number of distinct parameters in that category (see Theorem 4). With this more refined result in hand, we show that there are only $(1/\epsilon)^{O(\log\log(1/\epsilon))}$ many possible collections of multiplicities that need to be considered (see Lemma 5]).

Given this stronger structural characterization, our proper learning algorithm is fairly simple. We enumerate over the set of possible collections of multiplicities as described above. For each such collection, we set up a system of polynomial equations in the distinct parameters of $\mathbf{Q}$, so that solutions to the system will correspond to PBDs whose distinct parameters have the specified multiplicities which are also $\epsilon$-close to $\mathbf{P}$. For each system, we attempt to solve it using Renegar’s algorithm. Since there exists at least one PBD $\mathbf{Q}$ close to $\mathbf{P}$ with such a set of multiplicities, we are guaranteed to find a solution, which in turn must describe a PBD close to $\mathbf{P}$.

4 Related Work.

Distribution learning is a classical problem in statistics with a rich history and extensive literature (see e.g., [BBBB72, DG85, Sil86, Sco92, DL01]). During the past couple of decades, a body of work in theoretical computer science has been studying these questions from a computational complexity perspective; see e.g., [KMR+94, FM99, AK01, CGG02, VW02, FOS05, BS10, KMV10, MV10, DDS12, DDO+13, CDSS14a, CDSS14b, ADLS15].

We remark that the majority of the literature has focused either on non-proper learning (density estimation) or on parameter estimation. Regarding proper learning, a number of recent works in the statistics community have given proper learners for structured distribution families, by using a maximum likelihood approach. See e.g., [DR09, GW09, Wal09, DW13, CS13, KS14, BD14] for the case of continuous log-concave densities. Alas, the computational complexity of these approaches has not been analyzed. Two recent works [ADK15, CDGR15] yield computationally efficient proper learners for discrete log-concave distributions, by using an appropriate convex formulation. Proper learning has also been recently studied in the context of mixture models [FOS05, DK14, SOAJ14, LS15]. Here, the underlying optimization problems are non-convex, and efficient algorithms are known only when the number of mixture components is small.

5 Organization.

In Section 2, we prove our main structural result, and in Section 3, we describe our algorithm and prove its correctness. In Section 4, we conclude with some directions for future research.

Main Structural Result

In this section, we prove our main structural results thereby establishing Theorems 2 and 3. Our proofs rely on an analysis of the Fourier transform of PBDs combined with recent results from algebraic geometry on the solution structure of systems of symmetric polynomial equations. We show the following:

Theorem 4 implies that one needs to only consider $(1/\epsilon)^{O(\log\log(1/\epsilon))}$ different combinations of multiplicities:

For every $\mathbf{P}$ as in Theorem 4, there exists an explicit set $\mathcal{M}$ of multisets of triples $(m_{i},a_{i},b_{i})_{1\leq i\leq k}$ so that

For each element of $\mathcal{M}$ and each $i$, $[a_{i},b_{i}]$ is either one of the intervals $I_{i}$ or $J_{i}$ as in Theorem 4 or $$or$$.

For each element of $\mathcal{M}$, $k=O(\log(1/\epsilon))$.

This is proved in Appendix B.1 by a simple counting argument. We multiply the number of multiplicities for each interval, which is at most the maximum number of parameters to the power of the maximum number of distinct parameters in that interval, giving $(1/\epsilon)^{O(\log\log(1/\epsilon))}$ possibilities.

We now proceed to prove Theorem 4. We will require the following result from algebraic geometry:

As an immediate corollary, we obtain the following:

The basic idea of the proof will be to show that the Fourier transforms of $\mathbf{P}$ and $\mathbf{Q}$ are close to each other. In particular, we will need to make use of the following intermediate lemma:

The proof of this lemma, which is given in Appendix B.2, is similar to (part of) the correctness analysis of the non-proper learning algorithm in [DKS15b].

We proceed by means of Lemma 9. We need only show that for all $\xi$ with $|\xi|=O(\log(1/\epsilon))$ that $|\widehat{\mathbf{P}}(\xi)-\widehat{\mathbf{Q}}(\xi)|\ll\epsilon/\sqrt{\log(1/\epsilon)}.$ For this we note that

Taking a logarithm and Taylor expanding, we find that

A similar formula holds for $\log(\widehat{\mathbf{Q}}(\xi))$. Therefore, we have that

An application of Lemma 9 completes the proof. ∎

The basic idea of the proof is as follows. First, we will show that it is possible to modify $\mathbf{P}$ in order to satisfy (ii) without changing its mean, increasing its variance (or decreasing it by too much), or changing it substantially in total variation distance. Next, for each of the other intervals $I_{i}$ or $J_{i}$, we will show that it is possible to modify the parameters that $\mathbf{P}$ has in this interval to have the appropriate number of distinct parameters, without substantially changing the distribution in variation distance. Once this holds for each $i$, conditions (iii) and (iv) will follow automatically.

The expectation of $\mathbf{P}$ remains unchanged.

The total variation distance between the old and new distributions is $O(pp^{\prime})$, as is the change in variances between the distributions.

The variance of $\mathbf{P}$ is decreased.

Proper Learning Algorithm

In order to set up the necessary system of polynomial equations, we have the following theorem:

Consider another PBD with parameters $q_{i}$ of multiplicity $m_{i}$ contained in intervals $[a_{i},b_{i}]$ as described in Theorem 4. There exists an explicit system $\mathcal{P}$ of $O(\log(1/\epsilon))$ real polynomial inequalities each of degree $O(\log(1/\epsilon))$ in the $q_{i}$ so that:

Furthermore, such a system can be found with rational coefficients of encoding size $O(\log^{2}(1/\epsilon))$ bits.

Next, we need a low-degree polynomial to express the condition that Fourier coefficients of $\mathbf{Q}$ are approximately correct. To do this, we let $S$ denote the set of indices $i$ so that $[a_{i},b_{i}]\subset[0,1/2]$ and $T$ the set so that $[a_{i},b_{i}]\subset[1/2,1]$ and let $m=\sum_{i\in T}m_{i}$. We let

be an approximation to the logarithm of $\widehat{\mathbf{Q}}(\xi)$. We next define $\exp^{\prime}$ to be a Taylor approximation to the exponential function

We complete our system $\mathcal{P}$ with the final inequality:

In order for our analysis to work, we will need for $q_{\xi}$ to approximate $\widehat{\mathbf{Q}}(\xi)$. Thus, we make the following claim:

This is proved in Appendix C by showing that $g_{\xi}$ is close to a branch of the logarithm of $\widehat{\mathbf{Q}}(\xi)$ and that $|g_{\xi}+2\pi io_{\xi}|\leq O(\log(1/\epsilon))$, so $\exp^{\prime}$ is a good enough approximation to the exponential.

Hence, our system $\mathcal{P}$ is defined as follows:

$q_{i}$ for each distinct parameter $i$ of $\mathbf{Q}$.

Equations: Equations (4), (5), (6), (7), (8), and (9).

As we have defined it so far, the system $\mathcal{P}$ does not have rational coefficients. Equation (7) makes use of $e(\pm\xi/M)$ and $\pi$, as does Equation (8). To fix this issue, we note that if we approximate the appropriate powers of $(\pm 1\pm e(\pm\xi/M))$ and $q\pi i$ each to accuracy $(\epsilon/\sum_{i\in S}m_{i}))^{10}$, this produces an error of size at most $\epsilon^{4}$ in the value $g_{\xi}$, and therefore an error of size at most $\epsilon^{3}$ for $q_{\xi}$, and this leaves the above argument unchanged.

We then slightly modify Equation (8), replacing it by

Note that by our bound on $\sum_{i\in R}m_{i}$, this is of degree $O(\log(1/\epsilon))$.

We now need only prove the analogue of Claim 12 in order for the rest of our analysis to follow.

We prove this in Appendix C, by proving similar bounds to those needed for Claim 12. This completes the proof of our theorem in the second case. ∎

Our algorithm for properly learning PBDs is given in pseudocode below:

It is clear that the sample complexity of our algorithm is $O(\epsilon^{-2}\log^{2}(1/\epsilon))$. The runtime of the algorithm is dominated by Step 5. We note that by Lemma 5, $|\mathcal{M}|=(1/\epsilon)^{O(\log\log(1/\epsilon))}$. Furthermore, by Theorems 10 and 11, the runtime for solving the system $\mathcal{P}_{m}$ is $O(\log(1/\epsilon))^{O(\log(1/\epsilon))}=(1/\epsilon)^{O(\log\log(1/\epsilon))}$. Therefore, the total runtime is $(1/\epsilon)^{O(\log\log(1/\epsilon))}$.

Conclusions and Open Problems

A related open question concerns obtaining faster proper algorithms for learning more general families of discrete distributions that are amenable to similar techniques, e.g., sums of independent integer-valued random variables [DDO+13, DKS15b], and Poisson multinomial distributions [DKT15, DKS15a]. Here, we believe that progress is attainable via a generalization of our techniques.

The recently obtained cover size lower bound for PBDs [DKS15b] is a bottleneck for other non-convex optimization problems as well, e.g., the problem of computing approximate Nash equilibria in anonymous games [DP14b]. The fastest known algorithms for these problems proceed by enumerating over an $\epsilon$-cover. Can we obtain faster algorithms in such settings, by avoiding enumeration over a cover?

References

Appendix

Appendix A Sample Complexity Lower Bound for Parameter Estimation

Suppose that $n\geq 1/\epsilon$. Any learning algorithm that takes $N$ samples from an $n$-PBD and returns estimates of these parameters to additive error at most $\epsilon$ with probability at least $2/3$ must have $N\geq 2^{\Omega(1/\epsilon)}$.

We may assume that $n=\Theta(1/\epsilon)$ (as we could always make the remaining parameters all ) and demonstrate a pair of PBDs whose parameters differ by $\Omega(\epsilon)$, and yet have variation distance $2^{-\Omega(1/\epsilon)}$. Therefore, if such an algorithm is given one of these two PBDs, it will be unable to distinguish which one it is given, and therefore unable to learn the parameters to $\epsilon$ accuracy with at least $2^{\Omega(1/\epsilon)}$ samples.

In order to make this construction work, we take $\mathbf{P}$ to have parameters $p_{j}:=(1+\cos\left(\frac{2\pi j}{n}\right))/8$, and let $\mathbf{Q}$ have parameters $q_{j}:=(1+\cos\left(\frac{2\pi j+\pi}{n}\right))/8$. Suppose that $j=n/4+O(1)$. We claim that none of the $q_{i}$ are closer to $p_{j}$ that $\Omega(1/n)$. This is because for all $i$ we have that $\left(\frac{2\pi i+\pi}{n}\right)$ is at least $\Omega(1/n)$ from $\left(\frac{2\pi j}{n}\right)$ and $\left(\frac{2\pi(n-j)}{n}\right)$.

Appendix B Omitted Proofs from Section 2

For completeness, we restate the lemma below.

Lemma 5. For every $\mathbf{P}$ as in Theorem 4, there exists an explicit set $\mathcal{M}$ of multisets of triples $(m_{i},a_{i},b_{i})_{1\leq i\leq k}$ so that

For each element of $\mathcal{M}$ and each $i$, $[a_{i},b_{i}]$ is either one of the intervals $I_{i}$ or $J_{i}$ as in Theorem 4 or $$or$$.

For each element of $\mathcal{M}$, $k=O(\log(1/\epsilon))$.

For this choice of $\mathcal{M}$, (i) is automatically satisfied, and (iii) follows immediately from Theorem 4. To see (ii), we note that the total number of term in an element of $\mathcal{M}$ is at most

B.2 Proof of Lemma 9.

For completeness, we restate the lemma below.

The proof of this lemma is similar to the analysis of correctness of the non-proper learning algorithm in [DKS15b].

By Plancherel’s Theorem, the RHS above is

Appendix C Omitted Proofs from Section 3

In this section, we prove Claims 12 and 13 which we restate here.

Let $\mathbf{Q}^{\prime}$ be the PBD obtained from $\mathbf{Q}$ upon removing all parameters corresponding to elements of $R$. We note that

Therefore, it suffices to prove our claim when $R=\emptyset$.