On the Structure, Covering, and Learning of Poisson Multinomial Distributions

Introduction

In this paper, we advance our understanding of the structure and learnability of this fundamental family of distributions by studying the following questions:

Can we approximate PMDs via simpler distributions such as multi-dimensional Gaussians or Poissons? Do they always “behave as” discretized multi-dimensional Gaussians or Poissons? If not, what is the range of possible “behaviors” that PMDs may exhibit?

Given $n$, $k$ and $\varepsilon$, is there a small set of distributions that $\varepsilon$-cover, in total variation distance, the set of all $(n,k)$-PMDs? And, how does the size of the cover scale with $n$, $k$ and $\varepsilon$?

How many samples from a $(n,k)$-PMD do we need to learn its density to within $\varepsilon$ in total variation distance? What is the dependence of the learning complexity on the size $n^{O(k)}$ of their support?

In summary, known bounds show that a $(n,k)$-PMD can be approximated by simpler, $\text{poly}(k)$-parameter, distributions, but the quality of their approximation depends on the first few moments of the PMD or its summands. Our goal instead is to provide universal approximation theorems showing how to approximate a given $(n,k)$-PMD by simpler distributions for any desired approximation $\varepsilon$ and without assumptions about the moments of the PMD or its summands. Our main structural theorem is the following.

By introducing the independent $({\rm poly}(k/\varepsilon),k)$-PMD, our structural result side-steps the degradation of the CLT bound of [VV11] with $\log n$ and the smallest eigenvalue of the PMD’s covariance matrix, correcting it to any desired approximation $\varepsilon$. Interestingly, there may be directions where the variance of the discretized Gaussian used in our result may be arbitrarily far from that of the approximated PMD. The sparse PMD added to the Gaussian serves to correct the variance in those directions, but does so in a correlated manner across several directions. Moreover, while [VV11] discretize their approximating multidimensional Gaussian to the closest lattice point, our discretization is more faithful to the structure of its covariance matrix; see Definition 6. We provide more intuition about our structural result in Section 1.1, where we also outline its proof. A more detailed proof of Theorem 1 appears in Section 3 and a more detailed statement is given as Theorem 5.

Covers for PMDs

Building covers for $(n,k)$-PMDs was pursued in [DP08, DP15] as a means to develop approximation algorithms for Nash equilibria in anonymous games. These are games where $n$ players share the same action set, say $\{1,\ldots,k\}$, and each player’s utility depends on their own choice of action as well as the distribution of how many of the other players choose each of the available actions, but players’ utility functions may otherwise be different. It was shown that proper $\varepsilon$-covers, in total variation distance, of $(n,k)$-PMDsAn $\varepsilon$-cover ${\cal F}_{\varepsilon}$ of a set of distributions ${\cal F}$ is called proper iff ${\cal F}_{\varepsilon}\subseteq{\cal F}$. imply approximation algorithms for Nash equilibria in these games, whose complexity scales with the size of the cover. Intuitively, this is because switching from a mixed Nash equilibrium to a mixed strategy profile with the same distribution of how many players choose each action does not affect players’ payoffs by more than $\varepsilon$.

The covers for $(n,k)$-PMDs obtained in the anonymous games papers cited above have size:

Such covers are of theoretical interest, their interesting feature being that the size is polynomial in $n$. Indeed, the standard discretization of the parameters of a PMD’s constituent vectors results in covers of size exponential in $n$, so a more delicate “global” discretization is needed to obtain covers whose size is polynomial in $n$.

Besides providing an asymptotically smaller search space for Nash equilibria in anonymous games, or any other optimization problem over PMDs, the polynomial rather than exponential dependence of the cover size on $n$ has direct consequences to the learnability of these distributions; see Theorem 7 (from [DK14]) and [AJOS14] for a similar result, which improve a long line of similar results in the probability literature [DL01]. In particular, a cover of polynomial size implies directly that these distributions can be learned from a number of samples logarithmic in $n$, despite their support being polynomial in $n$. Motivated by such applications of covers to algorithms and learning we use our structural result to obtain an improved cover theorem.

We make a few remarks about our cover. First, the cover is non-proper, containing distributions that are of the form specified in Theorem 1, i.e. are convolutions of a discretized Gaussian and a PMD. Moreover, it is straightforward to see that any cover has size at least $n^{\Omega(k)}$ and at least $(1/\varepsilon)^{\Omega(k)}$. For the first lower bound, count the number of $(n,k)$-PMDs whose summands are deterministic. For the second, count the number of $(1,k)$-PMDs whose probabilities are integer multiples of $\varepsilon$. So, for fixed $k$, our bound has the right qualitative dependence on $n$ (namely polynomial), and a near-right dependence on $1/\varepsilon$ (namely quasi-polynomial rather than polynomial). Moreover, it obtains the same qualitative dependence on $n$ and $\varepsilon$ as the $k=2$ cover of [DP09, DP14], namely polynomial in $n$ and quasi-polynomial in $1/\varepsilon$.

Learning PMDs

In view of tools for hypothesis selection from a cover (see, i.e., Theorem 7), our cover theorem directly implies that $(n,k)$-PMDs can be learned from $O(k^{5k}\cdot\log n\cdot\log^{k+2}(1/\varepsilon)/\varepsilon^{2})$ samples. These are near-optimal in terms of $\varepsilon$, as $\Omega(k/\varepsilon^{2})$ samples are necessary even for learning a $(1,k)$-PMD. We show that the dependence on $n$ can be completely removed from the learner, generalizing the results on Poisson Binomial Distributions [DDS12].

samples from $X$, runs in timeWe work in the standard “word RAM” model in which basic arithmetic operations on $O(\log n)$-bit integers are assumed to take constant time.

Additional Results: Learning k𝑘k-SIIRVs

A $(n,k)$-SIIRV is the sum of $n$ independent (single-dimensional) random variables supported on $\{0,\ldots,k-1\}$. SIIRVs generalize Poisson Binomial distributions, which correspond to the case $k=2$. At the same time, SIIRVs can be viewed as projections of PMDs onto the vector $(0,1,\ldots,k-1)$. In particular, if $X$ is a $(n,k)$-SIIRV, there exists a $(n,k)$-Poisson multinomial random vector $Y$, such that $X=(0,1,\ldots,k-1)^{\rm T}\cdot Y$.

1 Approach

The multi-dimensional nature of PMDs poses challenges in understanding their structure. The projection of a $(n,k)$-Poisson multinomial random vector onto each standard basis vector is a $n$-Poisson Binomial random variable, i.e. distributed as the sum of $n$ independent indicators. Depending on our choice of $\varepsilon$, the latter may be $\varepsilon$-close (in total variation distance) to a discretized Normal distribution (“heavy projection”) or a distribution whose essential support is a length $O(1/\varepsilon^{3})$ subinterval of $\{0,\ldots,n\}$ (“light projection”) [DP14]. Intuitively, one would like to aggregate all heavy projections into a discretized multi-dimensional Gaussian and all light projections into a distribution of small support, independent of $n$. However, projections onto different standard basis vectors may be correlated, and they cannot be disentangled this simply.

In fact, even if all projections of a PMD onto the standard basis vectors are heavy—even if they have variance super-polynomial in $k/\varepsilon$, it is still unclear whether the PMD can always be well approximated by a discretized multi-dimensional Gaussian. In particular, the multi-dimensional CLT of Valiant and Valiant [VV11] (Theorem 6) does pay a penalty that scales with $\log n$.

Finally, projections onto non-standard basis vectors may behave more erratically. As we pointed out earlier, the projection of a $(n,k)$-PMD onto the vector $\vec{v}=(0,1,\ldots,k-1)$ is a $(n,k)$-SIIRV, which need not be log-concave or even unimodal, and could even exhibit “mod-structure” and be $n$-modal; think of the distribution of $Y+2\cdot Z$ where $Z$ is sampled from a Binomial$(n,0.5)$ and $Y$ is a Bernoulli$(1/3)$. Whichever simpler distribution we identify to approximate a given $(n,k)$-PMD thus needs to respect the potential mod-structure that the PMD’s projection onto $\vec{v}$, its permutations or other integral vectors may exhibit.

Our analysis sidesteps the difficulties identified above by showing that, for all $\varepsilon$, $n$, $k$, a $(n,k)$-Poisson multinomial random vector is $\varepsilon$-close to the sum of a discretized Gaussian and an independent $({\rm poly}(k/\varepsilon),k)$-Poisson multinomial random vector. Roughly speaking, the Gaussian absorbs the variance in the heavy dimensions, and explains the correlation between light and heavy dimensions, while the sparse PMD explains the remaining variance in the light dimensions. Of course, what dimensions are “light” and “heavy” in the above discussion depends on our desired approximation $\varepsilon$.

At the heart of our proof lies the aforecited CLT by Valiant and Valiant [VV11], approximating a Poisson Multinomial by a discretized Gaussian. There are several issues with its application here: the accuracy of the approximation cannot be made an arbitrary $\varepsilon$, but worse, it deteriorates (logarithmically) as we increase $n$ or decrease the minimum eigenvalue of the covariance matrix of the PMD. The main intuition behind our structural theorem and the main technical roadblock for its proof lies in avoiding paying these two penalties.

To mitigate the latter cost (corresponding to the smallest eigenvalue), we use a stripped down version of the trickle-down sampling procedure from [DP08] to round the parameters of our given PMD. This allows us to shift the parameters of the PMD’s constituent random vectors such that they are either equal to or $1$, or sufficiently far from or $1$. A coordinated “rounding” of these parameters combined with a coupling argument and single-dimensional Poisson approximations allow us to argue that the effect of the rounding is small in the total variation distance of the resulting PMD compared to the original PMD. Each constituent random vector in the resulting PMD now has decent variance in every axis direction where it has non-zero variance. Partitioning the PMD’s constituent vectors into sets based on the axis directions where they have non-zero variance, we get that the minimum eigenvalue of each resulting sub-PMD is large in the span of these directions; see Proposition 6.Again, as pointed out earlier, when we refer to the eigenvalues of the covariance matrix of a PMD spanning a certain subspace, we always project the PMD onto a subspace of one dimension less, as otherwise the covariance matrix always has a eigenvalue since the distribution does not have full-dimensional support. Details about this step are given in Section B.1.

To avoid paying the logarithmic cost in the value of $n$ (the number of summands) which appears in the CLT, we repeatedly partition and sort the random vectors into buckets. The sub-PMD corresponding to each bucket will have the property that the logarithm of the number of summands is negligible compared to the minimum eigenvalue of its covariance matrix, so that we can apply the central limit theorem from [VV11]. We note that there will be a small number of random vectors which do not fall into a bucket that has this property – these leftover vectors result in the sparse Poisson Multinomial component in our structural result. Details about this step are given in Section B.2.

The above approximations result in a distribution comprising several discretized Gaussians and a sparse Poisson multinomial. We subsequently merge all component discretized Gaussians into a single distribution. It is well-known that the sum of two Gaussians is another Gaussian whose parameters are equal to the sum of the parameters of its two components. The same is not true for discretized Gaussians, and we must quantify the error induced by this merging operation. More details are provided in Section B.3.

Our structural results are described further in Section 3.

Cover

We provide two covers for $(n,k)$-PMDs, which are advantageous for different regimes of $k$ and $\varepsilon$. The first cover follows directly from Theorem 5, which gives a structural characterization of a PMD as the sum of an appropriately discretized Gaussian and a $(\operatorname*{poly}(k/\varepsilon),k)$-PMD. We simply take an additive grid over all the parameters of this characterization to achieve a cover size which is polynomial in $n$ and exponential in $k$ and $1/\varepsilon$.

Similar to [DP14], we can reduce the dependence of the cover size to pseudo-polynomial in $1/\varepsilon$, albeit at an increased cost in $k$. This is done using a generalization of the moment matching techniques known for Poisson Binomial distributions. At a high level, this avoids the naive gridding over all $(\operatorname*{poly}(k/\varepsilon),k)$-PMDs by filtering out the ones with unique “moment profiles,” which describe the first several moments of the distribution. We prove that any two distributions with matching moment profiles will have small total variation distance by leveraging results by Roos on Krawtchouk approximations to PMDs [Roo02].

A further description of our cover results is provided in Section 4.

Learning

Our cover theorem (Theorem 2) directly implies (using Theorem 7) that $(n,k)$-PMDs can be learned from $O(\log N/\varepsilon^{2})$ samples, where $N$ is the size of our cover. Given that $N$ is polynomial in $n$, the resulting sample complexity is logarithmic in $n$. To remove the dependence on $n$ from our sample complexity, we need to exploit not just the size but also the structure of the cover.

In particular, we know from our structural characterization (Theorem 1) that any $(n,k)$-Poisson Multinomial random vector is $\varepsilon$-close to the sum of a discretized multi-dimensional Gaussian and an independent $({\rm poly}(k/\varepsilon),k)$-PMD. The dependence of the cover size on $n$ is due to enumerating over a cover of discretized multi-dimensional Gaussians, as enumerating over $({\rm poly}(k/\varepsilon),k)$-PMDs has no dependence on $n$. The challenge is this: given sample access to an unknown $(n,k)$-PMD can we zoom in to a smaller set of candidate discretized multi-dimensional Gaussians whose size is independent of $n$ and which suffice for the purposes of guaranteeing an approximation to the unknown PMD?

Let us start with an easier task. Suppose that our structural theorem decides that a $(n,k)$-PMD is $\varepsilon$-close in total variation distance to a discretized multi-dimensional Gaussian. In this case, is it possible to recover the Gaussian from ${\rm poly}(k/\varepsilon)$ samples from the PMD? Intuitively the answer should be “yes,” as learning a multi-dimensional Gaussian to within $\varepsilon$ in total variation distance is feasible from $O(k/\varepsilon^{2})$ samples. Only there are two complications. First, we are seeking to actually learn a discretized multi-dimensional Gaussian and, most importantly, we do not have sample access to the Gaussian, but a distribution that is $\varepsilon$-close to it in total variation distance. The first complication becomes an issue when the covariance matrix of the Gaussian has minimum eigenvalue that does not scale with some ${\rm poly}(k/\varepsilon)$, which may very well be the case. The second is more severe as it necessitates robust estimators for the moments of a (discretized) multi-dimensional Gaussian that are resilient to an arbitrary movement of $\varepsilon$ probability mass. We are not aware of such estimators even for a (continuous) multi-dimensional Gaussian.

Despite these apparent issues, even in the simple case we are considering, the saving grace comes from a closer examination of the proof of our structural result. When our structural theorem deems a $(n,k)$-PMD approximable by a discretized multi-dimensional Gaussian, we can argue that the covariance matrices $\Sigma$ of the former and $\Sigma_{G}$ of the latter are spectrally close, satisfying $|x^{\rm\tiny T}\Sigma x-x^{\rm\small T}\Sigma_{G}x|\leq\varepsilon\cdot x^{\rm\small T}\Sigma x$, for all $x$. So it suffices to learn the covariance matrix of the PMD to which we have direct sample access, thereby obviating the need for a robust estimator. Learning the covariance matrix of a PMD is feasible from ${\rm poly}(k/\varepsilon)$ samples by bounding the kurtosis of any projection of the PMD (Lemma 8).

The bigger challenge is generalizing the approach to when our structural theorem deems a $(n,k)$-Poisson Multinomial random vector $X$ approximable by the sum of a discretized multi-dimensional Gaussian $G$ and a $({\rm poly}(k/\varepsilon),k)$-Poisson Multinomial random vector $Y$. We can enumerate over the latter, but enumerating over the former is too expensive (i.e. will incur a dependence on $n$). So we have to learn it with sample access to $X$. Unfortunately, our spectral approximation is now much weaker. The covariance matrices $\Sigma$ of $X$ and $\Sigma_{G}$ of $G$ are now related as follows, for all $x$: $|x^{\rm\tiny T}\Sigma x-x^{\rm\small T}\Sigma_{G}x|\leq\varepsilon\cdot x^{\rm\small T}\Sigma x+{\rm poly}(k/\epsilon)$. Hence, for directions $x$ where the variance $x^{\rm\tiny T}\Sigma x$ of $X$ is small, this approximation is quite loose to just approximate $\Sigma_{G}$ with $\Sigma$.

Our approach is instead to use samples from $X$ to get a handle on the spectrum of $\Sigma_{G}$. As before, by bounding the kurtosis of any projection of the PMD, we can produce an estimate $\hat{\Sigma}$ that approximates $\Sigma$ spectrally: for all $x$, $|x^{\rm\tiny T}\Sigma x-x^{\rm\small T}\hat{\Sigma}x|\leq\varepsilon\cdot x^{\rm\small T}\Sigma x$ (Lemma 8). Then, using Courant minimax principle through the proof of our structural result, we can argue that the $i$-th eigenvalue $\lambda^{G}_{i}$ of $\Sigma_{G}$ and $\hat{\lambda}_{i}$ of $\hat{\Sigma}$ are related as follows: $|\lambda^{G}_{i}-\hat{\lambda}_{i}|\leq O(\varepsilon)\hat{\lambda}_{i}+{\rm poly}(k/\epsilon)$. So, using the eigenvalues of our learned $\hat{\Sigma}$, we can produce a small cover for the eigenvalues of $\Sigma_{G}$. Unfortunately, the corresponding eigenvectors of $\Sigma_{G}$ and $\hat{\Sigma}$ need not be as closely related, and it is not clear how to grid over those as the ratio of the smallest to the largest eigenvalue may be polynomial in $n$. We show how to use the knowledge of the eigenvalues and the spectral relation between $\hat{\Sigma}$ and $\Sigma_{G}$ to produce a small cover over matrices $\hat{\Sigma}_{G}$ (and not eigenvectors) such that at least one matrix in the cover spectrally approximates our target $\Sigma_{G}$. The details are provided in Section D.3. At this point, we have a small cover over possible distributions $Y$ and a small cover over possible discretized multi-dimensional Gaussians. So we can select among these hypotheses using Theorem 7.

Our learning algorithm is described in Section 5.

Preliminaries

Throughout this paper, we will repeatedly refer to three key parameters, $c=c(\varepsilon,k)=\operatorname*{poly}(\varepsilon/k)$, $t=t(\varepsilon,k)=\operatorname*{poly}(k/\varepsilon)$, and $\gamma=O(1)$. We set

for constants $\delta_{c},\delta_{t},\delta_{\gamma}>0$.

2 Definitions

We start by defining several of the distribution classes we will consider. First, and most importantly, we start with a formal definition of Poisson Multinomial Distributions.

A $k$-Categorical Random Variable ($k$-CRV) is a random variable that takes values in $\{e_{1},\dots,e_{k}\}$ where $e_{j}$ is the $k$-dimensional unit vector along direction $j$. $\pi(i)$ is the probability of observing $e_{i}$.

An $(n,k)$-Poisson Multinomial Distribution ($(n,k)$-PMD) is given by the law of the sum of $n$ independent but not necessarily identical $k$-CRVs. An $(n,k)$-PMD is parameterized by a nonnegative matrix $\pi\in^{n\times k}$ each of whose rows sum to $1$ is denoted by $M^{\pi}$, and is defined by the following random process: for each row $\pi(i,\cdot)$ of matrix $\pi$ interpret it as a probability distribution over the columns of $\pi$ and draw a column index from this distribution. Finally, return a row vector recording the total number of samples falling into each column (the histogram of the samples).

We note that a sample from an $(n,k)$-PMD is redundant – given $k-1$ coordinates of a sample, we can recover the final coordinate by noting that the sum of all $k$ coordinates is $n$. For instance, while a Binomial distribution is over a support of size $2$, a sample is $1$-dimensional since the frequency of the other coordinate may be inferred given the parameter $n$. With this inspiration in mind, we define the Generalized Multinomial Distribution, which is the primary object of study in [VV11].

A Truncated $k$-Categorical Random Variable is a random variable that takes values in $\{0,e_{1},\dots,e_{k-1}\}$ where $e_{j}$ is the $(k-1)$-dimensional unit vector along direction $j$, and is the $(k-1)$ dimensional zero vector. $\rho(0)$ is the probability of observing the zero vector, and $\rho(i)$ is the probability of observing $e_{i}$.

An $(n,k)$-Generalized Multinomial Distribution ($(n,k)$-GMD) is given by the law of the sum of $n$ independent but not necessarily identical truncated $k$-CRVs. A GMD is parameterized by a nonnegative matrix $\rho\in^{n\times(k-1)}$ each of whose rows sum to at most $1$ is denoted by $G^{\rho}$, and is defined by the following random process: for each row $\rho(i,\cdot)$ of matrix $\rho$ interpret it as a probability distribution over the columns of $\rho$ – including, if $\sum_{j=1}^{k}\rho(i,j)<1$, an “invisible” column – and draw a column index from this distribution. Finally, return a row vector recording the total number of samples falling into each column (the histogram of the samples).

For both $(n,k)$-PMDs and $(n,k)$-GMDs, we will refer to $n$ and $k$ as the size and dimension, respectively.

We note that a PMD corresponds to a GMD where the “invisible” column is the zero vector, and thus the definition of GMDs is more general than that of PMDs. However, whenever we refer to a GMD in this paper, it will explicitly have a non-zero invisible column.

While we will approximate the Multinomial distribution with Gaussian distributions, it does not make sense to compare discrete distributions with continuous distributions, since the total variation distance is always $1$. As such, we must discretize the Gaussian distributions. We will use the notation $\lfloor x\rceil$ to say that $x$ is rounded to the nearest integer (with ties being broken arbitrarily). If $x$ is a vector, we round each coordinate independently to the nearest integer.

As seen in the definition of an $(n,k)$-GMD, we have one coordinate which is equal to $n$ minus the sum of the other coordinates. We define a similar notion for a discretized Gaussian. However, we go one step further, to take care of when there are several such Gaussians which live in disjoint dimensions. By this, we mean that given two Gaussians, the set of directions in which they have a non-zero variance are disjoint. Without loss of generality (because we can simply relabel the dimensions), we assume all of a Gaussian’s non-zero variance directions are consecutive, i.e., the covariance matrix is all zeros, except for a single block on the diagonal. Therefore, when we add the covariance matrices, the result is block diagonal. The resulting distribution is described in the following definition.

The structure preserving rounding of a multidimensional Gaussian Distribution takes as input a multi-dimensional Gaussian $\mathcal{N}(\mu,\Sigma)$ with $\Sigma$ in block-diagonal form. It chooses one coordinate as a “pivot” in each block, samples from the Gaussian ignoring these pivots and rounds each value to the nearest integer. Finally, the pivot coordinate of each block is set by taking the difference between the sum of the means and the sum of the values sampled within the block.

Structure of PMDs

In this section, we show a structural result, stating that any $(n,k)$-PMD is close to the sum of an appropriately discretized Gaussian and a $(\operatorname*{poly}(k/\varepsilon),k)$-PMD:

For parameters $c$ and $t$ as described in Section 2.1, every $(n,k)$-Poisson multinomial random vector is $\varepsilon$-close to the sum of a Gaussian with a structure preserving rounding and a $(tk^{2},k)$-Poisson multinomial random vector. For each block of the Gaussian, the minimum non-zero eigenvalue of $\Sigma_{i}$ is at least $\frac{tc}{2k^{4}}$.

There are three main steps in the proof of this theorem.

First, we replace our $(n,k)$-PMD with one where all parameters are sufficiently far from and $1$, while still being close to the original in total variation distance. To motivate this operation, we introduce one of our main tools in our approach, the central limit theorem of Valiant and Valiant [VV11], which approximates an $(n,k)$-GMD by a discretized multivariate Gaussian.

Given a generalized multinomial distribution $G^{\rho}$, with $k$ dimensions and $n$ rows, let $\mu$ denote its mean and $\Sigma$ denote its covariance matrix, then

where $\sigma^{2}$ is the minimum eigenvalue of $\Sigma$.

We note that this has an error term which depends on the minimum eigenvalue of the covariance matrix of the GMD. If we perform this rounding procedure and ignore any zero coordinates, then we are given the guarantee that the minimum eigenvalue will be sufficiently large.

Recall that in Section 2.1 we have set $c=\operatorname*{poly}(\varepsilon/k)$. This lemma summarizes the result of the rounding procedure:

For any $c\leq\frac{1}{2k}$, given access to the parameter matrix $\rho$ for an $(n,k)$-PMD $M^{\rho}$, we can efficiently construct another $(n,k)$-PMD $M^{\hat{\rho}}$, such that, for all $i,j$, $\hat{\rho}(i,j)\not\in(0,c)$, and

The procedure starts by fixing two coordinates $i$ and $j$, and considers all CRVs with a parameter in $i$ which is close to , and has maximum parameter in coordinate $j$. We move some of the weight in this “heavy” coordinate either to or from the “light” coordinate, while approximately preserving the overall mean vector of the set of CRVs.

The analysis of this process uses a stripped-down version of the “trickle-down” process in [DP08]. This gives an approximate way to sample from a PMD, resulting in a distribution which is very close in total variation distance. While we postpone technical details to Section B.1, roughly speaking, it works as follows. First, take a sample from the PMD but disregard the values for its light coordinate $i$ and heavy coordinate $j$. Instead, sample a new value for coordinate $i$ according to a Poisson distribution with parameter $\mu_{i}$, the mean value for coordinate $i$. Finally, set coordinate $j$ to ensure that all coordinates of the sample sum to $n$. As mentioned before, the rounding process approximately preserves the value of $\mu_{i}$, and thus this alternate sampling procedure is closely coupled for the rounded and original PMD. Thus, by triangle inequality, the rounded and original PMDs are close in total variation distance.

We repeat this rounding procedure for each $i$ and $j$, eventually leading to all parameters either being equal to or far from and $1$. A full description and analysis of the rounding procedure are in Section B.1.

Now, we have a “massaged” $(n,k)$-PMD $M^{\hat{\rho}}$, with no parameters lying in the intervals $(0,c)$ or $(1-c,1)$. Next, we will show how to relate the massaged $(n,k)$-Poisson multinomial random vector to a sum of $k$ Gaussians with a structure preserving rounding plus a “sparse” $(\operatorname*{poly}(k/\varepsilon),k)$-PMD. The general roadmap is as follows. We start by partitioning the constituent $k$-CRVs into $k$ sets, $S_{1},\dots,S_{k}$, based on which basis vector we are most likely to observe. We work seperately for each set $S_{i}$ by considering the GMD formed by leaving out the coordinate $i$. Our goal is to use the CLT of Theorem 6 to bound the total variation distance between the corresponding GMD and a discretized Gaussian with the same mean and covariance matrix. We must be careful when applying Theorem 6, since the bound depends on the size of the GMD. Instead of applying the theorem directly, to get a useful bound, we further partition the set $S_{i}$ into smaller subsets and apply the theorem to each of the resulting subsets. We can then “merge” the resulting discretized Gaussians together using the following lemma whose proof is given in Section A.5:

Let $X_{1}\sim\mathcal{N}(\mu_{1},\Sigma_{1})$ and $X_{2}\sim\mathcal{N}(\mu_{2},\Sigma_{2})$ be $k$-dimensional Gaussian random variables, and let $\sigma=\min_{j}\max_{i}\sigma_{i,j}$ where $\sigma_{i,j}$ is the standard deviation of $X_{i}$ in the direction parallel to the $j$th coordinate axis. Then

In more detail, we partition each set $S_{i}$ into $2^{k-1}$ subsets, grouping together CRVs according to the dimensions they are non-zero in, i.e. set $S_{i}^{\mathcal{I}}$ contains all CRVs that are non zero in the coordinates given by set $\mathcal{I}\subseteq[k]\setminus\{i\}$. We then group these sets into buckets, where a set is assigned to a bucket depending on its cardinality; bucket $B^{l}$ gets all sets $S_{i}^{\mathcal{I}}$ with $|S_{i}^{\mathcal{I}}|\in[l^{\gamma}t,(l+1)^{\gamma}t)$, with $\gamma=O(1)$ and $t=\operatorname*{poly}(k/\varepsilon)$ as defined in Section 2.1. This bounds the ratio between the size and the minimum eigenvalue of the covariance of the GMD within every bucket other than $B^{0}$. This allows us to apply Theorem 6 and replace the CRVs within each bucket $B^{l}$ for $l\geq 1$ with a discretized Gaussian, leaving us with a $(\operatorname*{poly}(2^{k}/\varepsilon),k)$-GMD consisting of all the CRVs of bucket $B^{0}$. To reduce the number of remaining CRVs to polynomial in $k$, we show that by removing only $\operatorname*{poly}(k/\varepsilon)$ of these CRVs, we can apply Theorem 6 again to the rest and obtain another discretized Gaussian. In particular, in Section B.2, we prove the following lemma:

Let $G^{\hat{\rho}_{k}^{0}}$ be the $(|B^{0}|,k)$-GMD induced by the truncated CRVs in bucket $B^{0}$. Given $\hat{\rho}_{k}^{0}$, we can efficiently compute a partition of $B^{0}$ into $S$ and $\bar{S}$, where $|\bar{S}|\leq kt$. Letting $\mu_{S}$ and $\Sigma_{S}$ be the mean and covariance matrix of the $(|S|,k)$-GMD induced by $S$, and $G^{\hat{\rho}_{k}^{\bar{S}}}$ be the $(|\bar{S}|,k)$-GMD induced by $\bar{S}$,

Furthermore, the minimum non-zero eigenvalue of $\Sigma_{S}$ is at least $\frac{tc}{k}$.

After merging together all discretized Gaussians (at most one coming from each bucket $B^{l}$ for all $l\geq 0$) by iteratively applying Lemma 2, we are able to approximate each original set of CRVs $S_{i}$ as the sum of a single discretized Gaussian and a $(\operatorname*{poly}(k/\varepsilon),k)$-PMD. Combining the result from each of the sets $S_{i}$ of the initial partition, we obtain the sum of $k$ discretized Gaussians and a $(\operatorname*{poly}(k/\varepsilon),k)$-PMD. The details of this step are described in Section B.2.

The final step is to show that the $k$ discretized Gaussians can be merged into a single Gaussian with a structure preserving rounding. We note that we cannot apply Lemma 2 here, since each discretized Gaussian has a different pivot coordinate that has been left out. (Recall that by construction, the CRVs in set $S_{i}$ are approximated by a discretized Gaussian that leaves out coordinate $i$). We thus need a new tool to enable us to merge Gaussians defined in different dimensions. The main idea is that if two Gaussians with a structure preserving rounding overlap in some dimension, we can use the common dimension as the pivot. We then add the mean vectors and covariance matrices to merge the distributions. Iteratively repeating this process will merge all distributions which overlap in some coordinate. This leaves us with one or many discretized Gaussians that lie in completely disjoint coordinates which we can describe as a single Gaussian with a structure preserving rounding (defining blocks according to the coordinates spanned by each Gaussian). If these were (continuous) Gaussians, the swapping and merging operations would have no cost, but some care is required when dealing with discretized Gaussians. There are two costs which we must bound here. First, we must show that swapping the pivot of a PMD is inexpensive, and second, we need to bound the cost of repeatedly merging Gaussians.

We bound the cost of swapping the pivot by proving the following lemma:

$Y_{i}$ be the distribution in which we draw a sample $(x_{1},\dots,x_{k-1})\sim X_{i}$ and return

By applying Lemma 4, we can make two discretized Gaussians have the same left out coordinate and then merge them using Lemma 2 if at least one of them has large variance in every direction. While each of the $k$ discretized Gaussians starts with this property (for the dimensions in which it is non-deterministic), it is not clear whether this is true after a sequence of pivot swaps and merges.

In many cases, swapping the pivot decreases the minimum eigenvalue of the distribution’s covariance matrix by a factor of $\operatorname*{poly}(k)$. This is acceptable if we only perform a single swap, but naively applying this bound for a sequence of $k$ swaps and merges results in the minimum eigenvalue dropping by a factor of $k^{O(k)}$. We show that such a bad situation cannot occur, no matter how one performs the sequence of swaps and merges, by proving the following lemma:

$\Sigma^{(i)}$ has eigenvalue with corresponding eigenvector $\vec{1}$

There exists coordinate $j^{*}\in S^{(i)}$ such that $\Sigma^{(i)}_{S^{(i)}\setminus\{j^{*}\}}$ has minimum eigenvalue at least $\lambda$

Then, for all $j\in S$, the minimum eigenvalue of $\Sigma_{S\setminus\{j\}}$ is at least $\frac{\lambda}{2k^{3}}$.

The details of this step, the proofs of Lemma 4 and Lemma 5 as well as the proof of Theorem 5 are described in Section B.3.

Covers for PMDs

In this section, we describe a pair of covers for $(n,k)$-PMDs.

The first cover follows directly from Theorem 5, which gives a structural characterization of a $(n,k)$-Poisson multinomial random vector as the sum of an appropriately discretized Gaussian and an $(tk^{2},k)$-Poisson multinomial random vector. We grid over all possible mean vectors and covariance matrices for the Gaussian component, and all possible parameter values for the $(tk^{2},k)$-PMD. These are covered by sets of size $(n\cdot\operatorname*{poly}(k/\varepsilon))^{k^{2}}$ and $2^{\operatorname*{poly}(k/\varepsilon)}$ respectively, resulting in an overall cover of size $n^{k^{2}}\cdot 2^{\operatorname*{poly}(k/\varepsilon)}$.

The proof of this lemma is presented in Section C.1.

The second cover further sparsifies the cover for the $(tk^{2},k)$-PMD component, by using a multivariate generalization of the moment matching technique described in [DP14]. This reduces the cover size for this component to $2^{O(k^{5k}\log^{k+2}(1/\varepsilon))}$. In [Roo02], Roos shows that a PMD can be written as the weighted sum of partial derivatives of a regular multinomial distribution. He goes on to show that dropping the higher order derivatives in this sum results in a total variation approximation, where the quality of the approximation depends on the parameters of the PMD and the point at which we evaluate the derivatives. We take advantage of this tool to obtain an $\varepsilon$-approximation, through a careful partitioning of the CRVs and choice of point at which to evaluate the derivatives of the multinomial distributions. This implies that any two distributions which have matching “moment profiles” (which roughly describe the lower order derivatives of the distribution) are $\varepsilon$-close to each other, and thus only one representative element must be kept from each such equivalence class. The size of the cover follows by a counting argument on the number of moment profiles.

The proof of this lemma is given in Section C.2. We note that this cover can be efficiently enumerated over, using a dynamic program similar to that of [DP14].

By combining these two lemmas, we obtain Theorem 2.

Learning PMDs

As mentioned before, Theorem 2 combined with Theorem 7 below (taken from [DK14]) immediately implies that $(n,k)$-PMDs can be learned from $O(\log N/\varepsilon^{2})$ samples, where $N$ is the size of our cover.

Theorem 7 is using a tournament-style algorithm for hypothesis selection, which takes a set of candidate distributions and outputs one which is $O(\varepsilon)$-close to the unknown distribution (if such a distribution exists)We note that this tournament additionally requires a “PDF comparator,” which we describe for our setting in Section D.4. . Given that $N$ is polynomial in $n$, the resulting sample complexity is logarithmic in $n$. To remove the dependence on $n$ from our sample complexity, we need to exploit not just the size but also the Gaussian structure of the cover. Instead of trying all possible Gaussians that the cover could describe, we instead estimate the moments of the Gaussian directly.

Our strategy will not be to generate an $\varepsilon$-cover for all $(n,k)$-PMDs, but instead we take samples and select only distributions from our cover which are consistent with the data. Similar to before, we will apply Theorem 7 to do hypothesis selection but instead of applying it to the complete cover resulting from Theorem 2, we will apply it to a much smaller set of hypothesis that we obtain after making several “guesses” for the parameters of our distribution. At least one set of these parameters will be sufficiently accurate to obtain an $\varepsilon$ total variation distance guarantee and we will be able to determine a good candidate using Theorem 7.

The first step of our learning algorithm is to guess the block-diagonal structure of the Gaussian component of our distribution by guessing the partition of the coordinates and choosing an arbitrary pivot within each block. This requires at most $k^{k}$ guesses. Note that any choice of pivot in the partition is acceptable (as shown in Lemma 4 above).

Next, we estimate the mean and covariance of the Gaussian component for each block. We need to estimate them accurately enough in order to learn each block of the discretized Gaussians to within $O(\varepsilon/k)$ in total variation distance. A useful tool for showing this is the following proposition:

Proposition 1 implies that, in order to achieve the required bound in total variation distance, it suffices to get an estimate that approximately matches the mean and variance of the Gaussian component in every direction. In Section D.1, we prove Lemma 8 which shows that using $\operatorname*{poly}(k)/\varepsilon^{2}$ samples from the PMD, we can get an estimate of the mean and covariance matrix that achieves this guarantee in every direction. However, this estimate is with respect to the PMD we are sampling from and not with respect to the Gaussian component, which is the guarantee we desire.

Given sample access to a $(n,k)$-PMD $X$ with mean $\mu$ and covariance matrix $\Sigma$ (with minimum eigenvalue at least $1$), there exists an algorithm which can produce estimates $\hat{\mu}$ and $\hat{\Sigma}$ such that with probability at least $9/10$:

The sample and time complexity are $O(k^{4}/\varepsilon^{2})$.

In order to obtain a guarantee for the Gaussian component, we observe that there are two possible sources of errors in our estimation:

The first source of error comes from the rounding step. In proving our structural result, the real PMD had to be rounded so that no CRV has any probability that is in the range $(0,c)$, which affected the mean and covariance. In Section D.2, we show that this only affects the mean and variance in each direction up to a small multiplicative factor.

The second source of error is due to the existence of the sparse component creates an additional additive error in each direction. This error might be very significant in some directions as the variance of the Gaussian component can be very small compared to the number of sparse CRVs.

Understanding that our estimation is off by an additive error and a multiplicative error, we show how to efficiently correct this estimation by searching around it for the underlying covariance matrice of the Gaussian distribution. In particular, we obtain a cover of positive semidefinite matrices that are close to the estimated covariance matrix and which contains a good approximation to the covariance matrix of the underlying Gaussian. This is challenging because the above two sources of error might affect the spectrum of the covariance matrix significantly. However, we are able to tackle this issue by carefully guessing appropriate corrections to the eigenvectors and eigenvalues of the matrix. We prove Lemma 9 which states that this cover has cardinality at most $(k/\varepsilon)^{O(k^{2})}$, and thus we can get a very accurate estimate for the underlying Gaussian distribution by guessing different points in the cover.

Let $A$ be a symmetric $k\times k$ PSD matrix with minimum eigenvalue $1$ and let $S$ be the set of all matrices $B$ such that $|y^{T}(A-B)y|\leq\varepsilon_{1}y^{T}Ay+\varepsilon_{2}y^{T}y$ for all vectors $y$, where $\varepsilon_{1}\in[0,1/4)$ and $\varepsilon_{2}\in[0,\infty)$. Then, there exists an $\varepsilon$-cover $S_{\varepsilon}$ of $S$ that has size $|S_{\varepsilon}|\leq\left(\frac{k(1+\varepsilon_{2})}{\varepsilon}\right)^{O(k^{2})}$.

At this point, we have a collection of distributions such that at least one is close to the Gaussian component. We do the same for the sparse PMD component by simply enumerating over all the elements in the cover. By reading the corresponding term from the statement of Theorem 2, this requires $\min\{2^{\operatorname*{poly}(k/\varepsilon)},2^{O(k^{5k}\cdot\log^{k+2}(1/\varepsilon))}\}$ guesses.

In conclusion, using $\operatorname*{poly}(k)/\varepsilon^{2}$ samples, we have generated a set $\mathcal{S}$ of size

which contains a distribution which is $\varepsilon$-close to the true distribution with constant probability. In order to choose a “good” distribution from this set, we apply the hypothesis selection algorithm of Theorem 7 to obtain a distribution which is $O(\varepsilon)$-close to the unknown distribution with constant probability, which concludes the proof of Theorem 3. More details about the learning steps and complete proofs can be found in Section D.

Learning k𝑘k-SIIRVs

The proof uses the structural result of [DDO+13], which says that any $(n,k)$-SIIRV is close to either a low variance distribution with limited support, or a high variance distribution which enjoys certain Gaussian structural properties.

Let $S=X_{1}+\dots+X_{n}$ be a $(n,k)$-SIIRV for some positive integer $k$. Let $\mu$ and $\sigma^{2}$ be respectively the mean and variance of $S$. Then for all $\varepsilon>0$, the distribution of $S$ is $O(\varepsilon)$-close in total variation distance to one of the following:

a random variable supported on $\frac{k^{9}}{\varepsilon^{4}}$ consecutive integers with variance $\sigma^{2}\leq 15(k^{18}/\varepsilon^{6})\log^{2}(1/\varepsilon)$; or

the sum of two independent random variables $S_{1}+cS_{2}$, where $c$ is some positive integer $1\leq c\leq k-1$, $S_{2}$ is distributed according to $\lfloor\mathcal{N}(\mu,\sigma^{2})\rceil$, and $S_{1}$ is a $c$-IRV; in this case, $\sigma^{2}=\Omega\left(\frac{k^{18}}{\varepsilon^{6}}\log^{2}(1/\varepsilon)\right)$.

To cover the former case, we use the PMD cover of Theorem 2. In this setting, the SIIRV has a variance upper bounded by $\operatorname*{poly}(k/\varepsilon)$. By applying a rounding procedure, it can be shown that this can be approximated by an offset $(\operatorname*{poly}(k/\varepsilon),k)$-SIIRV. Recalling that any $(n,k)$-SIIRV can be expressed as the projection of an $(n,k)$-PMD onto the vector $(0,1,\dots,k-1)$ and applying our quasi-polynomial cover result in Theorem 2 covers this case with $2^{O(k^{5k}\cdot\log^{k+2}(1/\varepsilon))}$ candidates.

To cover the latter case, we first perform $k-1$ guesses for the value of $c\in[k-1]$. For each guess, we learn the two distributions $S_{1}$ and $S_{2}$ separately. To learn $S_{1}$, we use the same approach as [DDO+13], which uses the empirical distribution obtained after mapping the samples onto $\{0,1,\dots,c-1\}$ using their residue mod $c$. Our method for learning $S_{2}$ is novel – we first round the value of each sample down to the next multiple of $c$, and examine the distribution on this support, which will be close in total variation distance to $S_{2}$. We estimate the moments of this distribution using robust statistical tools, as in [DDO+13]. The empirical median is used to estimate the mean, and a rescaling of the interquartile range is used to estimate the standard deviation. Thus, we cover this case using only $k-1$ candidates, one for each guess of $c$.

References

Appendix A Useful Tools

To compare probability distributions, we will require the total variation and Kolmogorov distances:

The total variation distance between two probability measures $P$ and $Q$ on a $\sigma$-algebra $F$ is defined by

Unless explicitly stated otherwise, in this paper, when two distributions are said to be $\varepsilon$-close, we mean in total variation distance.

The Kolmogorov distance between two probability measures $P$ and $Q$ with CDFs $F_{P}$ and $F_{Q}$ is defined by

We note that Kolmogorov distance is, in general, weaker than total variation distance. In particular, total variation distance between two distributions is lower bounded by the Kolmogorov distance.

A.2 Probabilistic Tools

We will use the following form of Chernoff/Hoeffding bounds:

Let $Z_{1},\dots,Z_{m}$ be independent random variables with $Z_{i}\in$ for all $i$. Then, if $Z=\sum_{i=1}^{n}Z_{i}$ and $\gamma\in(0,1)$,

We note the Dvoretzky-Kiefer-Wolfowitz (DKW) inequality, which is a powerful tool, giving a generic algorithm for learning any distribution with respect to the Kolmogorov metric [DKW56].

We will use the Data Processing Inequality for total variation distance (see part (iv) of Lemma 2 of [Rey11] for the proof). This lemma says that taking any function of two random variables can only reduce their total variation distance. Our statement of the inequality is taken from [DDO+13].

Let $X,X^{\prime}$ be two random variables over a domain $\Omega$. Fix any (possibly randomized) function $F$ on $\Omega$ (which may be viewed as a distribution over deterministic functions on $\Omega$) and let $F(X)$ be the random variable such that a draw from $F(X)$ is obtained by drawing independently $x$ from $X$ and $f$ from $F$ and then outputting $f(x)$ (likewise for $F(X^{\prime})$). Then we have

Finally, we require a hypothesis selection algorithm. Roughly, given a set of $N$ distributions with the guarantee that at least one is $\varepsilon$-close to an unknown distribution $X$, we can choose a hypothesis which is $O(\varepsilon)$-close to $X$. The running time is near-linear in $N$ and the number of samples is logarithmic in $N$.

Let $H_{1}$ and $H_{2}$ be probability distributions over some set ${\cal D}$. A PDF comparator for $H_{1},H_{2}$ is an oracle that takes as input some $x\in{\cal D}$ and outputs $1$ if $H_{1}(x)>H_{2}(x)$, and otherwise.

A.3 Bounds for Distances Between Distributions

Given two $k$-dimensional Gaussians $\mathcal{N}_{1}=\mathcal{N}(\mu_{1},\Sigma_{1}),\mathcal{N}_{2}=\mathcal{N}(\mu_{2},\Sigma_{2})$ such that for all $i,j\in[k]$, $|\Sigma_{1}(i,j)-\Sigma_{2}(i,j)|\leq\alpha$, and the minimum eigenvalue of $\Sigma_{1}$ is at least $\sigma^{2}$,

Without loss of generality, assume $\mu=0$, $\Sigma=I$, and $\Sigma^{\prime}$ is diagonal. This can be done by setting $y=Q\Lambda^{-1/2}Q^{\prime T}x$, where $\Sigma=Q\Lambda Q^{T}$ and $\Sigma^{\prime}=Q^{\prime}\Lambda^{\prime}Q^{\prime T}$ are the eigendecompositions of $\Sigma$ and $\Sigma^{\prime}$.

This implies that we now have the following guarantees for all $i\in[k]$:

Since each coordinate is independent and noting that $\Sigma^{\prime}_{i,i}\geq 1-\varepsilon$, we can apply Proposition 2 to each coordinate direction to obtain a total variation distance of $2\varepsilon k$. ∎

Let $X_{1},\dots,X_{n}$ be independent random variables, with $E[X_{i}]=0,E[X_{i}^{2}]=\sigma_{i}^{2}>0,E[|X_{i}|^{3}]=\rho_{i}<\infty$, and define $X=\sum_{i=1}^{n}X_{i},\sigma^{2}=\sum_{i=1}^{n}\sigma_{i}^{2},\rho=\sum_{i=1}^{n}\rho_{i}$. Then for an absolute constant $C_{0}\leq 0.56$,

A.4 Covariance Matrices of Truncated Categorical Random Variables

First, recall the definition of a symmetric diagonally dominant matrix.

A matrix $A$ is symmetric diagonally dominant (SDD) if $A^{T}=A$ and $A_{ii}\geq\sum_{j\neq i}|A_{ij}|$ for all $i$.

As a tool, we will use this corollary of the Gershgorin Circle Theorem [Ger31] which follows since all eigenvalues of a symmetric matrix are real.

Given an SDD matrix $A$ with positive diagonal entries, the minimum eigenvalue of $A$ is at least $\min_{i}A_{ii}-\sum_{j\neq i}|A_{ij}|$.

The minimum eigenvalue of the covariance matrix $\Sigma$ of a truncated CRV is at least $\rho(0)\min_{i}\rho(i)$.

We note that $\Sigma$ is SDD, since $\sum_{j\neq i}|\Sigma_{ij}|=\rho(i)\sum_{j\neq i}\rho(j)=\rho(i)(1-\rho(i)-\rho(0))\leq\rho(i)(1-\rho(i))=\Sigma_{ii}$. Thus, applying Proposition 5, we see that the minimum eigenvalue of $\Sigma$ is at least $\min_{i}\rho(i)(1-\rho(i))-\rho(i)(1-\rho(i)-\rho(0))=\rho(0)\min_{i}\rho(i)$. ∎

A.5 Sums of Discretized Gaussians

In this section, we will obtain total variation distance bounds on merging the sum of discretized Gaussians. It is well known that the sum of multiple Gaussians has the same distribution as a single Gaussian with parameters equal to the sum of the components’ parameters. However, this is not true if we are summing discretized Gaussians – we quantify the amount we lose by replacing the distribution with a single Gaussian, and then discretizing afterwards.

As a tool, we will use the following result from [DDO+13]:

Let $X_{1}\sim\mathcal{N}(\mu_{1},\sigma_{1}^{2})$ and $X_{2}\sim\mathcal{N}(\mu_{2},\sigma_{2}^{2})$. Then

First, suppose without loss of generality that $\sigma_{1}\geq\sigma_{2}$.

The second inequality uses Proposition 7. ∎

The proof is by induction on $k$. The base case of $k=1$ is handled by Proposition 8. For general $k$, we use a standard hybridization argument. Denote the $j$th coordinate of $X_{i}$ as $x_{ij}$.

The first inequality is the triangle inequality, the second uses Lemma 13, and the third uses the induction hypothesis and Proposition 8. ∎

Appendix B Details from Section 3

Fix some coordinate $x$, and select all $k$-CRVs where the parameter in coordinate $x$ is in the range $(0,c)$. Partition this subset into $k-1$ sets, depending on which coordinate $y\neq x$ is the heaviest. We apply a rounding procedure separately to each of these sets. After this procedure, none of the parameters in coordinate $x$ will be in $(0,c)$. We repeat this for all $k$ possible settings of $x$. From the description below (and the restriction that $c\leq\frac{1}{2k}$), it will be clear that we will not “undo” any of our work and move probabilities back into $(0,c)$, so $O(k^{2})$ applications of our rounding procedure will produce the result claimed in the theorem statement.

Recall that the goal of this rounding procedure will be to shift probability mass either to or from coordinate $x$ to coordinate $y$, such that no parameter in coordinate $x$ lies in the interval $(0,c)$, while simultaneously approximately preserving the mean vector of the distribution. We are able to do this since coordinate $y$ is “heavy” and thus small additions will not affect the distribution in this coordinate much.

We fix some $x,y$ in order to describe and analyze the process more formally. Define $\mathcal{I}^{x}_{y}=\{i\,|\,0<\rho(i,x)<c\wedge y=\arg\max_{j}\rho(i,j)\}$ (breaking ties lexicographically), and let $M^{\rho_{I^{x}_{y}}}$ be the $(n,k)$-PMD induced by this set. For the remainder of this section, without loss of generality, assume that the indices selected by $\mathcal{I}^{x}_{y}$ are $1$ through $|\mathcal{I}^{x}_{y}|$.

Select an arbitrary set $\mathcal{R}\subseteq\mathcal{I}^{x}_{y}$ such that $|\mathcal{R}|=\left\lfloor\frac{\sum_{i^{\prime}\in I^{x}_{y}}\rho_{I^{x}_{y}}(i^{\prime},x)}{c}\right\rfloor$. Intuitively, this set will be the CRVs for which we set the parameter $\rho(\cdot,x)$ to be $c$, while $\mathcal{I}^{x}_{y}\setminus\mathcal{R}$ will have $\rho(\cdot,x)$ set to . We can perform the following rounding scheme to $\rho_{I^{x}_{y}}$ to obtain a new parameter matrix $\hat{\rho}_{I^{x}_{y}}$:

We define the process Fork, for sampling from a $k$-CRV $\rho(i,\cdot)$ in $\mathcal{I}^{x}_{y}$:

Let $X_{i}$ be an indicator random variable, taking $1$ with probability $\frac{1}{k}$ and otherwise.

If $X_{i}=1$, then return $e_{x}$ with probability $k\rho(i,x)$ and $e_{y}$ with probability $1-k\rho(i,x)$.

If $X_{i}=0$, then return $e_{j}$ with probability if $j=x$, $\frac{k}{k-1}(\rho(i,x)+\rho(i,y)-\frac{1}{k})$ if $j=y$, and $\frac{k}{k-1}\rho(i,j)$ otherwise.

The intuition behind this procedure is that we isolate the changes in our rounding procedure when $X_{i}=1$, as when $X_{i}=0$, the rounded and unrounded distributions are identical. We note that Fork is well defined as long as $\rho(i,x)\leq\frac{1}{k}$ and $\rho(i,x)+\rho(i,y)\geq\frac{1}{k}$. The former is true since $c\leq\frac{1}{k}$, and the latter is true since $y$ was chosen to be the heaviest coordinate. Additionally, by calculating the probability of any outcome, we can see that Fork is equivalent to the regular sampling process. Define the (random) set $\bm{X}=\{i\,|\,X_{i}=1\}$. We will use $\bm{\theta}$ to refer to a particular realization of this set. We define Fork for sampling from $\hat{\rho}(i,\cdot)$ in the same way, though we will denote the indicator random variables by $\hat{X}_{i}$ and $\bm{\hat{X}}$ instead. Note that, if $c\leq\frac{1}{k}$, the process will still be well defined after rounding. This is because $\hat{\rho}(i,x)\leq c\leq\frac{1}{k}$, and $\hat{\rho}(i,x)+\hat{\rho}(i,y)=\rho(i,x)+\rho(i,y)\geq\frac{1}{k}$. For the rest of this section, when we are drawing a sample from a CRV, we draw it via the process Fork.

The proof of Lemma 1 follows from the following three lemmata. Intuitively, the first states that the PMD induced by the CRVs for which $X_{i}=1$ gives a Poisson Binomial distribution with mean concentrated around its expected value, for both the rounded and unrounded PMDs. The second states that if this value is concentrated, then the two distributions are close in total variation distance. The proof relates the rounded and unrounded distributions by comparing the total variation distance between the Poisson distributions with the same means. The third lemma eliminates the condition on the second lemma by using the first lemma, which states that this condition is likely to hold.

If $\sum_{i\in I^{x}_{y}}\rho(i,x)\geq 3ck\log\left(\frac{1}{ck}\right)$, then

Suppose that, for some $\bm{\theta}$, the following hold:

Then, letting $Z_{i}$ be the Bernoulli random variable with expectation $k\rho_{I^{x}_{y}}(i,x)$ (and $\hat{Z}_{i}$ defined similarly with $k\hat{\rho}_{I^{x}_{y}}(i,x)$),

Since our final rounded $(n,k)$-PMD is generated after applying this rounding procedure $O(k^{2})$ times, Lemma 1 follows from our construction and Lemma 16 via the triangle inequality.

Proof of Lemma 14: Note that $\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,x)=\sum_{i\in I^{x}_{y}}\Omega_{i}$, where

We apply Lemma 11 to the rescaled random variables $\Omega_{i}^{\prime}=\frac{1}{ck}\Omega_{i}$, with $\gamma=\sqrt{\frac{3\log{\frac{1}{ck}}}{E[\sum_{i\in I^{x}_{y}}\Omega_{i}^{\prime}]}}$, giving

Applying the same argument to $\hat{\rho}_{I^{x}_{y}}$ gives

Since $\bm{X}\sim\bm{\hat{X}}$, by considering the joint probability space where $\bm{\theta}=\bm{X}=\bm{\hat{X}}$ and applying a union bound, we get

Proof of Lemma 15: Fix some $\bm{\theta}=\bm{X}=\bm{\hat{X}}$. Without loss of generality, assume E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}\geq E\Big{[}\sum_{i\in\bm{\hat{X}}}k\hat{\rho}_{I^{x}_{y}}(i,j)\Big{]}. There are two cases:

E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}\leq(ck)^{3/4}

From the first assumption in the lemma statement,

Similarly, by the second assumption in the lemma statement and since E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}\geq E\Big{[}\sum_{i\in\bm{\hat{X}}}k\hat{\rho}_{I^{x}_{y}}(i,j)\Big{]}, we also have that $\sum_{i\in\bm{\theta}}k\hat{\rho}_{I^{x}_{y}}(i,j)\leq g(c,k)$.

By Markov’s inequality, \Pr\Big{[}\sum_{i\in\bm{\theta}}Z_{i}\geq 1\Big{]}\leq\sum_{i\in\bm{\theta}}k\rho_{I^{x}_{y}}(i,j)\leq g(c,k), and similarly, \Pr\Big{[}\sum_{i\in\bm{\theta}}\hat{Z}_{i}\geq 1\Big{]}\leq g(c,k). This implies that

E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}\geq(ck)^{3/4}

We use the following proposition, which is a combination of a classical result in Poisson approximation [BHJ92] and Lemma 3.10 in [DP07].

For any set of independent Bernoulli random variables $\{Z_{i}\}_{i}$ with expectations $E[Z_{i}]\leq ck$,

We must now bound the distance between the two Poisson distributions. We use the following lemma from [DP08]:

If $\lambda=\lambda_{0}+D$ for some $D>0,\lambda_{0}>0$,

To bound this, we need the following proposition, which we prove below:

Thus, using the triangle inequality and this proposition, for sufficiently small $c$, we get

By comparing Cases 1 and 2, we see that the desired bound holds in both cases.

Proof of Proposition 10: By the definition of our rounding procedure, we observe that

By the assumptions of Lemma 15 and the assumption that E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}\geq E\Big{[}\sum_{i\in\bm{\hat{X}}}k\hat{\rho}_{I^{x}_{y}}(i,j)\Big{]},

From the assumption that E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,x)\Big{]}\geq(ck)^{3/4}, for sufficiently small $c$,

Combining this with the first assumption of Lemma 15,

Similarly, since E\Big{[}\sum_{i\in\bm{\hat{X}}}k\hat{\rho}_{I^{x}_{y}}(i,j)\Big{]}\geq E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,j)\Big{]}-c\geq(ck)^{3/4}-c, for $c$ sufficiently small,

where the last equality follows for $c$ sufficiently small because E\Big{[}\sum_{i\in\bm{X}}k\rho_{I^{x}_{y}}(i,x)\Big{]}\geq(ck)^{3/4}.

From (1) and (2), for $c$ sufficiently small,

from which the proposition statement follows. $\hfill\qed$

Throughout this proof, we will couple the two sampling processes such that $\bm{\theta}\coloneqq\bm{X}=\bm{\hat{X}}$, which is possible since $\bm{X}\sim\bm{\hat{X}}$. Let $\phi$ be the random event that $\bm{\theta}$ satisfies the following conditions:

Suppose that $\phi$ occurs, and fix a $\bm{\theta}$ in this probability space. We start by showing that for such a $\bm{\theta}$,

Let $M^{\rho_{I^{x}_{y}}^{\bm{\theta}}}$ and $M^{\rho_{I^{x}_{y}}^{\bm{\bar{\theta}}}}$ be the $(n,k)$-PMDs induced by the $k$-CRVs in $M^{\rho_{I^{x}_{y}}}$ with indices in $\bm{\theta}$ and not in $\bm{\theta}$, respectively. Define $M^{\hat{\rho}_{I^{x}_{y}}^{\bm{\theta}}}$ and $M^{\hat{\rho}_{I^{x}_{y}}^{\bm{\bar{\theta}}}}$ similarly. We can see

The first inequality is the triangle inequality, the second inequality is because the distributions for $k$-CRVs in $\bm{\bar{\theta}}$ are identical (since we do not change them in our rounding), and the third inequality is Lemma 15.

By the law of total probability for total variation distance,

where the inequality is obtained by applying Lemma 14 and the bound shown above pointwise for $\bm{\theta}$ which satisfy $\phi$.

B.2 Converting to a Discretized Gaussian using the Valiant-Valiant CLT

We will now apply a result by Valiant and Valiant [VV10]. We recall the aforementioned CLT by Valiant and Valiant, Theorem 6, which we restate for convenience.

As we can see from this inequality, there are two issues that may arise and lead to a bad approximation:

$G^{\rho}$ has small variance in some direction (cf. Proposition 6)

$G^{\rho}$ has a large size parameter $n$

We must avoid both of these issues simultaneously – we will apply this result to several carefully chosen sets, and then merge the resulting Gaussians into one using Lemma 2.

The first step is to partition our CRVs into several sets, and then convert the PMDs induced by each set into GMDs (with an appropriately chosen pivot). The original PMD can be sampled by sampling each of these GMDs and then adding their results. In other words, the probability mass function of the PMD is the convolution of the probability mass functions of these GMDs.

We start by partitioning the $k$-CRVs into $k$ sets $S_{1},\dots,S_{k}$, where $S_{j^{\prime}}=\{i\,|\,j^{\prime}=\arg\max_{j}\hat{\pi}(i,j)\}$ and ties are broken by lexicographic ordering. This defines $S_{j^{\prime}}$ to be the set of indices of $k$-CRVs in which $j^{\prime}$ is the heaviest coordinate. Let $M^{\hat{\pi}_{j^{\prime}}}$ be the $(|S_{j^{\prime}}|,k)$-PMD induced by taking the $k$-CRVs in $S_{j^{\prime}}$. For the remainder of this section, we will focus on $S_{k}$, the other cases follow symmetrically.

We convert each CRV in $S_{k}$ into a truncated $k$-CRV by omitting the $k$th coordinate, giving us a $(|S_{k}|,k)$-GMD $G^{\hat{\rho}_{k}}$. Since the $k$th coordinate was the heaviest, we can make the following observation:

$\hat{\rho}_{k}(i,0)\geq\frac{1}{k}$ for all $i\in S_{k}$.

If we tried to apply Theorem 6 to $G^{\hat{\rho}_{k}}$, we would obtain a vacuous result. For instance, if there exists a $j$ such that $\hat{\rho}_{k}(i,j)=0$ for all $i$, the variance in this direction would be and Theorem 6 would give us a trivial result. Therefore, we further partition $S_{k}$ into $2^{k-1}$ sets indexed by $2^{[k-1]}$, where each set contains the elements of $S_{k}$ which are non-zero on its indexing set and zero otherwise. More formally, $S_{k}^{\mathcal{I}}=\{i\,|\,(i\in S_{k})\wedge(\hat{\rho}(i,j)\geq c\ \forall j\in\mathcal{I})\wedge(\hat{\rho}(i,j)=0\ \forall j\not\in\mathcal{I})\}$. For each of these sets, due to our rounding procedure, we know that the variance is non-negligible in each of the non-zero directions. Naively, we would apply the CLT separately to each of these sets. The issue is that merging the resulting $2^{k}$ Gaussians would be costly. Roughly, merging two Gaussians into one incurs a cost proportional to the inverse of the minimum standard deviation of either Gaussian. In order to avoid the cost of merging exponentially many Gaussians with similar variances, before applying the CLT, we group sets $S_{k}^{\mathcal{I}}$ of similar variance together. The resulting collection of Gaussians have variances which increase rapidly, and by merging them in the correct order, we can minimize this cost.