Testing Conditional Independence of Discrete Distributions

arXiv:1711.11560v1 [] 30 Nov 2017

Clément L. Canonne∗
Stanford University
ccanonne@stanford.edu
Daniel M. Kane‡
University of California, San Diego
dakane@cs.ucsd.edu

Ilias Diakonikolas†
University of Southern California
diakonik@usc.edu
Alistair Stewart
University of Southern California
alistais@usc.edu

December 1, 2017

Abstract
We study the problem of testing conditional independence for discrete distributions. Specifically,
given samples from a discrete random variable (X, Y, Z) on domain [`1 ] × [`2 ] × [n], we want to distinguish, with probability at least 2/3, between the case that X and Y are conditionally independent given
Z from the case that (X, Y, Z) is ε-far, in `1 -distance, from every distribution that has this property.
Conditional independence is a concept of central importance in probability and statistics with a range of
applications in various scientific domains. As such, the statistical task of testing conditional independence has been extensively studied in various forms within the statistics and econometrics communities
for nearly a century. Perhaps surprisingly, this problem has not been previously considered in the framework of distribution property testing and in particular no tester with sublinear sample complexity is
known, even for the important special case that the domains of X and Y are binary.
The main algorithmic result of this work is the first conditional independence tester with sublinear
sample complexity for discrete distributions over [`1 ] × [`2 ] × [n]. To complement our upper bounds,
we prove information-theoretic lower bounds establishing that the sample complexity of our algorithm
is optimal, up to constant factors, for a number of settings. Specifically, for the prototypical setting when
`1 , `2 = O(1), we show that the sample complexity of testing conditional independence (upper bound
and matching lower bound) is




Θ max n1/2 /ε2 , min n7/8 /ε, n6/7 /ε8/7
.
To obtain our tester, we employ a variety of tools, including (1) a suitable weighted adaptation of
the flattening technique [DK16], and (2) the design and analysis of an optimal (unbiased) estimator
for the following statistical problem of independent interest: Given a degree-d polynomial Q : Rn →
R and sample access to a distribution p over [n], estimate Q(p1 , . . . , pn ) up to small additive error.
Obtaining tight variance analyses for specific estimators of this form has been a major technical hurdle
in distribution testing (see, e.g., [CDVV14]). As an important contribution of this work, we develop a
general theory providing tight variance bounds for all such estimators. Our lower bounds, established
using the mutual information method, rely on novel constructions of hard instances that may be useful
in other settings.
∗
Supported by a Motwani Postdoctoral Fellowship. Some of this work was performed while visiting USC, and a graduate
student at Columbia University.
†
Supported by NSF Award CCF-1652862 (CAREER) and a Sloan Research Fellowship.
‡
Supported by NSF Award CCF-1553288 (CAREER) and a Sloan Research Fellowship.

1
1.1

Introduction
Background

Suppose we are performing a medical experiment. Our goal is to compare a binary response (Y ) for two
treatments (X), using data obtained at n levels of a possibly confounding factor (Z). We have a collection
of observations grouped in strata (fixed values of Z). The stratified data are summarized in a series of 2 × 2
contingency tables, one for each strata. One of the most important hypotheses in this context is conditional
independence of X and Y given Z. How many observations (X, Y, Z) do we need so that we can confidently
test this hypothesis?
The above scenario is a special case of the following statistical problem: Given samples from a joint
discrete distribution (X, Y, Z), are the random variables X, Y independent conditioned on Z? This is the
problem of testing conditional independence — a fundamental statistical task with applications in a variety
of fields, including medicine, economics and finance, etc. (see, e.g., [MH59, SGS00, WH17] and references
therein). Formally, we have the following definition:
Definition 1.1 (Conditional Independence). Let X, Y, Z be random variables over discrete domains X , Y, Z
respectively. We say that X and Y are conditionally independent given Z, denoted by (X ⊥ Y ) | Z, if for all
(i, j, z) ∈ X ×Y ×Z we have that: Pr[ X = i, Y = j | Z = z ] = Pr[ X = i | Z = z ]·Pr[ Y = j | Z = z ].
Conditional independence is an important concept in probability theory and statistics, and is a widely
used assumption in various scientific disciplines [Daw79]. Specifically, it is a central notion in modeling
causal relations [SGS00] and of crucial importance in graphical modeling [Pea88]. Conditional independence is, in several settings, a direct implication of economic theory. A prototypical such example is the
Markov property of a time series process. The Markov property is a natural property in time series analysis and is broadly used in economics and finance [EO87]. Other examples include distributional Granger
non-causality [Gra80] — which is a particular case of conditional independence — and exogeneity [BH07].
Given the widespread applications of the conditional independence assumption, the statistical task of
testing conditional independence has been studied extensively for nearly a century. In 1924, R. A. Fisher [Fis24]
proposed the notion of partial correlation coefficient, which leads to Fisher’s classical z-test for the case
that the data comes from a multivariate Gaussian distribution. For discrete distributions, conditional independence testing is one of the most common inference questions that arise in the context of contingency
tables [Agr92]. In the field of graphical models, conditional independence testing is a cornerstone in the context of structure learning and testing of Bayesian networks (see, e.g., [Nea03, TBA06, NUU17, CDKS17]
and references therein). Finally, conditional independence testing is a useful tool in recent applications of
machine learning involving fairness [HPS16].
One of the classical conditional independence tests in the discrete setting is the Cochran–Mantel–
Haenszel test [Coc54, MH59], which requires certain strong assumptions about the marginal distributions.
When such assumptions do not hold, a common tester used is a linear combination of χ2 -squared testers
(see, e.g., [Agr92]). However, even for the most basic case of distributions over {0, 1}2 × [n], no finite
sample analysis is known. (Interestingly enough, we note that our tester can be viewed as an appropriately
weighted linear combination of χ2 -squared tests.) A recent line of work in econometrics has been focusing on conditional independence testing in continuous settings [LG96, DM01, SW07, SW08, Son09, GS10,
Hua10, SW14, ZPJS11, BT14, dMASdBP14, WH17]. The theoretical results in these works are asymptotic
in nature, while the finite sample performance of their proposed testers is evaluated via simulations.
In this paper, we will study the property of conditional independence in the framework of distribution
testing. The field of distribution property testing [BFR+ 00] has seen substantial progress in the past decade,
1

see [Rub12, Can15, Gol17] for two recent surveys and books. A large body of the literature has focused on
characterizing the sample size needed to test properties of arbitrary distributions of a given support size. This
regime is fairly well understood: for many properties of interest there exist sample-efficient testers [Pan08,
CDVV14, VV14, DKN15b, ADK15, CDGR16, DK16, DGPP16, CDS17, Gol17, DGPP17]. Moreover,
an emerging body of work has focused on leveraging a priori structure of the underlying distributions
to obtain significantly improved sample complexities [BKR04, DDS+ 13, DKN15b, DKN15a, CDKS17,
DP17, DDK18, DKN17].

1.2

Our Contributions

Rather surprisingly, the problem of testing conditional independence has not been previously considered in
the context of distribution property testing. In this work, we study this problem for discrete distributions
and provide the first conditional independence tester with sublinear sample complexity. To complement our
upper bound, we also provide information-theoretic lower bounds establishing that the sample complexity
of our algorithm is optimal for a number of important regimes. To design and analyze our conditional
independence tester we employ a variety of tools, including an optimal (unbiased) estimator for the following
statistical task of independent interest: Given a degree-d polynomial Q : Rn → R and sample access to a
distribution p over [n], estimate Q(p1 , . . . , pn ) up to small additive error.
In this section, we provide an overview of our results. We start with some terminology. We denote by
∆(Ω) the set of all distributions over domain Ω. For discrete sets X , Y, Z, we will use PX ,Y|Z to denote
the property of conditional independence, i.e.,
PX ,Y|Z := { p ∈ ∆(X × Y × Z) : (X, Y, Z) ∼ p satisfies (X ⊥ Y ) | Z } .
We say that a distribution p ∈ ∆(X × Y × Z) is ε-far from PX ,Y|Z , if for every distribution q ∈ PX ,Y|Z
we have that dTV (p, q) > ε. We study the following hypothesis testing problem:
T (`1 , `2 , n, ε): Given sample access to a distribution p over X × Y × Z, with |X | = `1 , |Y| = `2 ,
|Z| = n, and ε > 0, distinguish with probability at least 2/3 between the following cases:
• Completeness: p ∈ PX ,Y|Z .
• Soundness: dTV (p, PX ,Y|Z ) ≥ ε.
Even though the focus of this paper is on testing under the total variation distance metric (or equivalently
the `1 -distance), we remark that our techniques yield near-optimal algorithms under the mutual information
metric as well. The interested reader is referred to Appendix A for a short description of these implications.
The property of conditional independence captures a number of other important properties as a special
case. For example, the n = 1 case reduces to the property of independence over [`1 ] × [`2 ], whose testing
sample complexity was resolved only recently [DK16]. Arguably the prototypical regime of conditional
independence corresponds to the other extreme. That is, the setting that the domains X , Y are binary (or,
more generally, of small constant size), while the domain Z is large. This regime exactly captures the wellstudied and practically relevant setting of 2×2×n contingency tables (mentioned in the motivating example
of the previous section). For the setting where X , Y are small, our tester and our sample complexity lower
bound match, up to constant factors. Specifically, we prove:
Theorem 1.1. There exists a computationally efficient tester for T (2, 2, n, ε) with sample complexity






O max n1/2 /ε2 , min n7/8 /ε, n6/7 /ε8/7
2



.

Moreover, this sample
upper
bound is tight, up to constant


factors. That is, any tester for T (2, 2, n, ε)
1/2
2
7/8
6/7
8/7
requires at least Ω max n /ε , min n /ε, n /ε
samples.
To the best of our knowledge, prior to our work, no o(n) sample algorithm was known for this problem.
Our algorithm in this regime is simple: For every fixed value of z ∈ [n], we consider the conditional distribution pz . Note that pz is a distribution over X × Y. We construct an unbiased estimator Φ of the squared
`2 -distance of any distribution on X × Y from the product of its marginals. Our conditional independence
tester uses this estimator in a black-box manner for each of the pz ’s. In more detail, our tester computes
a weighted linear combination of Φ(pz ), z ∈ [n], and rejects if and only if this exceeds an appropriate
threshold.
To obtain the required unbiased estimator of the squared `2 -distance, we observe that this task is a
special case of the following more general problem of broader interest: For a distribution p = (p1 , . . . , pn )
and a polynomial Q : Rn → R, obtain an unbiased estimator for the quantity Q(p1 , . . . , pn ). We prove the
following general result:
Theorem 1.2. For any degree-d polynomial Q : Rn → R and distribution p over [n], there exists a unique
and explicit unbiased estimator UN for Q(p) given N ≥ d samples. Moreover, this estimator is linear in Q
and its variance is at most
X

n
Y
p si
i

s∈Nn
1≤ksk≤d

i=1

!

∂ ksk Q(p)
∂X1s1 . . . ∂Xnsn

!2 

(N − ksk)!
Q
N ! ni=1 si !



,

which itself can be further bounded as a function of Q+ , the degree-d polynomial obtained by taking the
absolute values of all the coefficients of Q, and its partial derivatives.
We note that Theorem 1.2 can be appropriately extended to the setting where we are interested in estimating Q(p, q), where p, q are discrete distributions over [n] and Q is a real degree-d polynomial on 2n
variables. In addition to being a crucial ingredient for our general conditional independence tester, we believe that Theorem 1.2 is of independent interest. Indeed, in a number of distribution testing problems,
we need unbiased estimators for some specific polynomial Q of a distribution p (or a pair of distributions
p, q). For example, the `2 -tester of [CDVV14] (which has been used as a primitive to obtain a wide range of
sample-optimal testers [DK16]) is an unbiased estimator for the squared `2 -distance between two distributions p, q over [n]. While the description of such unbiased estimators may be relatively simple, their analyses
are typically highly non-trivial. Specifically, obtaining tight bounds for the variance of such estimators has
been a major technical hurdle in distribution testing. As an important contribution of this work, we develop
a general theory providing tight variance bounds for all such estimators.
The conditional independence tester Theorem 1.1 straightforwardly extends to larger domains X , Y,
alas its sample complexity becomes at least linear in the size of these sets. To obtain a sublinear tester for
this general case, we require a number of additional conceptual and technical ideas. Our main theorem for
conditional independence testing for domain [`1 ] × [`2 ] × [n] is the following:
Theorem 1.3. There exists a computationally efficient tester for T (`1 , `2 , n, ε) with sample complexity
1/4 1/4

2/7 2/7

n7/8 `1 `2 n6/7 `1 `2
O max min
,
ε
ε8/7

1/2 1/2

!

2/3 1/3

1/2 1/2

n3/4 `1 `2 n2/3 `1 `2 n1/2 `1 `2
,
,
,
ε
ε2
ε4/3

where we assume without loss of generality that `1 ≥ `2 .
3

!!

,

(1)

The expression of the sample complexity in Theorem 1.3 may seem somewhat unwieldy. In an attempt
to interpret this bound, we consider several important special cases of interest:
• For `1 = `2 = O(1), (1) reduces to the binary case for X, Y , recovering the tight bound of Theorem 1.1.
• For n = 1, our problem reduces to the task of testing independence of a distribution over [`1 ] × [`2 ],
which has been extensively studied [BFF+ 01, LRR11, ADK15,
(1) recovers
the


DK16]. In this case,
2/3 1/3 4/3 √
2
[DK16].
optimal sample complexity of independence testing, i.e., Θ max `1 `2 /ε , `1 `2 /ε
• For `1 = `2 = n (and ε = Ω(1)), the sample complexity of (1) becomes O(n7/4 ). In Theorem 1.4
below, we show that this bound is optimal as well.
We conclude with the aforementioned tight sample lower bound for constant values of ε, in the setting
where all three coordinates are of approximately the same cardinality:
Theorem 1.4. Any tester for T (n, n, n, 1/20) requires Ω(n7/4 ) samples.

1.3
1.3.1

Our Techniques
Conditional Independence Tester for Binary X , Y

In the case where X and Y are binary, for each bin z ∈ Z we will attempt to estimate the squared `2 -distance
of the corresponding conditional distribution and the product of its conditional marginals. In particular, if
X = Y = {0, 1} the square of p00 p11 − p01 p10 , where pij is the probability that X = i and Y = j, for
Z = z, is proportional to this difference. Since this square is a degree-4 polynomial in the samples, there
is an unbiased estimator of this quantity that can be computed for any value z ∈ Z from which we have at
least 4 samples. Furthermore, for values of z ∈ Z for which we have more than 4 samples, the additional
samples can be used to reduce the error of this estimator. The final algorithm computes a weighted linear
combination of these estimators (weighted so that the more accurate estimators from heavier bins are given
more weight) and compares it to an appropriate threshold. The correctness of this estimator requires a rather
subtle analysis. Recall that there are three different regimes of ε versus n in the optimal sample complexity
and the tester achieves this bound without a case analysis. As usual, we require a bound on the variance of
our estimator and a lower bound on the expectation in the soundness case.
On the one hand, a naive bound on the variance for our estimator for an individual bin turns out to be
insufficient for our analysis. In particular, let p be a discrete probability distribution and Q(p) a polynomial
in the individual probabilities of p. Given m ≥ deg(Q) independent samples from p, it is easy to see that
there is a unique symmetric, unbiased estimator for Q(p), which we call Um Q. Our analysis will rely on
obtaining tight bounds for the variance of Um Q. It is not hard to show that this variance scales as O(1/m),
but this bound turns out to be insufficient for our purposes. In order to refine this estimate, we show that
Var(Um Q) = R(p)/m + O(1/m2 ), for some polynomial R for which we devise a general formula. From
this point on, we can show that for our polynomial Q (or in general any Q which is the square of a lower
degree polynomial) Var(Um Q) = O(Q(p)/m + 1/m2 ). This gives a much sharper estimate on the variance
of our estimator, except in cases where the mean is large enough that the extra precision is not necessary.
Another technical piece of our analysis is relating the mean of our estimator to the total variation distance
of our distribution from being conditionally independent. In particular, our estimator is roughly the sum
(over the Z-bins with enough samples) of the squared `2 distance that the conditional distribution is from
being independent. When much of the distance from conditionally independence comes from relatively
heavy bins, this relation is a more or less standard `1 /`2 inequality. However, when the discrepancy is

4

concentrated on very light bins, the effectiveness of our tester is bounded by the number of these bins which
obtain at least four samples, and a somewhat different analysis is required. In fact, out of the different cases
in the performance of our algorithm, one of the boundaries is determined by a transition between the hard
cases involving discrepancies supported on light bins to ones where the discrepancy is supported on heavy
bins.
If the variables X and Y are no longer binary, our estimates for the discrepancy of an individual bin must
be updated. In particular, we similarly use an unbiased estimator of the `2 distance between the conditional
distribution and the product of its conditional marginals. We note however that variance of this estimator is
large if the marginal distributions have large `2 norms. Therefore, in bins for which we have a large number
of samples, we can employ an idea from [DK16] and use some of our samples to artificially break up the
heavier bins, thus flattening these distributions. We elaborate on this case, and the required ingredients it
entails, in the next subsection.
1.3.2

General Conditional Independence Tester

Assuming that we take at least four samples from any bin z ∈ Z, we can compute an unbiased estimator
for the squared `2 distance between pz , the conditional distribution, and qz the product of its conditional
marginals. It is easy to see that this expectation is at least ε2z /(|X | |Y|), where εz is the `1 distance between
the conditional distribution and the closest distribution with independent X and Y coordinates. At a high
level, our algorithm takes a linear combination of these bin-wise estimators (over all bins from which we
got at least 4 samples), and compares it to an appropriate threshold. There is a number of key ideas that are
needed so that this approach gives us the right sample complexity.
Firstly, we use the idea of flattening, introduced in [DK16]. The idea here is that the variance of the `2
estimator is larger if the `2 norms of p and q are large. However, we can reduce this variance by artificially
breaking up the heavy bins. In particular, if we have m samples from a discrete distribution of support
size n, we can artificially add m bins and reduce the `2 norm of the distribution (in expectation) to at most
√
O(1/ m). We note that it is usually not a good idea to employ this operation for m  n, as it will
substantially increase the number of bins. Nor do we want to use all of our samples for flattening (since we
need to use some for the actual tester). Trading off these considerations, using min(m/2, n) of our samples
to flatten is a reasonable choice. We also remark that instead of thinking of p and q as distributions over
|X | |Y| bins, we exploit the fact that q is a two-dimensional product distribution over |X |×|Y|. By flattening
these marginal distributions independently, we can obtain substantially better variance upper bounds.
Secondly, we need to use appropriate weights for our bin-wise estimator. To begin with, one might
imagine that the weight we should use for the estimator of a bin z ∈ Z should be proportional to the
probability mass of that bin. This is a natural choice because heavier bins will contribute more to the final
`1 error, and thus, we will want to consider their effects more strongly. The probability mass of a bin is
approximately proportional to the number of samples obtained from that bin. Therefore, we might want to
weight each bin by the number of samples drawn from it. However, there is another important effect of
having more samples in a given bin. In particular, having more samples from a bin allows us to do more
flattening of that bin, which decreases the variance of the corresponding bin-wise estimator. This means that
we will want to assign more weight to these bins based on how much flattening is being done, as they will
give us more accurate information about the behavior of that bin.
Finally, we need to analyze our algorithm. Let m be the number of samples we take and n be the
domain of Z. If the
p bin weights are chosen appropriately, we show that the final estimator A has variance
O(min(n, m) + min(n, m)E[A] + E[A]3/2 ), with high probability over the number of samples falling
in each bin as well as the flattening we perform for each bin. (This high-probability statement, in turn,
5

is enough for us to apply Chebyshev’s inequality to our final estimator in the end.) Furthermore, in the
completeness case, we have that E[A] = 0. In order to be able to distinguish between completeness and
soundness, we
pneed it to be the case that for all distributions ε-far from conditional independence it holds
that E[A]  min(n, m). We know that if we are ε-far from conditional independence, we must have that
P
z εz wz  ε, where wz is the probability that Z = z. In order to take advantage of this fact, we will need
to separate the Z-bins into four categories based on the size of the wz . Indeed, if we are far from conditional
independence, then for at least one of these cases the sum of εz wz over bins of that type only will be  ε.
Each of these four cases will require a slightly different analysis:
• Case 1: wz < 1/m. In this case, the expected number of samples from bin z is small. In particular,
the probability of even seeing 4 samples from the bin might will be small. Here, the expectation is
dominated by the probability that we see enough samples from the bin.
• Case 2: 1/m < wz < |X | /m: In this case, we are likely to get our 4 samples from the bin, but
probably will get fewer than |X |. This means that our flattening will not saturate either of the marginal
distributions and we can reduce the squared `2 norm of q by a full factor of mz (where mz is the
number of samples from this bin).
• Case 3: |X | /m < wz < |Y| /m. In this case, we are likely to saturate our flattening over the Xmarginal but not the Y -marginal. Thus,
p our flattening only decreases the `2 norm of the conditional
distribution on that bin by a factor of |X | mz .
• Case 4: |Y| /m < wz : Finally, in thispcase we saturate both the X- and Y -marginals, so our flattening
decreases the `2 norm by a factor of |X | |Y|.
Within each sub-case, the expectation of A is a polynomial in m, |X | , |Y| multiplied by the sum over z ∈ Z
P
of some polynomial in εz and wz . We need to bound this
from
below
given
that
z εz wz  ε, and then set
p
m large enough so that this lower bound is more than min(n, m). We note that only in Case 1 is the case
where m < n relevant. Thus, our final bound will be a maximum over the 4 cases of the m required in the
appropriate case.
1.3.3

Sample Complexity Lower Bound Construction for Binary X , Y

We begin by reviewing the lower bound methodology we follow: In this methodology, a lower bound is
shown by adversarially constructing two distributions over pseudo-distributions. Specifically, we construct a
pair of ensembles D and D0 of pairs of nearly-normalized pseudo-distributions such that the following holds:
(1) Pseudo-distributions drawn from D satisfy the desired property and pseudo-distribution drawn from D0
are ε-far from satisfying the property with high probability, and (2) Poisson(s) samples are insufficient to
reliably determine from which ensemble the distribution was taken from, unless s is large enough.
To formally prove our lower bounds, we will use the mutual information method, as in [DK16]. In this
section, we provide an intuitive description of our sample complexity lower bound for testing conditional
independence, when X = Y = {0, 1} and Z = [n]. (Our lower bound for the regime X = Y = Z = [n] is
proved using the same methodology, but relies on a different construction.) We construct ensembles D and
D0 — where draws from D are conditionally independent and draws from D0 are ε-far from conditionally
independent with high probability — and show that s samples from a distribution on (X, Y, Z) are insufficient to reliably distinguish whether the distribution came from D or D0 , when s is small. We define D
and D0 by treating each bin z ∈ [n] of Z independently. In particular, for each possible value z ∈ [n] for
Z, we proceed as follows: (1) With probability min(s/n, 1/2), we assign the point Z = z probability mass
max(1/s, 1/n) and let the conditional distribution on (X, Y ) be uniform. Since the distribution is conditionally independent on these bins and identical in both ensembles, these “heavy” bins will create “noise” to
6

confuse an estimator. (2) With probability 1 − min(s/n, 1/2), we set the probability that Z = z to be ε/n,
and let the conditional distribution on (X, Y ) be taken from either C or C 0 , for some specific ensembles C
and C 0 . In particular, we pick C and C 0 so that a draw from C is independent and a draw from C 0 is far
from independent. These bins provide the useful information that allows us to distinguish between the two
ensembles D and D0 . The crucial property is that we can achieve the above while guaranteeing that any
third moment from C agrees with the corresponding third moment from C 0 . This guarantee implies that if
we draw 3 (or fewer) samples of (X, Y ) from some bin Z = z, then the distribution on triples of (X, Y )
will be identical if the conditional was taken from C or if it was taken from C 0 . That is, all information
about whether our distribution came from D or D0 will come from bins of type (2) (for which we have at
least 4 samples) of which there will be approximately n(sε/n)4 . On the other hand, there will be about
min(s, n) bins of type (1) with 4 samples in random
configuration adding “noise”. Thus, we will not be
p
able to distinguish reliably unless n(sε/n)4  min(s, n), as otherwise the noise’ due to the heavy bins
will drown out the signal of the light ones.
To define C and C 0 , we find appropriate vectors p, q over {0, 1}2 so that p + q and p + 3q each are
distributions with independent coordinates, but p, p+2q, p+4q are not. We let C return p+q and p+3q each
with probability 1/2, and let C 0 return p, p + 2q or p + 4q with probability 1/8, 3/4, 1/8 respectively. If we
wish to find the probability that 3 samples from a distribution r come in some particular pattern, we get f (r)
for some degree-3 polynomial f . If we want the difference in these probabilities for r a random draw from
C and a random draw from C 0 , we get f (p + q)/2 + f (p + 3q)/2 − f (p)/8 − f (p + 2q)(3/4) − f (p + 4q)/8.
We note that this is proportional to the fourth finite difference of a degree-3 polynomial, and is thus 0.
Therefore, any combination of at most 3 samples are equally likely to show up for some Z-bin from D as
from D0 .
To rigorously analyze the above sketched construction, we consider drawing Poisson(s) samples from
a random distribution from either D or D0 , and bound the mutual information between the set of samples
and the ensemble they came from. Since the samples from each bin are conditionally independent on the
ensemble, this is at most n times the mutual information coming from a single bin. By the above, the
probabilities of seeing any triple of samples are the same for either D or D0 and thus contribute nothing to
the mutual information. For sets of 4 or more samples, we note that the difference in probabilities comes
only from the case where 4 samples are drawn from a bin of type (2), which happens with probability at
most O(sε/n)4 . However, this is counterbalanced by the fact that these sample patterns are seen with much
higher frequency from bins of type (1) (as they have larger overall mass). Thus, the mutual information for
a combination including m ≥ 4 samples will be (O(sε/n)m )2 / min(s/n, 1/2) · Ω(1)m . The contribution
from m > 4 can be shown to be negligible, thus the total mutual information summed over all bins is
O(min(s, n) · (sε/n)8 ). This must be Ω(1) in order to reliably distinguish, and this proves our lower bound.

1.4

Organization

After setting up the required preliminaries in Section 2, we give our testing algorithm for the case of constant
|X | , |Y| in Section 3. In Section 4, we develop our theory for polynomial estimation. Section 5 leverages this
theory, along with several other ideas, to obtain our general testing algorithm for conditional independence.
Section 6 gives our information-theoretic lower bound for the setting of binary |X | , |Y|. In Section 7,
we give an information-theoretic lower bound matching the sample complexity of our algorithm for the
regime where |X | = |Y| = |Z|. In Appendix A, we discuss the implications of our results for conditional
independence testing with regard to the conditional mutual information.

7

2

Preliminaries and Basic Facts

We begin with some standard notation and definitions that we shall use throughout the paper. For m ∈ N,
we write [m] for the set {1, . . . , m}, and log for the binary logarithm.
Distributions and Metrics A probability distribution over discrete domain Ω is a function p : Ω → [0, 1]
P
such that kpk1 := ω∈Ω p(ω) = 1. Without the requirement that the total mass be one, p is said to be a
pseudo-distribution. We denote by ∆(Ω) the set of all probability distributions over domain Ω. For two
probability distributions p, q ∈ ∆(Ω), their total variation distance (or statistical distance) is defined as
P
dTV (p, q) := supS⊆Ω (p(S) − q(S)) = 21 ω∈Ω |p(ω) − q(ω)| , i.e., dTV (p, q) = 12 kp − qk1 , and their `2
distance is the distance kp − qk2 between their probability mass functions. Given a subset P ⊆ ∆(Ω) of
distributions, the distance from p to P is then defined as dTV (p, P) := inf q∈P dTV (p, q). If dTV (p, P) > ε,
we say that p is ε-far from P; otherwise, it is ε-close. For a distribution p we write X ∼ p to denote
that the random variable X is distributed according to p. Finally, for p ∈ ∆(Ω1 ) , q ∈ ∆(Ω2 ), we let
p ⊗ q ∈ ∆(Ω1 × Ω2 ) be the product distribution with marginals p and q.
Property Testing We work in the standard setting of distribution testing: a testing algorithm for a property
P ⊆ ∆(Ω) is an algorithm which, granted access to independent samples from an unknown distribution
p ∈ ∆(Ω) as well as distance parameter ε ∈ (0, 1], outputs either accept or reject, with the following
guarantees:
• If p ∈ P, then it outputs accept with probability at least 2/3.
• If dTV (p, P) > ε, then it outputs reject with probability at least 2/3.
The two measures of interest here are the sample complexity of the algorithm (i.e., the number of samples it
draws from the underlying distribution) and its running time.

2.1

Conditional Independence

We record here a number of notations definitions regarding conditional independence. Let X, Y, Z be random variables over discrete domains X , Y, Z respectively. Given samples from the joint distribution of
(X, Y, Z), we want to determine whether X and Y are conditionally independent given Z, denoted by
(X ⊥ Y ) | Z, versus ε-far in total variation distance from every distribution of random variables (X 0 , Y 0 , Z 0 )
such that (X 0 ⊥ Y 0 ) | Z 0 . For discrete sets X , Y, Z, we will denote by PX ,Y|Z the property of conditional
independence, i.e., PX ,Y|Z := { p ∈ ∆(X × Y × Z) : (X, Y, Z) ∼ p satisfies (X ⊥ Y ) | Z }. We say
that a distribution p ∈ ∆(X × Y × Z) is ε-far from PX ,Y|Z , if for every distribution q ∈ PX ,Y|Z we have
that dTV (p, q) > ε. Fix a distribution q ∈ PX ,Y|Z of minimum total variation distance to p. Then the
marginals of q on each of the three coordinates may have different distributions. We will also consider
testing conditional independence with respect to a different metric, namely the conditional mutual information [Dob59, Wyn78]. For three random variables X, Y, Z as above, the conditional mutual information
of X and Y with respect to Z is defined as I(X; Y |Z) := EZ [(I (X; Y ) | Z)], i.e., as the expected (with
respect to Z) K-L divergence between the distributions of (X, Y ) | Z and the product of the distributions of
(X | Z) and (Y | Z). In this variant of the problem (considered in Appendix A), we will want to distinguish
I(X; Y |Z) = 0 from I(X; Y |Z) ≥ ε.
Notation. Let p ∈ ∆(X × Y × Z). For z ∈ Z, we will denote by pz ∈ ∆(X × Y) the distribution
defined by pz (i, j) := Pr(X,Y,Z)∼p [ X = i, Y = j | Z = z ] and by pZ ∈ ∆(Z) the distribution pZ (z) :=
8

Pr(X,Y,Z)∼p [ Z = z ]. By definition, for any p ∈ ∆(X × Y × Z), we have that p(i, j, z) = pZ (z) · pz (i, j).
For z ∈ Z, we will denote by pz,X ∈ ∆(X ) the distribution pz,X (i) = Pr(X,Y,Z)∼p [ X = i | Z = z ] and
pz,Y ∈ ∆(Y) the distribution pz,Y (j) = Pr(X,Y,Z)∼p [ Y = j | Z = z ].
We can now define the product distribution of the conditional marginals:
Definition 2.1 (Product of Conditional Marginals). Let p ∈ ∆(X × Y × Z). For z ∈ Z, we define the
product of conditional marginals of p given Z = z to be the product distribution qz ∈ ∆(X × Y) defined
by qz := pz,X ⊗ pz,Y , i.e., qz (i, j) = pz,X (i) · pz,Y (j). We will also denote by q the mixture of product
P
distributions q := z∈Z pZ (z)qz ∈ PX ,Y|Z , i.e., q(i, j, z) := pZ (z) · qz (i, j).

2.2

Basic Facts

We start with the following simple lemma:
Lemma 2.1. Let p, p0 ∈ ∆(X × Y × Z). Then we have that
dTV p, p0 ≤


X

pZ (z) · dTV pz , p0z + dTV pZ , p0Z ,




(2)

z∈Z

with equality if and only if pZ = p0Z .
Using Lemma 2.1, we deduce the following useful corollary:
Fact 2.1. If p ∈ ∆(X × Y × Z) is ε-far from PX ,Y|Z , then, for every p0 ∈ PX ,Y|Z , either (i) dTV (pZ , p0Z ) >
P
ε/2, or (ii) z∈Z pZ (z) · dTV (pz , p0z ) > ε/2.
Proof. Let p0 ∈ PX ,Y|Z . Since p is ε-far from PX ,Y|Z we have that dTV (p, p0 ) > ε. By Lemma 2.1, we this
P
obtain that z∈Z pZ (z) · dTV (pz , p0z ) + dTV (pZ , p0Z ) > ε, which proves the fact.
The next lemma shows a useful structural property of conditional independence that will be crucial for
our algorithm. It shows that if a distribution p ∈ ∆(X × Y × Z) is close to being conditionally independent, then it is also close to an appropriate mixture of its products of conditional marginals, specifically
distribution q from Definition 2.1:
Lemma 2.2. Suppose p ∈ ∆(X × Y × Z) is ε-close to PX ,Y|Z . Then, p is 4ε-close to the distribution
P
q = z∈Z pZ (z)qz .

2.3

Flattening Distributions

We now recall some notions and lemmata from previous work, regarding the technique of flattening of
discrete distributions:
Definition 2.2 (Split distribution [DK16]). Given a distribution p ∈ ∆([n]) and a multi-set S of elements
of [n], define the split distribution pS ∈ ∆([n + |S|]) as follows: For 1 ≤ i ≤ n, let ai denote 1 plus
P
the number of elements of S that are equal to i. Thus, ni=1 ai = n + |S|. We can therefore associate
the elements of [n + |S|] to elements of the set BS := { (i, j) : i ∈ [n], 1 ≤ j ≤ ai }. We now define a
distribution pS with support BS , by letting a random sample from pS be given by (i, j), where i is drawn
randomly from p and j is drawn uniformly from [ai ].
9

Fact 2.2 ([DK16, Fact 2.5]). Let p, q ∈ ∆([n]), and S a given multi-set of [n]. Then: (i) We can simulate
a sample from pS or qS by taking a single sample from p or q, respectively. (ii) It holds dTV (pS , qS ) =
dTV (p, q).
We will also require the analogue of [DK16, Lemma 2.6] (how flattening reduces the `2 -norm of a distribution) for the non-Poissonized setting, i.e., when exactly m samples are drawn (instead of Poisson(m)).
The proof of this lemma is similar to that of [DK16, Lemma 2.6], and we include it in Appendix B for
completeness:
Lemma 2.3. Let p ∈ ∆([n]). Then: (i) For any multi-setsh S ⊆ Si0 of [n], kpS 0 k2 ≤ kpS k2 , and (ii) If S is
1
obtained by taking m independent samples from p, then E [kpS k22 ≤ m+1
.
Remark 2.3. Given S and (ai )i∈[n] as in Definition 2.2, it is immediate that for any p, q ∈ ∆([n]) it holds
kpS − qS k22 =
for p, q.

2.4

Pn

i=1

(pi −qi )2
ai

so that an `22 statistic for pS , qS can be seen as a particular rescaled `2 statistic

Technical Facts on Poisson Random Variables

We state below some technical result on moments of truncated Poisson random variables, which we will use
in various places of our analysis. The proof of these claims are deferred to Appendix B.
Claim 2.1. There exists an absolute constant C > 0 such that, for N ∼ Poisson(λ),
h

i

Var[N 1{N ≥4} ] ≤ CE N 1{N ≥4} .
Moreover, one can take C = 4.22.
Claim 2.2. There exists an absolute constant C > 0 such that, for X ∼ Poisson(λ) and a, b ≥ 0,
 q

q



Var[X min(X, a) min(X, b)1{X≥4} ] ≤ CE X min(X, a) min(X, b)1{X≥4} .
Claim 2.3. There exists an absolute constant C > 0 such that, for X ∼ Poisson(λ) and integers a, b ≥ 2,
 q



q

E X min(X, a) min(X, b)1{X≥4} ≥ C min(λ min(λ, a) min(λ, b), λ4 ) .

3

Conditional Independence Tester: The Case of Constant |X |, |Y|

Let p ∈ ∆(X × Y × z). In this section, we present and analyze our conditional independence tester for the
case that |X |, |Y| = O(1). Specifically, we will present a tester for this regime whose sample complexity
is optimal, up to constant factors. Our tester uses as a black-box an unbiased estimator for the `22 -distance
between a 2-dimensional distribution and the product of its marginals. Specifically, we assume that we
have access to an estimator Φ with the following performance: Given N samples s = (s1 , . . . , sN ) from a
distribution p ∈ ∆(X × Y), Φ satisfies:
E[Φ(s)] = kp − pX ⊗ pY k22


E[Φ(s)]
1
Var[Φ(s)] ≤ C
+ 2 ,
N
N

(3)
(4)

for some absolute constant C > 0. Such an estimator follows as a special case of our generic polynomial
estimators of Section 4.
10

Notation Let p ∈ ∆(X × Y). We denote its marginal distributions by pX , pY . That is, we have that
pX ∈ ∆(X ) with pX (x) := Pr(X,Y )∼p [ X = x ], x ∈ X , and similarly for pY . Let p ∈ ∆(X × Y × Z).
For z ∈ Z, we will denote by qz the product distribution pz,X ⊗ px,Y .
Let M be a Poisson(m) random variable representing the number of samples drawn from p ∈ ∆(X × Y × Z).
Given the multi-set S of M samples drawn from p, let Sz := { (x, y) : (x, y, z) ∈ S } denote the multi-set
of pairs (x, y) ∈ X × Y corresponding to samples (x, y, z) ∈ S, i.e., the multi-set of samples coming from
the conditional distribution pz . For convenience, we will use the notation σz := |Sz |. Let
Az := σz · Φ(Sz ) · 1{σz ≥4} ,
for all z ∈ Z. Our final statistic (that we will compare to a suitable threshold in the eventual test) is
A :=

X

Az .

z∈Z

We set ε0 := √

ε
|X ||Y|

= Θ(ε), and choose
m ≥ β max

√

2



n/ε0 , min n7/8 /ε0 , n6/7 /ε0

8/7



,

(5)

for a sufficiently large absolute constant β > 0.
Interestingly enough, there are three distinct regions for this expression, based on the relation between
n and ε, as illustrated in the following figure:

Figure 1: The three regimes of the sample complexity (for n = 10000 and ε ∈ (0, 1]) (log-log plot).
Our conditional
tester
 independence

 outputs “accept” if A ≥ τ and “reject” otherwise, where τ is
2 m4 ε0 4
0
selected to be Θ max mε , n3
. A detailed pseudo-code for the algorithm is given in Algorithm 1.
11

Algorithm 1 T EST C OND I NDEPENDENCE
Require: Parameter n := |Z|, `1 :=
(0, 1], and sample access to p ∈ ∆(X × Y × Z).
 |X |, `2 := |Y|, ε ∈
√
√
2
8/7
0
7/8
0
6/7
0
1: Set m ← β max
n/ε , min n /ε , n /ε
, where ε0 := ε/ `1 `2 . β ≥ 1 is a sufficiently
large absolute constant


4 04
2
γ
5
, where γ := 1 − 2e
2: Set τ ← 2 max mε0 , mnε3
.
. Threshold for accepting
3: Draw M ∼ Poisson(m) samples from p and let S be the multi-set of samples.
4: for all z ∈ Z do
5:
Let Sz ⊆ X × Y be the multi-set Sz := { (x, y) : (x, y, z) ∈ S }.
6:
if |Sz | ≥ 4 then
. Enough samples to call Φ
7:
Compute Φ(Sz ).
8:
Set Az ← |Sz | · Φ(Sz ).
9:
else
10:
Set Az ← 0.
11:
end if
12: end for
P
13: if A :=
z∈Z Az ≥ τ then
14:
return accept
15: else
16:
return reject
17: end if

3.1

Proof of Correctness

In this section, we prove correctness of Algorithm 1. Specifically, we will show that: (1) If p ∈ PX ,Y|Z
(completeness), then Algorithm 1 outputs “accept” with probability at least 2/3, and (2) If dTV (p, PX ,Y|Z ) >
ε, then Algorithm 1 outputs “reject” with probability at least 2/3. The proof proceeds by analyzing the
expectation and variance of our statistic A and using Chebyshev’s inequality. We note that β, γ are absolute
constants defined in the algorithm pseudo-code.
3.1.1

Analyzing the Expectation of A

The main result of this subsection is the following proposition establishing the existence of a gap in the
expected value of A in the completeness and soundness cases:


Proposition 3.1. We have the following: (a) If p ∈ PX ,Y|Z , then E[A] = 0. (b) If dTV p, PX ,Y|Z


then E[A] > γ min mε0 2 ,

m4 ε0 4
8n3





> ε,

p
≥ β·γ
min(n, m).
8 ·

The rest of this subsection is devoted to the proof of Proposition 3.1. We start by providing a convenient
lower bound on the expectation of A. We prove the following lemma:
Lemma 3.1. For z ∈ Z, let δz := kpz − qz k2 and αz := m · pZ (z). Then, we have that:
E[A] ≥ γ ·

X

δz2 min(αz , αz4 ) .

(6)

z∈Z

Proof. Conditioned on σz = |Sz |, Eq. (3) gives that E[ Az | σz ] = σz δz2 1{σz ≥4} . Therefore, for σ :=
P
(σz )z∈Z , we can write E[ A | σ ] = z∈Z σz δz2 1{σz ≥4} . Using the fact that the σz ’s are independent and
12

σz ∼ Poisson(αz ), we obtain the following closed-form expression for the expectation:
E[A] = E[E[ A | σ ]]σ =

X

h

i

X

δz2 E σz 1{σz ≥4} =

z∈Z

δz2 · f (αz ) ,

(7)

z∈Z
k

3

x
x
−x
2
where f : R+ → R is the function f (x) = e−x ∞
k=4 k k! = x − e (x + x + 2 ) . Let g : R+ → R be
4
defined by g(x) = min(x, x ) . It is not hard to check that the function f (x)/g(x) achieves its minimum
5
at x = 1, where it takes the value γ := 1 − 2e
> 0. That is, f (αz ) ≥ γ · g(αz ) and the lemma follows
from (7).

P

Given (7), the first statement of Proposition 3.1 is immediate. Indeed, if p is conditionally independent,
then all δz ’s are zero. To establish the second statement, we will require a number of intermediate
lemmata.


We henceforth focus on the analysis of the soundness case, i.e., we will assume that dTV p, PX ,Y|Z > ε.
We require the following useful claim:




Claim 3.1. If dTV p, PX ,Y|Z > ε, then

P

> 2mε0 .

z∈Z δz αz

Proof. We use the identity dTV (p, q) = z∈Z pZ (z) · dTV (pz , qz ), which follows from Lemma 2.1 noting
that qZ = pZ . By assumption, we have that dTV (p, q) > ε. We can therefore write
P

X

δz αz = m

X

kpz − qz k2 · pZ (z) ≥ p

z∈Z

z∈Z

X
2m
m
kpz − qz k1 pZ (z) > p
·
ε,
|X | |Y| z∈Z
|X | |Y|

where the last inequality is Cauchy–Schwarz.
Lemma 3.1 suggests the existence of two distinct regimes: the value of the expectation of our statistic
is dominated by (1) the “heavy” elements z ∈ Z for which αz > 1, or (2) the “light” elements z ∈ Z for
which αz ≤ 1. Formally, let ZH := { z ∈ Z : αz > 1 } and ZL := { z ∈ Z : αz ≤ 1 }, so that
X

δz2 min(αz , αz4 ) =

z∈Z

X

δz2 αz +

X

δz2 αz4 .

(8)

z∈ZL

z∈ZH

By Claim 3.1, at least one of the following holds: (1) z∈ZH δz αz > mε0 or (2) z∈ZL δz αz > mε0 . We
analyze each case separately.
P
(1) Suppose that z∈ZH δz αz > mε0 . We claim that E[A] > γ·mε0 2 . Indeed, this follows from Lemma 3.1
and (8) using the following chain of (in)-equalities:
P

P

P
X

δz2 αz ≥

z∈ZH

P

z∈ZH

δz αz

2

z∈ZH αz

2

> mε0 ,

where the first inequality is Cauchy–Schwarz, and the second follows using that
(2) Suppose that z∈ZL δz αz > mε0 . We claim that E[A] > γ·
and (8) using the following chain of (in)-equalities:
P

X
z∈ZL

δz2 αz4 ≥

m4 ε0 4
8n3

13

z∈ZH

αz ≤ m.

. Indeed, this follows from Lemma 3.1

4
1 m4 ε0 4
1  X
δ
α
>
.
z
z
8n3 z∈Z
8 n3
L

P

The first inequality essentially follows by an application Jensen’s inequality as follows: Let δ :=
P
2/3
z∈ZL δz . By Jensen’s inequality we have:
 X

(δz2/3 /δ) · δz1/3 αz

4

X

≤

z∈ZL

(δz2/3 /δ) · δz4/3 αz4 ,

z∈ZL

or
 X

δz αz

4

≤ δ3

z∈ZL

X

δz2 αz4 .

z∈ZL

Since δz ≤ 2 for all z ∈ Z, it follows that δ ≤ 2n, which completes the proof of the claim.
We have thus far established that
m4 ε0 4
E[A] > γ min mε ,
8n3
02

!

,

giving the first inequality of Proposition 3.1 (b). To complete the proof of the proposition, it suffices to show
that
!
q
4 04
02 m ε

min mε ,
min(n, m) .
n3
We show this below by considering the following cases:

√
• If n ≥ βm, we must be in the range 1/n1/8 ≤ ε0 ≤ 1 where max( n/ε0 2 , min(n7/8 /ε0 , n6/7 /ε0 8/7 )) =
√
4 04
4 04
4 04
n6/7 /ε0 8/7 . We get mnε3 ≤ β 3 mnε3 ≤ mε0 4 ≤ mε0 2 , and then since mnε3 ≥ β 7/2 m by our choice
of m in Eq. (5), we get that
m 4 ε0 4
min mε ,
n3

!

02

≥

q

min(n, m) ,

as desired assuming that β ≥ 1.
• If βm ≥ n, we must be in the range 0 < ε0 ≤ 1/n1/8 , and therefore min(n7/8 /ε0 , n6/7 /ε0 8/7 ) =
√
√
4 04
n7/8 /ε0 . Since mε0 2 ≥ β n and mnε3 ≥ β 4 n by our choice of m in Eq. (5), we get that
m 4 ε0 4
min mε ,
n3
02

!

≥

q

min(n, m) ,

as desired assuming that β ≥ 1.
This completes the proof of Proposition 3.1.
3.1.2

Analyzing the Variance of A

The main result of this subsection is the following proposition establishing an upper bound on the variance
of A as a function of its expectation:
Proposition 3.2. We have that
Var[A] ≤ C 00 (min(n, m) + E[A]) ,
for some absolute constant C 00 .
14

(9)

The rest of this subsection is devoted to the proof of Proposition 3.2. By the law of total variance, we
have that:
Var A = E[Var[A | σ]] + Var E[ A | σ ] .
We will proceed to bound each term from above, which will give the proof. We start with the first term.
Conditioned on σz = |Sz |, Eq. (4) gives that
Var[Az | σz ] ≤

Cσz2

!

δz2
1
+ 2
σz
σz

1{σz ≥4} = C (1 + E[ Az | σz ]) 1{σz ≥4} .

Therefore, for σ := (σz )z∈Z , we can write
Var[A | σ] ≤ C (min(n, M ) + E[ A | σ ]) ,
where we used the inequality
ately get

P

z∈Z

1{σz ≥4} ≤

P

z∈Z

(10)

1{σz ≥1} ≤ min(n, M ). From Eq. (10), we immedi-

E[Var[A | σ]] ≤ C (min(n, m) + E[A]) ,
as desired.
P
We now proceed to bound the second term. As shown in Lemma 3.1, E[ A | σ ] = z∈Z σz δz2 1{σz ≥4} .
By the independence of the σz ’s, we obtain that
X

Var [E[ A | σ ]] =

δz4 Var[σz 1{σz ≥4} ] .

(11)

z∈Z

From (11) and Claim 2.1, recalling that δz ≤ 2, z ∈ Z, we get that
Var [E[ A | σ ]] ≤ 4C 0

X

h

i

δz2 E σz 1{σz ≥4} = 4C 0 E[A] .

z∈Z

This completes the proof of Proposition 3.2.
3.1.3

Completing the Proof
4 04

Recall that the threshold of the algorithm is defined to be τ := γ2 max(mε0 2 , mnε3 ).
In the completeness case, by Proposition 3.1 (a), we have that E[A] = 0. Proposition 3.2 then gives that
Var[A] ≤ C 00 · min(n, m). Therefore, by Chebyshev’s inequality we obtain
Pr[ A ≥ τ ] ≤

1
Var[A]
4
min(n, m)
≤ 2 C 00
≤ ,
4 ε0 4
2
2
m
τ
γ
3
max(mε0 , 3 )
n

where the last inequality follows by choosing the constant β to be sufficiently large (compared to γ, C 00 ).
In the soundness case, by Chebyshev’s inequality we get:
Var[A]
Pr[ A < τ ] ≤ Pr[ |A − E[A]| ≥ E[A] /2 ] ≤ 4
≤ 4C 00
E[A]2

min(n, m)
1
+
2
E[A]
E[A]

!

≤

1
,
3

where the third inequality uses Proposition 3.2 and the fourth inequality uses Proposition 3.1 (b), assuming
β is sufficiently large. This completes the proof of correctness.

15

4

Estimating a Polynomial of a Discrete Distribution

In this section, we consider the following general problem: Given a degree-d n-variate polynomial Q ∈
Rd [X1 , . . . , Xn ] and access to i.i.d. samples from a distribution p ∈ ∆([n]), we want to estimate the quantity
Q(p) = Q(p1 , . . . , pn ) to within an additive error ε. In this section, we analyze an unbiased estimator for
Q(p) and provide quantitative bounds on its variance.
The structure of this section is as follows: In Section 4.1, we describe the unbiased estimator and establish its basic properties. In Section 4.2, we bound from above the variance of the estimator. Finally, Section 4.3 applies the aforementioned results to the setting that is relevant for our conditional independence
tester.

4.1

Unbiased Estimator and its Properties

We start by noting that the general case can be reduced to the case that the polynomial Q is homogeneous.
Remark 4.1 (Reduction to homogeneous polynomials). It is sufficient to consider, without loss of generality,
the case where Q ∈ Rd [X1 , . . . , Xn ] is a homogeneous polynomial, i.e., a sum of monomials of total
degree exactly d. This is because otherwise one can multiply any monomial of total degree d0 < d by
0
P
P
( ni=1 Xi )d−d : since ni=1 pi = 1, this does not affect the value of Q(p).
We henceforth assume Q is a homogeneous polynomial of degree d. Before stating our results, we will
need to set some notation. Given a multi-set S of independent samples from a distribution p ∈ ∆([n]), we let
P
ΦS denote the fingerprint of S, i.e., the vector (ΦS,1 , . . . , ΦS,n ) ∈ Nn of counts: ni=1 ΦS,i = |S|, and ΦS,i
is the number of occurrences of i in S. Moreover, for a vector α = (α1 , .. . , αn ) ∈ Nn , we write X α for
P
Q
for the multinomial coefficient
the monomial X α := ni=1 Xiαi , kαk for the `1 -norm ni=1 αi , and kαk
α
kαk!
α1 !···αn ! . Finally, for any integer d ≥ 0, we denote by Hd ⊆ Rd [X1 , . . . , Xn ] the set of homogeneous
degree-d n-variate polynomials.
The estimators we consider are symmetric, that is only a function of the fingerprint ΦS . We first focus
on the special case N = d.
Lemma 4.1. There exists an unbiased symmetric linear estimator for Q(p), i.e., a linear function Udd : Rd [X1 , . . . , Xn ] →
Rd [X1 , . . . , Xn ] such that
h
i
E Udd Q(ΦS ) = Q(p) ,
where S is obtained by drawing d independent samples from p.
Proof. For any ordered d-tuple T ∈ [n]d , by independence of the samples in S, we see that Pr[ S =
T]=

ΦT . For any α ∈ Nn with kαk = d, the number of T ∈ [n]d with fingerprint α is d . Thus,
p
=
p
T
i=1
i
α


Q
we have that Pr[ ΦS = α ] = αd pα . Noting that since kαk = kΦS k, ni=1 ΦαS,i
=
δ
,
we
can
define
Φ
,α
S
i
Qn

Udd X α (ΦS )
h

:=

d
α

!−1

1{ΦS =α} =

d
α

!−1 n
Y
i=1

ΦS,i
αi

!

.

(12)

i

Then we have E Udd X α (ΦS ) = pα . We extend this linearly to all Q ∈ Hd . By linearity of expectation, we
obtain an unbiased estimator for any such Q(p).
We can generalize this to N ≥ d, by taking the average over all subsets of size d of S of the above
estimator.
16

Proposition 4.1 (Existence). For N ≥ d and Q ∈ Hd written in terms of monomials as Q(X) =
the symmetric linear estimator
UN Q(ΦS ) :=

X

cα

α

!−1 n
Y

N
α, N − kαk

Xi
αi

i=1

P

α cα X

α,

!

(13)

is an unbiased estimator for Q(p).
Proof. For the case N = d, this follows from Lemma 4.1 and Eq. (12). For any set of d indices I ⊆
[N ], |I| = d, the subset SI = { Si : i ∈ I } is a set of d independent samples from p, thus Udd Q(ΦSI ) is
an unbiased estimator for Q(p). To get a symmetric unbiased estimator (and to reduce the variance), we can
take the average over all subsets SI of S of size d. We claim that this estimator is UN as defined above, i.e.,
that
!−1
X
N
UN Q(ΦS ) =
Udd Q(ΦS 0 ) .
(14)
d
0
0
S ⊆S,|S |=d
By linearity of expectation, the RHS is an unbiased estimator for Q(p), and so (14) suffices to show the
proposition. By linearity of UN and Udd , we need to show (14) for each
monomial X α . Noting that the
Qn
ΦS,i 
0
number of subsets S of S of size d that have fingerprint α is i=1 αi , we have
N
d

!−1
X

Udd X α (ΦS 0 )

=

S 0 ⊆S,|S 0 |=d

N
d

!−1
X
S 0 ⊆S,|S 0 |=d

d
α

!−1

1{ΦS0 =α}

!−1

=

=

N
α, N − kαk

X

!−1 n
Y

N
α, N − kαk

1{ΦS0 =α}

S 0 ⊆S,|S 0 |=d

i=1

ΦS,i
αi

!

α

= UN X (ΦS ) .
This completes the proof.
Proposition 4.2 (Uniqueness). The unbiased estimator UN Q(ΦS ) of (13) is unique among symmetric estimators. That is, for every N ≥ d, for any symmetric estimator VN : [n]N → R satisfying E[VN (ΦS )] =
Q(p), where S is a multiset of N samples drawn from p, one must have VN (ΦS ) = UN Q(ΦS ) for all S.
Proof. We first show that it is sufficient to establish uniqueness only for the case d = N , i.e., to show
that Udd maps polynomials to singletons. To argue this is enough, suppose N > d, and with have two
different N -sample estimators VN , WN for a homogeneous degree-d polynomial Q. Considering R :=
P
( ni=1 Xi )N −d Q which is homogeneous of degree N and agrees with Q on every probability distribution
p, we obtain two different N -sample estimators VN , WN for a homogeneous degree-N polynomial.
N from polynomials to estimators that gives an unbiased estimator for
When N = d, we have a map UN
P
the polynomial. By (12), for Q(X) = α cα X α , this is given by
Udd Q(ΦS )

=

X

cα

α

17

d
α

!−1

1{ΦS =α} .

Any symmetric estimator on N samples can by written as a linear combination of 1{ΦS =α} . Hence, given
an estimator VN , we can find a unique polynomial QVn with Udd QVn (ΦS ) = Vn by choosing cα to match the
coefficients in this linear combination, i.e., Udd is a bijection between polynomials and symmetric estimators.
Thus, if we have two different N -sample estimators VN , WN for a homogeneous degree-N polynomial Q,
at least one of them is Udd R for some homogeneous degree-N polynomial R.
Now we have an estimator VN that is unbiased for two different homogeneous degree-N polynomials
Q and R. So we get that for every p ∈ ∆([n]), Q(p) = ES [VN (ΦS )] = R(p). Hence, their difference D :=
Q − R is a non-zero homogeneous degree-N polynomial which vanishes on every point (x1 , . . . , xn ) ∈ Nn
P
with ni=1 xi = 1. By homogeneity, for every non-zero x = (x1 , . . . , xn ) ∈ Rn+ ,
D(x) =

kxkd1 D

x
kxk1

!

= kxkN
1 ·0=0,

and therefore D vanishes on the whole non-negative quadrant Rn+ = { x ∈ Rn : xi ≥ 0 for all i }. Being
identically zero on an open set, D must be the zero polynomial, leading to a contradiction.
The above shows existence and uniqueness of an unbiased estimator, provided the number of samples
N is at least the degree d of the polynomial (in p) we are trying to estimate. The proposition below shows
this is necessary: if N < d, there is no unbiased estimator in general.
Proposition 4.3. Let Q ∈ Hd be a homogeneous n-variate polynomial such that nk=1 X does not divide
Q. Then, there exists no unbiased estimator for Q(p) from N samples unless N ≥ d.
P

Proof. Suppose by contradiction that, for such a Q ∈ Hd , there exists an unbiased estimator for Q(p) with
N < d samples. Then, since UN
N (with the notation of the proof of Proposition 4.2) is invertible, this
estimator is also an unbiased estimator for some homogeneous degree-N polynomial R ∈ HN . Therefore,
P
it is also an unbiased estimator for the degree-d homogeneous polynomial R0 := R · ( nk=1 Xk )d−N ∈ Hd .
P
But by Proposition 4.2, one must then have Q = R0 , which is impossible since nk=1 X does not divide
Q.

4.2

Bounding the Variance of the Unbiased Estimator

Having established existence and uniqueness of our unbiased estimator, it remains to bound its variance:
Theorem 4.2. Fix N ≥ d, and let UN : Rd [X1 , . . . , Xn ] → Rd [X1 , . . . , Xn ] be as above. Then, for every
Q ∈ Hd ,
!
ksk Q(p) 2
h
i
X
(N − d)!2
2
s d
Q
E (UN Q(ΦS )) =
p
,
(15)
s
dX
N !(N − 2d + ksk)! ni=1 si !
s∈Nn
ksk≤d

where the expectation is over S obtained by drawing N independent samples from p.
Proof.
In order to establish ithe identity, we first consider monomials: for α, β ∈ Nn , we will analyze
h
E UN X α (ΦS ))UN X β (ΦS )) , before extending it to Q ∈ Hd , relying on the linearity of UN . First, note
that by definition of UN (in Eq. (13)),
α

β

UN X UN X =

1

n
Y



N
N
α,N −kαk β,N −kβk i=1

We will use the following fact:
18

Xi
αi

!

!

Xi
.
βi

(16)

Claim 4.1. For 0 ≤ a, b ≤ n, we have
n
a

!

n
b

min(a,b)

!

!

n
a+b−s

X

=

s=0

!

a+b−s
.
a − s, b − s, s

n n

Proof. The left-hand-side a b is the number of subsets A,B of [n] with |A| = a,|B| = b. We can group
the set of such pairs of subsets by the size of their intersection s = |A ∩ B|. Summing the size of these
classes gives that
n
a

!

n
b

min(a,b)

!

X

=

s=0

!

n
a − s, b − s, s, n − (a + b) + s

,

which is easily seen to be equivalent to the claim by multiplying out the factorials.
1

We can then rewrite, combining Eq. (16) and Claim 4.1, and setting πα,β :=

(α,NN−kαk)(β,NN−kβk)

for

convenience, that
α

β

UN X UN X = πα,β

n min(α
Y
Xi ,βi )
i=1

s=0

X

= πα,β

!

Xi
αi + βi − s

n
Y

i=1
s∈Nn
s≤min(α,β)

Xi
αi + βi − si

!

αi + βi − s
αi − s, βi − s, s
!

!

αi + βi − si
.
αi − si , βi − si , si

Taking the expectation over an N -sample multiset S, we obtain
h

α

i

β

X

E UN X (ΦS )UN X (ΦS ) = πα,β

n
Y

i=1
s∈Nn
s≤min(α,β)

n
Y
αi + βi − si
ΦS,i
E
αi − si , βi − si , si
αi + βi − si
i=1

! "

!#

.

Recalling the proof of Proposition 4.1, we have
" n
Y

E

i=1

ΦS,i
αi + βi − si

!#

!

=

h
i
N
E UN X α+β−s (ΦS )
α + β − s, N − (kαk + kβk − ksk)

=

N
pα+β−s ,
α + β − s, N − (kαk + kβk − ksk)

!

leading to
E[UN X α (ΦS ) UN X β (ΦS )

i
n
Y
N
αi + βi − si
α + β − s, N − (kαk + kβk − ksk) i=1 αi − si , βi − si , si

!

= πα,β

X

p

α+β−s

p

α+β−s

s∈Nn
s≤min(α,β)

= πα,β

X
s∈Nn
s≤min(α,β)

N
α − s, β − s, s, N − (kαk + kβk − ksk)

19

!

.

!

To get a better hold on this expression and extend the analysis to general homogeneous polynomials (instead
of monomials), we first massage the expression above under the additional constraint that kαk = kβk = d.
E[UN X α (ΦS ) UN X β (ΦS )

n
n
Y
(N − d)!2 Y
αi !
βi !
N!
i=1
i=1

=

(N − d)!2
N!

=

Recalling that

i

dksk X α
dX s

=

X
s∈Nn
s≤min(α,β)

ps

s∈Nn
s≤min(α,β)

pα+β−s
(N − 2d + ksk)! i=1 (αi − si )!(βi − si )!si !
Qn

n
n
Y
Y
pα+β−s
βi !
αi !
Qn
.
(N − 2d + ksk)! i=1 si ! i=1 (αi − si )! i=1 (βi − si )!

αi −si
αi !
i=1 (αi −si )! Xi

Qn

X

for s ≤ α, we have

n
n
Y
Y
dksk pα dksk pβ
αi !
βi !
α+β−s
=
p
s
s
dX
dX
(αi − si )! i=1 (βi − si )!
i=1

for s ≤ min(α, β), from which
h

i

E UN X α (ΦS )UN X β (ΦS ) =

(N − d)!2
N!

X
s∈Nn
s≤min(α,β)

ps
dksk pα dksk pβ
Qn
.
(N − 2d + ksk)! i=1 si ! dX s dX s

By linearity of U and differentiation, this implies that, for any Q, R ∈ Hd ,
E[UN Q(ΦS )UN R(ΦS )] =

(N − d)!2 X
dksk Q(p) dksk R(p)
ps
Qn
.
s
s
s
!
N!
(N
−
2d
+
ksk)!
dX
dX
i
n
i=1
s∈N
ksk≤d

Choosing R = Q yields Eq. (15).
By the previous theorem, in order to analyze the variance Var [UN Q(ΦS )] = E (UN Q(ΦS ))2 −
E[UN Q(ΦS )]2 , one needs to bound the different terms of


h

2

E (UN Q(ΦS ))

i

=

d
X
X

p

s

h=0 s∈Nn
ksk=h

dh Q(p)
dX s

!2



d
X
(N − d)!2
Qn
=
Th (Q, p, d, N ) ,
N !(N − 2d + h)! i=1 si ! h=0

letting Th (Q, p, d, N ) denote the inner sum for a given 0 ≤ h ≤ d. Next, we provide some useful bounds
on some of these terms. We show that the first term will be mostly taken care of in the variance by the
subtracted squared expectation, E[UN Q(ΦS )]2 = Q(p)2 . This allows us to get a bound on the variance
directly:
Corollary 4.1. For h ≥ 0,
(N − h)! X s
Th (Q, p, d, N ) ≤
p
N!
s∈Nn

dh Q(p)
dX s

!2

1
Qn

i=1 si !

,

ksk=h

and so
Var UN Q(ΦS ) ≤

d
X
(N − h)! X
h=1

N!

s∈Nn
ksk=h

20

ps

dh Q(p)
dX s

!2

1
Qn

i=1 si !

.

Proof. We have
d−h
Y N − 2d + h + i h−1
Y
(N − d)!2
1
(N − h)!
=
≤
N !(N − 2d + h)!
N
−
d
+
i
N
−
j
N!
i=1
j=0

which gives the bound on Th (Q, p, d, N ). For h = 0, this gives
that T0 (Q, p, d, N ) ≤ Q(p)2 = E[UN UN Q(ΦS )]2

2
and so if we expand Var UN Q(ΦS ) = E UN UN Q(ΦS ) − E[UN UN Q(ΦS )]2 , the T0 (Q, p, d, N ) term is
at least cancelled by the −E[UN UN Q(ΦS )]2 .
In view of bounding the rest of the terms, let Q+ ∈ Hd denote the polynomial obtained from Q by
P
P
making all its coefficients non-negative: that is, if Q = kαk=d cα X α , then Q+ := kαk=d |cα | X α . Then,
we show the following:
Lemma 4.2. Fix any 0 ≤ g ≤ d. Then,
d
X

dh Q(p)
1
d +
.
2
Q
(p)
max
Th (Q, p, d, N ) = O
Ng
dX s
s:ksk≥g
h=g




Proof. For g as above, we have
d
X

d
X

1
Th (Q, p, d, N ) =
O
Nh
h=g
h=g


1
=O
Ng


1

 X
s∈Nn
ksk=h

Qn

i=1 si !

X
d
X
h=g

s∈Nn

1
Qn

i=1 si !

p

dh Q(p)
dX s

!2

s

!2

ps

dh Q(p)
dX s

.

ksk=h
h

A useful observation is that since Q is homogeneous of degree d, then so is X s ddXQs for every s. Consider a
term cα X α in Q; a term X α will appear in X s
n
Y

dh Q
dX s

if and only if α ≥ s, in which case this term will be

n
n
Y
Y
αi !
αi
α
α
si !
.
|cα | X = |cα | X
si
(α
−
s
)!
i
i
i=1
i=1
i=1

!

!

Therefore, summing over all s, we get
d
X
X
h=g

s∈Nn

1
Qn

i=1 si !

X

s

d
n
X
X
Y
dh (cα X α )
αi
α
=
|c
|
X
α
s
dX
si
i=1
h=g s≤α

ksk=h

!

≤

X

|cα | X

s≤α

α

n
Y
i=1

αi
si

!

= 2d |cα | X α ,

ksk=h

where the inequality is an abuse of notation, assuming X is a non-negative vector. For the last equality, we
relied on the facts that kαk = d and
n
XY
s≤α i=1

αi
si

!

=

n X
Y
i=1 s:s≤αi

21

αi
s

!

=

n
Y
i=1

2αi = 2kαk .

By linearity and the definition of Q+ , this yields
d
X
X

1

h=g s∈Nn
ksk=h

Qn

i=1 si !

dh Q(p)
≤ 2d Q+ (p) ,
dX s

ps

and thus
d
X

1
Th (Q, p, d, N ) = O
Ng
h=g




≤O

X
d
X

1

h=g s∈Nn
ksk=h

Qn

s

i=1 si !

p

dh Q(p)
dX s

!2

dh Q(p)
1
d +
.
2
Q
(p)
max
Ng
dX s
s:ksk≥g


This completes the proof.
Remark 4.3 (Poissonized case). Frequently, in distribution testing we analyze Poissonized statistics. That
is, instead of S being a set of N samples,
we consider
a set S of Poisson(N ) samples. In this case, ΦS,i is
i
hQ
Qn N
Si 
n
α
independent for different i’s and E i=1 αi = p
i=1 αi ! . Thus, we can define an unbiased estimator
Si
0 Q for a polynomial Q(p) by taking linear combinations of U 0 X α (Φ ) = N −kαk n
for UN
S
i=1 αi αi !. The
N
theory in the Poissonized setting is a little different: this estimator is not unique
and is unbiased
for any

 0 α
N > 0, including non-integral N and N < d. However, the expression for E UN
X (ΦS )2 is very similar,
(N −d)!2
and is obtained by an analogous proof. The difference is that we obtain a term N −h instead of N !(N
−2d+h)! .
The bound on the variance in Corollary 4.1 holds for the unbiased estimators in both the Poissonized and
non-Poissonized cases.

Q

4.3



Case of Interest: `2 -Distance between p and pX ⊗ pY

We now instantiate the results of the previous subsections to a case of interest, the polynomial Q corresponding to the `2 distance between a bivariate discrete distribution and the product of its marginals. In more detail,
for any distribution p ∈ ∆(X × Y), where |X | = `1 , |Y| = `2 , we let pΠ := pX ⊗ pY ∈ ∆(X × Y) be the
product of its marginals. Moreover, let Q ∈ R4 [X1,1 , X2,1 , . . . , X`1 ,1 , X`1 ,2 , . . . , X`1 ,`2 ] be the degree-4
(`1 `2 )-variate polynomial defined as
Q(X1,1 , . . . , X`1 ,`2 ) :=

`1 X
`2
X
i=1 j=1

2


Xi,j

XX
i0 6=i j 0 6=j

Xi0 ,j 0 −

X
i0 6=i

Xi0 ,j

X

Xi,j 0  .

(17)

j 0 6=j

An explicit expression for its unbiased estimator UN Q(ΦS ) will be given in Eq. (18). Specifically, we shall
prove the following result:
Proposition 4.4. Let Q be as in Eq. (17), and suppose that b ≥ max(kpk22 , kpΠ k22 ). Then, for N ≥ 4,
!
√
Q(p) b
b
Var [UN Q(ΦS )] = O
+ 2 .
N
N
For consistency of notation with the previous section, we let n := `1 `2 in what follows.
22

Claim 4.2. For any p over X × Y, we have Q(p) = kp − pΠ k22 .
Proof. Unraveling the definitions, we can write
kp − pΠ k22 =

`1 X
`2
X

(p(i, j) − pΠ (i, j))2 =

i=1 j=1

=

`1 X
`2
X

=

`2
`1 X
X

`2
X

p(i, j) −

`1
X

X

2

p(i0 , j)

i0 =1

2



p(i, j) − p(i, j) +

0

p(i, j ) p(i, j) +

j 0 6=j



p(i, j 0 )

j 0 =1



X

0

p(i , j)

i0 6=i



2



p(i, j) 1 − p(i, j) −

X

p(i0 , j) −

i0 6=i

i=1 j=1

=



i=1 j=1



i=1 j=1
`1 X
`2
X

`1 X
`2
X

X

p(i, j 0 ) −

j 0 6=j

X

p(i0 , j)

i0 6=i

X

p(i, j 0 )

j 0 6=j

2


p(i, j)

XX

p(i0 , j 0 ) −

i0 6=i j 0 6=j

i=1 j=1

X

p(i0 , j)

i0 6=i

X

p(i, j 0 ) = Q(p) ,

j 0 6=j

as claimed.
Firstly, we compute UN Q explicitly. By linearity of UN , we can compute the unbiased estimator for each
P`1 P`2
P
P
2
term separately, after writing Q(X) = i=1
i0 6=i
j 0 6=j Xi0 ,j 0 −
j=1 ∆ij (X) , where ∆ij (X) := Xi,j
P
P
P`1 P`2
2
2
i0 6=i Xi0 ,j
j 0 6=j Xi,j 0 . Now UN Q =
i=1
j=1 UN ∆ij and we want to compute UN ∆ij . Note that
the sums in ∆ij (X) are over disjoint sets of Xi,j ’s whose union is every Xi,j . We can consider ∆ij as a
polynomial over the probabilities of a distribution with support of size 4, which consists of the events given
by whether the marginal X is equal to i, and whether the marginal Y is equal to j. By uniqueness of the
unbiased estimator, UN ∆2ij is the same on this distribution of support 4 as on the original `1 `2 -size support
distribution. Formally, we will write
∆ij (X) := Xi,j X−i,−j − Xi,−j X−i,j ,
where X−i,−j :=

P

i0 6=i

P

j 0 6=j

Xi0 ,j 0 , X−i,j :=

P

i0 6=i Xi0 ,j ,

and Xi,−j :=

P

j 0 6=j

Xi,j 0 . Squaring gives

2
2
2
2
∆ij (X)2 = Xi,j
X−i,−j
+ Xi,−j
X−i,j
− Xi,j X−i,−j Xi,−j X−i,j ,

and it remains to apply UN to each of these terms. We see that
N!
UN Xi,j X−i,−j Xi,−j X−i,j = ΦS,i,j ΦS,−i,−j ΦS,i,−j ΦS,−i,j ,
(N − 4)!
N!
2
2
UN Xi,j
X−i,−j
= ΦS,i,j (ΦS,i,j − 1)ΦS,−i,−j (ΦS,−i,−j − 1) ,
(N − 4)!
and

N!
2
2
UN X−i,j
Xi,−j
= ΦS,−i,j (ΦS,−i,j − 1)ΦS,i,−j (ΦS,i,−j − 1) .
(N − 4)!

23

These counts are similarly summed so that, for example, ΦS,i,−j =
we get that:
`

P

j 0 6=j

ΦS,i,j 0 . Adding these together,

`

1 X
2
X
N!
N!
UN ∆ij (ΦS )2
UN Q(ΦS ) =
(N − 4)!
(N − 4)! i=1 j=1

=

`1 X
`2
X

(ΦS,i,j (ΦS,i,j − 1)ΦS,−i,−j (ΦS,−i,−j − 1)

i=1 j=1

+ ΦS,−i,j (ΦS,−i,j − 1)ΦS,i,−j (ΦS,i,−j − 1) − 2ΦS,i,j ΦS,−i,−j ΦS,i,−j ΦS,−i,j )
=

`1 X
`2 
X

(ΦS,i,j ΦS,−i,−j − ΦS,−i,j ΦS,−i,−j )2

i=1 j=1

+ ΦS,i,j ΦS,−i,−j (1 − ΦS,i,j − ΦS,−i,−j ) + ΦS,−i,j ΦS,i,−j (1 − ΦS,−i,j − ΦS,i,−j )) ,
(18)
P

P

P

where ΦS,−i,−j := i0 6=i j 0 6=j ΦS,i0 ,j 0 , ΦS,−i,j := i0 6=i ΦS,i0 ,j , and ΦS,i,−j :=
the explicit formula for our unbiased estimator of Q(p).

P

j 0 6=j

ΦS,i,j 0 . This yields

We then turn to bounding its variance. From Theorem 4.2, we then have that, for N ≥ 4,
h

2

E (UN Q(ΦS ))

i

4
X
X

=

!

h=0 s∈Nn
ksk=h

h s
p
s

dh Q(p)
dX s

!2

N −4
4−h

!

N
h, 4 − h, N − 4

!−1

1
.
h!2

(19)

The rest of this section is devoted to bounding this
quantity. For h ∈ {0, . . . , 4}, we let Th (N ) be the inner
 P
sum corresponding to h, so that E (UN Q(ΦS ))2 = 4h=0 Th (N ).
P 1 P`2
2
For clarity, we (re-)introduce some notation: that is, we write Q(X) = `i=1
j=1 ∆ij (X) , where
P
P
P
P
∆ij (X) := Xi,j i0 6=i j 0 6=j Xi0 ,j 0 − i0 6=i Xi0 ,j j 0 6=j Xi,j 0 as before. Each ∆ij is a degree-2 polynomial,
with partial derivatives


Xi,j
if k 6= i, ` 6= j



P
P


if k = i, ` = j
∂∆ij
i0 6=i
j 0 6=j Xi0 ,j 0
=
P
∂Xk,` 
if k = i, ` 6= j
− i0 6=i Xi0 ,j




− P 0 X 0
if k 6= i, ` = j
j 6=j i,j
and

∂ 2 ∆ij
= (δik − δik0 )(δj` − δj`0 ) .
∂Xk,` ∂Xk0 ,`0
• The first contribution, for h = 0, is O Q(p)2 /N by Corollary 4.1, so we have T0 under control.
Indeed,
√
Q(p) ≤ 2 b


√
by the triangle inequality and the definition of b. So, T0 (N ) − Q(p)2 = O Q(p) b/N .


• The second, h = 1, contributes
T1 (N ) =

X
s∈Nn
ksk=1

s

p



dQ(p)
dX s

2

N −4
3

!

!−1

N
1, 3, N − 4

24

=4


2
N −4 X
s dQ(p)
3
p
N
dX s
s∈Nn
4
ksk=1

Since

N −4 N 
3 / 4

= O(1/N ), it is enough to consider the other factor,
X

p

dQ(p)
dX s



s

s∈Nn
ksk=1

2

=

X

dQ(p)
dXk,`

pk,`

k,`

!2

.

Recalling the expression of the derivatives of ∆ij , we have that
1 dQ
1X
d∆ij
=
2∆ij
2 dXk,`
2 i,j
dXk,`
=

XX

Xi,j ∆ij (X) + ∆k` (X)

i6=k j6=`

XX
i6=k j6=`

Having this sum of four terms A1 , A2 , A3 , A4 for
4(A21

A22

+
+
cross factors.

A23

+

Xi,j −

A24 ),

X

∆kj (X)

j6=`

dQ
dXk,` ,

X

Xi,j −

i6=k

X

∆i` (X)

i6=k

by Cauchy–Schwarz it holds that

X

Xi,j .

j6=`



2
dQ
dXk,`

≤

and so we can bound each of the square of these terms separately, ignoring

– For the first, we have (again by Cauchy–Schwarz)


2

2


√
X
X
XX
X

pi,j ∆ij (p) ≤  pi,j ∆ij (p) ≤  p2i,j   ∆ij (p)2  ≤ bQ(p) ≤ bQ(p) ,

so

P

k,` pk,`

P

i6=k

2

P

– For the second, since ∆k` (p)

pk,` ∆k` (p)

k,`

which is equal to

P

i6=k

P

j6=` pi,j

2

≤ ∆k` (p)2 , we have

2


X

≤ bQ(p).

j6=` pi,j ∆ij (p)



i,j

i,j

i,j

i6=k j6=`

XX

pi,j  ≤

i6=k j6=`

√

X

pk,` ∆k` (p)2 ≤

k,`

sX
k,`

bQ(p).

25

p2k,`

sX
k,`

v
2
u
√ u
X
u
∆k` (p)4 ≤ bt ∆k` (p)2  ,
k,`

– For the third and fourth term (similarly handled by symmetry),
X


2

2
X
X
X
X
X
pk,`  ∆kj (p)
pk,`  ∆kj (p)
pi,j  ≤
pi,j 

k,`

j6=`

j

k,`

i6=k

i6=k


2
X
X X
X

pk,`
=
∆kj (p)
pi,j 
j

k

`

i6=k



2 
X X
X X
X

≤
∆kj (p)2
pk,`
pi,j  



j

k

j

`

i6=k

(Cauchy–Schwarz)


≤

X

!2 
X



∆kj (p)2

X X

j

k

j

=

j

≤

pi,j

j

i

·

X

X

∆kj (p)2 

j

k

pk,`

`


X



i

XX

X




!2
X X

pi,j

pk,`

`

v
2 X  X
2
2 u
uX  X
∆kj (p)2
pi,j t
pk,`

i

j

k

k

`

(Cauchy–Schwarz)
v
v
uX  X
2 X  X
2 u
2
u
uX  X
pi,j
∆kj (p)2 ,
≤t
pk,` t
j

where the last step relies on

P

j

P

i

i pi,j

2

k

`

≤ 1 (since it is the squared `2 -norm of a proba-

bility distribution, that of the first marginal of p) to write

P

P

j

k,`


2
X
X
pk,`  ∆kj (p)
pi,j  ≤
j6=`

j

≤

≤

qP

2
i vi

≤

P

P

P

j

i |vi |

i pi,j
2

2

.

,

v
uX  X
2 X  X
2 X X
u
t
pi,j
∆kj (p)2
pk,`

i6=k

=

2

i pi,j

Continuing from there, and using monotonicity of `p norms to write
X

j

k

sX

i

pY (k)2

k

X

j

√

`

pX (j)2 Q(p) =

j

k

sX

k

pΠ (k, j)2 Q(p)

k,j

bQ(p) ,



√
and so T1 (N ) = O Q(p) b/N .

Gathering these four terms, and by the above discussion, we obtain
T1 (N ) = 4

N −4 X
3
pk,`
N
k,`
4

dQ(p)
dXk,`

!2

≤4

N −4
3
N
4

√
· 8 · 4 bQ(p) = O

√

!

bQ(p)
.
N

• Finally, for the rest of the contributions (h ≥ 2), we invoke Lemma 4.2. Specifically, we first observe
26

that, for any distribution p ∈ ∆(X × Y),
Q+ (p) =

`1 X
`2
X

2


pi,j

XX

pi0 ,j 0 +

i0 6=i j 0 6=j

i=1 j=1

X

pi0 ,j

i0 6=i

X

pi,j 0  ≤


X

j 0 6=j

pi,j +

`1
X

pi0 ,j

2

pi,j 0 

j 0 =1

i0 =1

i,j

`2
X

2 
`2
`1


X
X
X

2
pi,j 0   ≤ 2 kpk22 + kpΠ k22 ≤ 4b .
≤2
pi0 ,j  
pi,j + 


2 



j 0 =1

i0 =1

i,j

Next, we need to bound from above the high-order derivatives of Q. By Leibniz’s rule, for h ≥ 2 and
ksk = h, we can write:
0

0

n
XXY
dh Q X dh ∆2ij
s` dks k ∆ij dksk−ks k ∆ij
=
=
dX s
dX s
dX s0
dX s−s0
s0`
i,j
i,j s0 ≤s `=1

!

≤

X
s0 ≤s

v
!u
s` u
uX
t
0

n
Y

s`

i=`

i,j
0

≤ max
0

X

s ≤s

i,j

dks0 k ∆ij
dX s0

dks k ∆ij
dX s0

!2

!2
i,j

n
XY
s`
s0 ≤s

dks−s0 k ∆ij
dX s−s0

X

(Cauchy–Schwarz)
0

!

= 2h max
0

s0`

i=`

!2

X

s ≤s

i,j

dks k ∆ij
dX s0

!2

.

Since ∆ij has degree 2, to bound this maximum we have to consider three cases: first,
Q(p) ≤ 4. Second, recalling the partial derivatives of ∆ij we computed earlier,
X
i,j

d∆ij (p)
dXk,`

=

XX

p2k,` + 

Third,
X
i,j

d2 ∆ij (p)
dXk,` dXk0 ,`0

=

pi0 ,j 0  +

X

p2i0 ,` +

i0 6=k

X

p2k,j 0 ≤ 4 .

j 0 6=`

X

(δik − δik0 )2 (δj` − δj`0 )2 ≤ 4 .

i,j



Th (N ) = O

h=2


d0 ∆ij (p) 2
dX 0

!2

Combining all of the above cases results in
4
X

X X

i0 6=k j 0 6=`

i6=k j6=`



i,j

2



!2

P

dh Q
dX s

1
N2

≤ 24 · 4 for any h ≥ 2 and ksk = h, and from there


· 24 · 4 · 4b = O



b
.
N2


Accounting for all the terms, we can thus bound the variance as
!
√
b
Q(p) b
Var UN Q(ΦS ) = (T0 (N ) − Q(p) ) + T1 (N ) +
Th (N ) = O
+ 2 ,
N
N
h=2
4
X

2

concluding the proof of Proposition 4.4.
Remark 4.4 (Estimating a Polynomial under Poisson Sampling). We observe that analogues of our theorems
hold under Poisson sampling (instead of multinomial sampling as treated in Section 4). We defer these
results, which follow from a straightforward (yet slightly cumbersome) adaptation of the proofs of this
section, to an updated version of this paper.
27

=

5

The General Conditional Independence Tester

In this section, we present and analyze our general algorithm for testing conditional independence. The
structure of this section is as follows: In Section 5.1, we begin by describing how we flatten the marginals of
the distribution pz , for each bin z for which we receive enough samples. After this flattening is performed,
in Section 5.2 we explain how we use the remaining samples for each such bin z to compute a statistic A
as an appropriate weighted sum of bin-wise statistics Az . Before going further, we discuss in Section 5.3
the eventual result our analysis yields and comment on the sample complexity bound of our algorithm.
In Section 5.4, we explain the three different sources of randomness involved in our estimator, in order to
clarify what will follow – as we will crucially later condition on part of this randomness to obtain bounds on
some of its conditional expectations and variances. Section 5.5 then details how the analysis of our statistic
A is performed (Sections 5.5.1 and 5.5.2 respectively contain the analysis of the expectation and variance
of A, conditioned on some of the randomness at play). Finally, Section 5.6 puts everything together and
derives the correctness guarantee of our overall algorithm.

5.1

Flattening X , Y for any Given Bin z

Given a multiset S of N ≥ 4 independent samples from p ∈ ∆(X × Y), where |X | = `1 , |Y| = `2 , we
perform the following. Losing at most three samples, we can assume N = 4 + 4t for some integer t. Let
t1 := min(t, `1 ) and t2 := min(t, `2 ). We divide S into two disjoint multi-sets SF , ST of size t1 + t2 and
2t + 4 respectively, where the subscripts F and T stand for Flatten and Test.
• We use SF to flatten X × Y, as per Definition 2.2. Namely, first we partition it into two multi-sets SF1 ,
SF2 of size t1 , t2 . Looking at the projections πX SF1 , πY SF2 of SF1 , SF2 onto X and Y respectively, we
have two multi-sets of t1 and t2 elements. We then let T ⊆ X ×Y obtained by, for each (x, y) ∈ X ×Y,
adding in T ax,y copies of (x, y), where
1 + ax,y := (1 + ax )(1 + a0y )

(20)

with ax (resp. a0y ) being the number of occurrences of x in πX SF1 (resp. of y in πY SF2 ). Note that
|T | + `1 `2 = (|X | + t1 )(|Y| + t2 ), and that for all (x, y) ∈ X × Y, by a similar proof as that
of Lemma 2.3 (using the fact that ax and a0y are independent),
"

#

 "

#

1
1
1
1
E
=E
E
≤
,
1 + ax,y
1 + ax
1 + ay
(1 + t1 )(1 + t2 )pX (x)pY (x)


and so, letting qT denote the product of the marginals of pT ,
h

i

E kqT k22 ≤

1
.
(1 + t1 )(1 + t2 )

(21)

• Next, we use the 2t + 4 ≥ 4 samples from ST to estimate the squared `2 -distance between pT and
qT , as per Section 4. Here, Remark 2.3 will come in handy, as it allows us to do it implicitly without
having to actually map p to pT . Indeed, recalling that the polynomial Q for which we wish to estimate
Q(pT ) is of the form
X
Q(X) =
∆ij (X)2 ,
(i,j)∈X ×Y

28

we will instead estimate RT (p), where R is defined as
X

RT (X) :=

ci,j ∆ij (X)2

(i,j)∈X ×Y

with ci,j :=

1
1+ai,j

for all (i, j) ∈ X × Y. From Remark 2.3, it is immediate that RT (p) = Q(pT ) =

qT k22 ,

and further by inspection of the proof of Proposition 4.4 it is not hard to see that the
kpT −
variance of our estimator UN RT on p is the same as that of UN Q on pT .
Let B := kqT k22 . Note that B is a random variable, determined by the choice of SF . The first
observation is that, while the statement of Proposition 4.4 would be with regard to the maximum of
kpT k22 , kqT k22 , we would like to relate it to B. To do so, observe that




kpT k22 ≤ (kqT k2 + kpT − qT k2 )2 ≤ 2 kqT k22 + kpT − qT k22 = 2 (B + Q(pT ))
so we can use B 0 := 2B + 2Q(pT ) instead of our original bound B.
Therefore, our bound B can be used in the statement of Proposition 4.4, leading to a variance for our
estimator of
!
!
√
√
Q(pT ) B 0
Q(pT ) B Q(pT )3/2
B0
B
Var [UN RT ] = O
+ 2 =O
+
+ 2 .
(22)
N
N
N
N
N
Now, recall that by Lemma 2.3 (more precisely, Eq. (21)), we only have a handle on the expectation
of B. We could try to first obtain instead a high-probability bound on its value by proving sufficiently
strong concentration followed by a union bound over all estimators that we may run (i.e., all n bins in
Z). However, this would lead to a rather unwieldy argument. Instead, as outlined in Section 5.5, we
will analyze our estimators by carefully conditioning on some of the randomness (the one underlying
the flattening we perform for each bin), and only convert the bounds obtained into high-probability
statements at the end, by a combination of Markov’s and Chebyshev’s inequalities.

5.2

From Flattening to an Algorithm

We now explain how the guarantees established above are sufficient to use in our algorithm. We will use
the same notations as above, but now specifying the bin z ∈ Z: that is, we will write pz , qz , Tz , pz,Tz , qz,Tz
instead of p, q, T, pT , qT to make the dependence on the bin we condition on explicit. In what follows, we
write σ = (σz )z∈Z , T = (Tz | σz )z∈Z .
We let
Az := σz · ωz · Φ(Sz ) · 1{σz ≥4} ,
for all z ∈ Z, where ωz :=

p

min (σz , `1 ) min (σz , `2 ). Our final statistic is
A :=

X

Az .

z∈Z

That is, compared to algorithm of Section 3, we now re-weight the statistics by σz ωz instead of σz (since,
√
intuitively, the flattening is done with “t1,z , t2,z ” samples for which t1,z t2,z = Θ(ωz ) samples, we multiply
the weight by the “flattening amount”).

29

Recalling that `1 ≥ `2 without loss of generality, we set
1/4 1/4

2/7 2/7

1/2 1/2

!

1/2 1/2

1/2 3/2

n7/8 `1 `2 n6/7 `1 `2 n`1 `2
n3/4 `1 `2 `21 `22 n`1 `2
m ≥ ζ max min
,
,
,
min
, 4 ,
ε
ε
ε
ε
ε
ε8/7
!
!
√
√
√
2/3 1/3
3/2 1/2
n`1 `2 n`1 `2
n`1 `2 `1 `2 
n2/3 `1 `2 `1 `2
min
,
,
,
,
,
min
, 4
ε4
ε2
ε
ε2
ε
ε4/3


!

,
(23)

for some sufficiently big absolute constant ζ ≥ 1. The resulting pseudo-code is given in Algorithm 2.
Algorithm 2 T EST C OND I NDEPENDENCE G ENERAL
Require: Parameter n := |Z|, `1 := |X |, `2 := |Y|, ε ∈ (0, 1], and sample access to p ∈ ∆(X × Y × Z).
1: Set m as in Eq. (23)
. ζ ≥ 1 is an absolute constant
p
2: Set τ ← ζ 1/4 min(n, m).
. Threshold for accepting
3: Draw M ∼ Poisson(m) samples from p and let S be the multi-set of samples.
4: for all z ∈ Z do
5:
Let Sz ⊆ X × Y be the multi-set Sz := { (x, y) : (x, y, z) ∈ S }.
6:
if |Sz | ≥ 4 then
. Enough samples to call Φ
7:
Set Nz ← 4 b(|Sz | − 4)/4c, and let Sz0 be the multi-set of the first Nz elements of Sz .
.
Nz = 4 + 4tz for some integer tz .
8:
Set t1,z ← min(tz , `1 ), t2,z ← min(tz , `2 ), and divide Sz0 into disjoint SF0 ,z , ST0 ,z of size
t1,z + t2,z and σz := 2tz + 4, respectively.
(z)
9:
(ax,y )(x,y)∈X ×Y ← I MPLICIT F LATTENING(SF0 ,z ) . Flatten X × Y using SF0 ,z as explained in
the first bullet of Section 5.1, by calling Algorithm 3
(z)
10:
Φz ← U NBIASED E STIMATOR((ax,y )(x,y)∈X ×Y , SF0 ,z )
. Compute Φ(SF0 ,z ), the unbiased
estimator of Q as defined in the second bullet
p of Section 5.1, by calling Algorithm 4
11:
Set Az ← σz ωz · Φz , where ωz ← min (σz , `1 ) min (σz , `2 ).
12:
else
13:
Set Az ← 0.
14:
end if
15: end for
P
16: if A :=
z∈Z Az ≥ τ then
17:
return accept
18: else
19:
return reject
20: end if

5.3

Discussion of the Sample Complexity

The expression of our sample complexity in Eq. (23) may seem rather complicated. We argue here that it
captures at least some of the regimes of our four parameters in a tight way:
• For `1 = `2 = 2, we fall back to the case X = Y = {0, 1}, for which we had proven a tight bound
in Section 3. Note that in this case the expression of m in Eq. (23) reduces to



√
O max min(n7/8 /ε, n6/7 /ε8/7 ), n/ε2 ,
matching the bounds of Section 3.
30

Algorithm 3 I MPLICIT F LATTENING
Require: Multi-set S ⊆ X × Y.
1: . This simulates the construction of the “flattening set” as per Section 5.1; by Remark 2.3, it is actually
sufficient to compute the corresponding normalization coefficients ax,y , which we perform below.
2:
. All bx and cy are initialized to 0
3: for all (x, y) ∈ S do
4:
bx ← bx + 1
5:
cy ← cy + 1
6: end for
7:
. Note that the step below can be done more efficiently by only looping through elements (x, y) for
which either bx or cy is positive
8: for all (x, y) ∈ X × Y do
9:
ax,y ← (1 + bx )(1 + cy ) − 1
. Implement Eq. (20)
10: end for
11: return (ax,y )(x,y)∈X ×Y
Algorithm 4 U NBIASED E STIMATOR
Require: Set of coefficients (ax,y )(x,y)∈X ×Y , multi-set of samples S ⊆ X × Y.
1:
. This computes the unbiased estimator UN RT for Q(pT ) = RT (p) from the samples in S, as
explained in Section 5.1: where
RT (X) =

X
(x,y)∈X ×Y

1
∆x,y (X)2
1 + ax,y

2:

Let N ← |S|.

3:

. Recall that ΦS,x,y denotes the count of occurrences of (x, y) in the multi-set S
for all (x, y) ∈ X × Y do
. Compute for UN ∆x,y (ΦS )2 , from Eq. (18)
P
P
ΦS,−x,−y ← x0 6=x y0 6=y ΦS,x0 ,y0
P
ΦS,−x,y ← x0 6=x ΦS,x0 ,y
P
ΦS,x,−y ← y0 6=y ΦS,x,y0
Ci,j ← (ΦS,i,j ΦS,−i,−j − ΦS,−i,j ΦS,−i,−j )2 + ΦS,i,j ΦS,−i,−j (1 − ΦS,i,j − ΦS,−i,−j ) +
ΦS,−i,j ΦS,i,−j (1 − ΦS,−i,j − ΦS,i,−j )
end for
P
1
return (NN−4)!
(x,y)∈X ×Y 1+ax,y Ci,j
!

4:
5:
6:
7:
8:
9:
10:

• For n = 1 (and `1 ≥ `2 as before) we fall back to the independence testing
[BFF+ 01, LRR11, 
 problem

2/3 1/3 4/3 √
ADK15, DK16], for which the tight sample complexity is known to be Θ max `1 `2 /ε , `1 `2 /ε2 [DK16].



√
2/3 1/3
It is easy to see that, with these parameters, Eq. (23) reduces to O max `1 `2 /ε4/3 , `1 `2 /ε2
as well.
7/4
• For
 `1 =
 `2 = n (and ε not too small), the choice of m reduces to O(n /ε). This matches the
Ω n7/4 lower bound of Section 7 for ε = 1/20.
Remark 5.1. We further note that the expression of Eq. (23), which emerges from the analysis, can be simplified by a careful accounting of the regimes of the parameters. Namely, one can show that it is equivalent
31

to
1/4 1/4

2/7 2/7

n7/8 `1 `2 n6/7 `1 `2
m ≥ β max min
,
ε
ε8/7

5.4

!

1/2 1/2

2/3 1/3

1/2 1/2

n3/4 `1 `2 n2/3 `1 `2 n1/2 `1 `2
,
,
,
ε
ε2
ε4/3

!

.

(24)

The Different Sources of Randomness

As the argument will heavily rely on conditioning on some of the randomness at play and analyzing the
resulting conditional expectations and variances, it is important to clearly state upfront what the different
sources of randomness are and how we refer to them.
In what follows, we will use the following notations: for each bin z ∈ Z,
• σz is the number of samples from p we obtain with the Z coordinate falling in bin z;
• Tz is the randomness corresponding to the flattening of X , Y for the corresponding bin z (as described
in Section 5.1);
• Rz is the randomness of the estimator ΦSz on bin z.
Accordingly, we will write σ = (σz )z∈Z , T = (Tz )z∈Z , and R = (Rz )z∈Z for the three sources of
randomness (over all bins).

Analyzing A

5.5

The goal of this subsection is to show that, with high probability over σ, T , the following holds:
• If p is indeed conditionally independent, then E[ A | σ, T ] = 0 and Var[A | σ, T ] = O(min(n, m)).
p
• If p is far from conditionally independent, then E[ A | σ, T ] = Ω( min(n, m)) and Var[A | σ, T ] is
“not too big”compared to min(n, m) and E[ A | σ, T ].
This high-probability guarantee will allow us to use Chebyshev’s inequality in Section 5.6 to conclude
that, by comparing A to a suitably chosen threshold, we can distinguish between the two cases with high
probability (both over σ, T and R).
The reason for which we only obtain the above guarantees “with high probability over σ, T ” is, roughly
speaking, that we need to handle the complicated dependencies between A and T , which prevent us from
analyzing E[A] and Var[A] directly. To do so, we introduce an intermediate statistic, D (which itself only
depends on σ and R, but not on the flattening randomness T ), and relate it to E[ A | σ, T ]. This enables
us to analyze the expectation and variance of D instead of E[ A | σ, T ], before concluding by Markov’s
and (another application of) Chebyshev’s inequality that these bounds carry over to E[ A | σ, T ] with high
probability over σ, T .
5.5.1

The Expectation of A

We have that
E[ Az | σz , Tz ] = σz ωz kpTz − qTz k22 1{σz ≥4}
but Q(pTz ) = kpTz − qTz k22 depends on Tz . To get around this, we will start by analyzing D :=
where
ε2
Dz := σz ωz z 1{σz ≥4} ,
`1 `2

32

(25)
P

z∈Z

Dz ,

and εz := dTV (pz , qz ). Note that now D only depends on R and σ (and no longer on T ). For simplicity,
we will often write ε0z = √ε` z` . Next we show that whatever flattenings T = (Tz )z∈Z we use given
1 2
σ = (σz )z∈Z , D is a lower bound for the conditional expectation of A:
Lemma 5.1. E[ A | σ, T ] ≥ D.
Proof. Using that kpTz − qTz k2 ≥ √
D=

2εz
(`1 +t1,z )(`2 +t2,z )

≥

√εZ ,
`1 `2

we have that

X
ε2
σz ωz √ z 1{σz ≥4} ≤
σz ωz kpTz − qTz k22 1{σz ≥4} = A .
`
`
1 2
z∈Z
z∈Z
X

We will require the following analogue of Lemma 3.1 for D:
Lemma 5.2. For z ∈ Z, let αz := m · pZ (z). Then, we have that:
X

E[D] ≥ γ ·

ε2z min(αz βz , αz4 )

(26)

z∈Z

for some absolute constant γ > 0, where βz :=

p

min(αz , `1 ) min(αz , `2 ).

Proof. From the definition of D, we obtain that its expectation is:
E[D] =

X

h

i

Eσ σz ωz 1{σz ≥4} ε02
z .

z

Now
E[D] =

X

h

i

ε02
z Eσ σz ωz 1{σz ≥4} = Ω(1)

X

4
ε02
z min(αz βz , αz ) ,

z

z

h p

i

using the fact (Claim 2.3) that, for a Poisson random variable X with parameter λ, E X min(X, a) min(X, b)1{X≥4} ≥
γ min(λ min(λ, a) min(λ, b), λ4 ), for some absolute constant γ > 0.
p

We will leverage this lemma to show the following lower bound on the expectation of D:


Proposition 5.1. If dTV p, PX ,Y|Z
definition of m).








 p

> ε, then E[D] = Ω ζ min(n, m) (where ζ the constant in the

√mε
Proof. Since dTV p, PX ,Y|Z > ε, we have that z∈Z ε0z αz ≥ 2√1` `
z∈Z εz αz > `1 `2 .
1 2
We once again divide Z into heavy and light bins, ZH := z : αz3 ≥ βz and ZL := Z \ ZH .
P
P
By the above, we must have z∈ZH εz αz > mε or z∈ZL εz αz > mε. We proceed as in the proof
of Proposition 3.1 to handle these two cases.

P

P

• In the first case, we want to lower bound z∈ZH ε02
z αz βz . We consider three sub-cases, partitioning
ZH in 3: (1) ZH,1 := { z ∈ Z : `2 ≤ `1 < αz }, (2) ZH,2 := { z ∈ Z : `2 ≤ αz ≤ `1 }, and (3)
ZH,3 := { z ∈ Z : αz < `2 ≤ `1 }. By a similar argument, at least one of these sets is such that
P
1
0
z∈ZH,i εz αz > 3 mε .
P

33

– In the first sub-case:
P

ε02
z αz βz =

X

p

`1 `2

z∈ZH,1

ε02
z αz ≥

X

p

0
z∈ZH,1 εz αz

`1 `2

P

z∈ZH,1

z∈ZH,1

by Cauchy–Schwarz and recalling that z∈ZH,1 αz ≤
inequality after taking expectations on both sides,
P

ε02
z αz βz ≥

p

0
z∈ZH,1 εz αz

`1 `2

>

m

z∈ZH,1

P

P

≥

αz

p

`1 `2

0
z∈ZH,1 εZ αz

2

m

αz = m; and again by Jensen’s

z∈Z

2

P
X

2

p

`1 `2

1
1 mε2
2
√
mε0 =
.
36
36 `1 `2

(27)

– In the second sub-case:
2

P
X

ε02
z αz βz =

3/2
ε02
≥
z αz

X

p

`2

0
z∈ZH,2 εz αz
P
√
z∈ZH,2 αz

p

`2

z∈ZH,2

z∈ZH,2

2

P

≥

p

0
z∈ZH,2 εz αz
√
`2
√
min( mn, m/ `1 )

p
√ P
√
√
αz = m z∈ZH,2 pZ (z) ≤ mn, and
q
p
√ P
√
P
√
also that by definition of ZH,2 we have z∈ZH,2 αz ≤ m z∈ZH,2 pZ (z) ≤ m `m1 `m1 =
√m . Again by Jensen’s inequality after taking expectations on both sides,
`

by Jensen’s inequality and then using that

P

z∈ZH,2

1

2

P

X

ε02
z αz βz ≥

z∈ZH,2

p

0
p 1
z∈ZH,2 εz αz
m2 ε0 2
√
√
`2
> `2
√
√
36 min( mn, m/ `1 )
min( mn, m/ `1 )



1 m3/2 ε2
1
√ max √ ,
=
36 `1 `2
n

s



`1 
.
m

(28)

√ P
√
√
P
However, note that since z∈ZH,2 ε0z αz ≤ 2 z∈ZH,2 αz ≤ 2 |ZH,2 | `1 ≤ 2n`1 (as αz ≤
√
`1 for z ∈ ZH,2 ), the second sub-case cannot happen if mε0 ≥ 2 2n`1 .
– In the third sub-case:
P
X

ε02
z αz βz =

z∈ZH,3

X

2
ε02
z αz ≥

0
z∈ZH,3 εz αz

P

z∈ZH,3

z∈ZH,3

1

2

P

≥

0
z∈ZH,3 εz αz

2

min(n, m)

by Jensen’s inequality and recalling that |ZH,3 | ≤ min(n, m); and again by Jensen’s inequality
after taking expectations on both sides,
P
X
z∈ZH,3

ε02
z αz ≥

0
z∈ZH,3 εz αz

min(n, m)

2

!

1
m2
1
m2
ε2
2
>
ε0 =
max
,m
.
36 min(n, m)
36
n
`1 `2

(29)

√ P
√
√
P
However, note that since z∈ZH,3 δz αz ≤ 2 z∈ZH,3 αz ≤ 2 |ZH,3 | `2 ≤ 2n`2 (as αz <
√
`2 for z ∈ ZH,3 ), the third sub-case cannot happen if mε0 ≥ 2 2n`2 .

34

4
• In the second case, we want to lower bound z∈ZL ε02
z αz . We then use the same chain of (in)equalities as in the second case of Proposition 3.1, to obtain

P

P
X
z∈ZL

and recall that ε0z =
X

√εz
2 `1 `2

4
ε02
z αz ≥

z∈ZL

≤

√1
`1 `2

0
z∈ZL εz αz

4
ε02
z α z ≥ P

4


02/3 3

,

z∈ZL εz

to conclude

4
X
`1 `2  X 0 4 `1 `2  1
1 m 4 ε4
√
>
ε
α
.
=
ε
α
z
z
z
z
4n3 z∈Z
4n3 2 `1 `2 z∈Z
8 n3 `1 `2

(30)

L

L

√
√ P
√
P
However, note that since z∈ZL δz αz ≤ 2 z∈ZL αz ≤ 2 |ZL | ≤ 2n (as αz ≤ 1 for z ∈ ZL ),
√
the second case cannot happen if mε0 ≥ 2 2n.
p

It remains to use Eqs. (27) to (30) and our setting of m to show that E[D] ≥ C min(n, m) (where the
constant C > 0 depends on the choice of the constant in the definition of m).
√
• From Eq. (27) and the fact that m ≥ ζ min(`1 `2 /ε4 , n`1 `2 /ε2 ), we get
ε02
z αz βz 

X

q

ζ min(n, m)

z∈ZH,1

in the first sub-case of the first case.
√ √
3/2 1/2
2/3 1/3
• From Eq. (28) and the fact that m ≥ ζ min(n2/3 `1 `2 /ε4/3 , `1 `2 /ε4 , n`1 `2 /ε2 , n`1 `2 /ε),
we get
q
ε02
z αz βz 

X

ζ min(n, m)

z∈ZH,2





in the second sub-case of the first case (depending on whether min(n, m) min n, `m1 is equal to n2 ,
m2 /`1 , or mn). (The last term in the min enforcing the condition that this sub-case can only happen
whenever mε0 = O(n`1 ).)
1/2 1/2

1/2 3/2

• From Eq. (29) and the fact that m ≥ ζ min(n3/4 `1 `2 /ε, `21 `22 /ε4 , n`1 `2 /ε), we get
ε02
z αz βz 

X

q

ζ min(n, m)

z∈ZH,3

in the third sub-case of the first case (depending on whether
n3/2 /m2

p

min(n, m) min



1
n
m , m2



is equal to

1/m1/2 ).

or
(The last term in the min enforcing the condition that this sub-case can only
happen whenever mε0 = O(n`2 ).)
1/4 1/4

2/7 2/7

1/2 1/2

• From Eq. (30) and the fact that m ≥ ζ min(n7/8 `1 `2 /ε, n6/7 `1 `2 /ε8/7 , n`1 `2 /ε), we get
X

q

2
ε02
min(n, m)
z αz βz  ζ

z∈ZL

in the second case (depending on whether min(n, m) is equal to n or m). (The last term in the min
enforcing the condition that this sub-case can only happen whenever mε0 = O(n).)
This completes the proof of Proposition 5.1.
35

5.5.2

Variances of D and A

First we bound the variance of D:
Lemma 5.3.
Var[D] ≤ O(E[D]) .
Proof. Recall that D = z∈Z σz ωz ε02
z 1{σz ≥4} . Since the σz ’s are independent and Dz is a function of σz ,
the Dz ’s are independent as well and so
P

X

Var[D] =

Var[σz ωz ε02
z 1{σz ≥4} ] =

z∈Z

X

ε04
z Varσ [σz ωz 1{σz ≥4} ] .

z∈Z

As σz is distributed as Poisson(αz ), we can use Claim 2.2 to bound this
Var[D] ≤ C 0

X

ε04
z Eσ [σz ωz 1{σz ≥4} ]

z∈Z

≤C

0

X

ε02
z Eσ [σz ωz 1{σz ≥4} ]

z∈Z
0

= C E[D] ,
for some absolute constant C 0 > 0.
Since our statistic A is a linear combination of the ΦSz ’s and all the Sz ’s are independent by Poissonization, we get the analogue of Proposition 3.2:
Proposition 5.2. Let E :=

P

z∈Z

ωz2 BTz 1{σz ≥4} . Then,


Var[A | σ, T ] ≤ C E + E 1/2 E[ A | σ, T ] + E[ A | σ, T ]3/2



,

(31)

where E[ E | σ ] = O(min(n, M )) and C > 0 is some absolute constant.
Proof. Since Var[Az | σz , Tz ] = σz2 ωz2 1{σz ≥4} Var[Φ(Sz ) | σz , Tz ], we have by Eq. (22) that, for some
absolute constant C > 0,
Var[Az | σz , Tz ] ≤ C

σz2 ωz2

kpTz − qTz k22 BTz
BT
kpTz − qTz k32
+ 2z +
σz
σz
σz
p

!

!

1{σz ≥4}

.

We will handle the three terms of the RHS separately. First, by Cauchy–Schwarz and monotonicity of `p
norms we get that
X

σz ωz2 kpTz

−

qTz k22

q

!1/2
X

BTz 1{σz ≥4} ≤

z∈Z

ωz2 BTz 1{σz ≥4}

z∈Z

X

σz ωz kpTz −

2
qTz k22 1{σz ≥4}

!1/2

z∈Z

!1/2
X

≤

ωz2 BTz 1{σz ≥4}

z∈Z
1/2

=E

X

σz ωz kpTz − qTz k22 1{σz ≥4}

z∈Z

E[ A | σ, T ] ,

the last equality from Eq. (25).
36

(32)

B

Moreover, for the second term σz2 ωz2 σT2z 1{σz ≥4} = ωz2 BTz 1{σz ≥4} , it is immediate that summing over
z
P
all bins we get z∈Z ωz2 BTz 1{σz ≥4} = E.
Let us now turn to the last term of our upper bound on the variance. We can write
kpTz
σz2 ωz2

− qTz k32
1{σz ≥4} = σz ωz2 kpTz − qTz k32 1{σz ≥4} =
σz


≤ σz ωz kpTz − qTz k22 1{σz ≥4}

r

3/2

ωz 3/2 3/2
σz ωz kpTz − qTz k32 1{σz ≥4}
σz
= E[ Az | σ, T ]3/2

recalling that ωz ≤ σz by definition. We may use the inequality between the `1 and `3/2 norms to conclude
P
that z∈Z E[ Az | σ, T ]3/2 ≤ E[ A | σ, T ]3/2 , which leads by the above to
#

"

E

X

σz ωz2 kpTz

−

qTz k32 1{σz ≥4}

σ, T

. ≤ E[ A | σ, T ]3/2

(33)

z∈Z

Since the Az ’s are independent conditioned on σz and Tz , we have
Var[A | σ, T ] =

X

Var[A | σz , Tz ] ,

z∈Z

and therefore by Eqs. (32) and (33) and the definition of E we obtain




Var[A | σ, T ] ≤ O E 1/2 E[ A | σ, T ] + E + E[ A | σ, T ]3/2 .

(34)

It remains to establish the further guarantee that E[ E | σ ] = O(min(n, M )). To do so, observe that we can
write, as ωz only depends on the randomness σ,
h

i

E ωz2 BTz 1{σz ≥4} σ = ωz2 E[ BTz | σ ] 1{σz ≥4} ≤

ωz2
1
(1 + t1,z )(1 + t2,z ) {σz ≥4}

by Eq. (21). Recalling that ti,z = min((σz − 4)/4, `i ) by the definition of the flattening (Section 5.1) and
that ωz2 = min(σz , `1 ) min(σz , `2 ), this leads to
h

i

E ωz2 BTz 1{σz ≥4} σ ≤ O(1) · 1{σz ≥4} .
In particular, by summing over all bins z this implies that
E[ E | σ ] ≤ O(1)

X

1{σz ≥4} ≤ O(1)

z∈Z

X

1{σz ≥1} = O(min(n, M )) ,

(35)

z∈Z

as claimed.




Lemma 5.4 (Soundness). If dTV p, PX ,Y|Z > ε, then with probability at least 99/100 over σ, T we have
p

simultaneously E[ A | σ, T ] = Ω



ζ min(n, m) and



Var[A | σ, T ] ≤ O min(n, m) +

q

min(n, m)E[ A | σ, T ] + E[ A | σ, T ]

37

3/2



.

Proof. By Lemma 5.1, we have that D is a lower
pbound on E[ A | σ, T ] for all σ, T . Since Proposition 5.1
and Lemma 5.3 further ensures that E[D] ≥ Ω( ζ min(n, m)) and Var[D] ≤ O(E[D]), applying Chebyshev’s inequality on D results in
Var[D]
Pr E[ A | σ, T ] < κ ζ min(n, m) ≤ Pr [ D < O(E[D]) ] = O
σ,T
σ,T
E[D]2




q

1
=O
E[D]




1
=O p
ζ min(n, m)

!

!

≤

1
200

for some absolute constant κ > 0. This gives the first statement. For the second, we start from Eq. (31) and
we apply Markov’s inequality to E: as E[ E | σ ] = O(min(n, M )), with probability at least 399/400 we
P
have E ≤ 400E[ E | σ ] = O(min(n, M )). Moreover, recalling that M = z∈Z σz is a Poisson random
variable with parameter m, we have Pr[ M > 2m ] ≤ 399/400. Therefore, by a union bound


0

Pr Var[A | σ, T ] ≥ κ

σ,T



min(n, m) +

q

3/2



min(n, m)E[ A | σ, T ] + E[ A | σ, T ]

≤

1
1
1
+
=
400 400
200

again for some absolute constant κ0 > 0. This gives the second statement. A union bound over both events
concludes the proof.
Lemma 5.5 (Completeness). If p ∈ PX ,Y|Z , then with probability at least 99/100 over σ, T we have
simultaneously E[ A | σ, T ] = 0 and Var[A | σ, T ] ≤ O(min(n, m)).
Proof. The first statement is obvious by the definition of A as sum of the Az ’s, since εz = 0 for all z ∈ Z.
For the second, the proof is identical as that of Lemma 5.4, but having only to deal with the term E in the
bound on the variance (as the others are zero).

5.6

Completing the Proof

Let our threshold τ be set to ζ 1/4 min(n, m). Gathering the above pieces we obtain the following:
p

Lemma 5.6 (Soundness). If p is ε-far from conditionally independent, then Pr[ A < τ ] ≤

1
50 .

Proof. We apply Chebyshev’s inequality once more, this time to A0 := (A | σ, T ) and relying on the
bounds on its expectation and variance established in Lemma 5.4. Specifically, let E denote the event that
both bounds of Lemma 5.4 hold simultaneously; then
Pr[ A ≤ τ ] = Pr A0 ≤ τ ≤ Pr A0 ≤ τ E + Pr[E]








1
1
≤ Pr A0 − E[ A | σ, T ] ≥ E[ A | σ, T ] E +
2
100




where the second line is because, conditioned on E, E[ A | σ, T ] ≥ Ω

38

p



ζ min(n, m) ≥ ζ 1/4 min(n, m) =2τ .
p

It only remains to bound the first term:


0



Pr A − E[ A | σ, T ] ≥ E[ A | σ, T ] E ≤ O

min(n, m) +

p

min(n, m)E[ A | σ, T ] + E[ A | σ, T ]3/2
E[ A | σ, T ]2
p

min(n, m)
min(n, m)
1
=O
2 + E[ A | σ, T ] +
E[ A | σ, T ]
E[ A | σ, T ]1/2


1
1
≤ O 1/4 ≤
100
ζ
for the choice of a sufficiently large constant ζ in the definition of m.
Lemma 5.7 (Completeness). If p is conditionally independent, then Pr[ A ≥ τ ] ≤

!

1
50 .

Proof. Analogously to the proof in the soundness case, we apply Chebyshev’s inequality to A0 := (A | σ, T )
and relying on Lemma 5.5. Specifically, let E] denote the event that the bound of Lemma 5.5 holds; then
Pr[ A ≥ τ ] = Pr A0 ≥ τ ≤ Pr A0 ≥ τ E] + Pr[E 0 ] .








To conclude, we bound the first term:
Var[A0 | E 0 ]
min(n, m)
≤O
2
τ
τ2
 
1
1
≤
≤O
ζ
100
again for the choice of a sufficiently large constant ζ in the definition of m.


Pr A0 ≥ τ E 0 ] ≤


6





Sample Complexity Lower Bounds: The Case of Constant |X |, |Y|

In this section, we prove our tight sample complexity lower bound of


√
Ω max(min(n6/7 /ε8/7 , n7/8 /ε), n/ε2 )
for testing conditional independence in the regime that X = Y = {0, 1} and Z = [n]. This matches the
sample complexity of our algorithm in Section 3, up to constant factors. In the main body of this section,
we prove each lower bound separately.
The following expression for the total variation distance will be useful in the analysis of the lower bound
constructions:
Fact 6.1. For any p ∈ ∆(X × Y × Z) for X = Y = {0, 1} and Z = [n]. we have that:
dTV (pz , qz ) = 2 |Cov[(X | Z = z), (Y | Z = z)]| = kpz − qz k2 .
Proof. We have the following:
2dTV (pz , qz ) = |pz (1, 1) − (pz (1, 0) + pz (1, 1)) · (pz (0, 1) + pz (1, 1))|
+ |pz (1, 0) − (pz (1, 0) + pz (1, 1)) · (pz (1, 0) + pz (0, 0))|
+ |pz (0, 1) − (pz (0, 1) + pz (0, 0)) · (pz (0, 1) + pz (1, 1))|
+ |pz (0, 0) − (pz (0, 1) + pz (0, 0)) · (pz (1, 0) + pz (0, 0))|
= 4 |pz (0, 0) · pz (1, 1) − pz (0, 1) · pz (1, 0)|
= 4 |Cov[(X | Z = z), (Y | Z = z)]| .

39

(36)

!





First Lower Bound Regime: Ω n6/7 /ε8/7 for ε >

6.1

1
n1/8

Assume that we are in the regime where
√
max(min(n6/7 /ε8/7 , n7/8 /ε), n/ε2 ) = n6/7 /ε8/7 ,
i.e., ε > 1/n1/8 . Suppose there is an algorithm for ε-testing conditional independence drawing m ≤
cn6/7 /ε8/7 samples from p, for some sufficiently small universal constant c > 0. Note that in this regime,
we have m  n, i.e., we can assume that m < c0 n for some small constant c0 > 0.
The yes-instance A pseudo-distribution p is drawn from the yes-instances as follows: Independently for
each value z ∈ [n], we set:
1
1
• With probability m
n , we select pZ (z) = m and pz (i, j) = 4 for all i, j ∈ {0, 1}. In other words, we
select uniform marginals for the conditional distributions pz (i, j).
ε
• With probability 1 − m
n , we select pZ (z) = n and pz (i, j) us defined by the 2 × 2 matrix:

– With probability 1/2,

16 24
24 36

!

1
100

– With probability 1/2,

36 24
24 16

!

1
100

:= Y1 ,
:= Y2

It is easy to see that the resulting distribution p satisfies
" n
X

#

m 1
m
pZ (z) = n
·
+ 1−
E
n m
n
z=1






ε
·
n



∈ [1, 1 + ε] ,

i.e., the marginal for Z has mass roughly 1 in expectation, and that p ∈ P{0,1},{0,1}|[n] .
The no-instance A pseudo-distribution p is drawn from the no-instances as follows: Independently for
each value z ∈ [n], we set
1
1
• With probability m
n , we set pZ (z) = m and pz (i, j) = 4 for all i, j ∈ {0, 1}.
ε
• With probability 1 − m
n , we set pZ (z) = n and pz (i, j) be defined by the 2 × 2 matrix:

– With probability 1/8,

6 24
24 46

!

1
100

– With probability 1/8,

46 24
24 6

!

1
100

– With probability 3/4,

26 24
24 26

!

1
100

:= N1 ,
:= N2 ,
:= N3 .

Similarly, we have that E[ nz=1 pZ (z)] ∈ [1, 1 + ε]. Furthermore, the expected total variation distance
P
between such a p and the corresponding q := nz=1 pZ (z)qz ∈ P{0,1},{0,1}|[n] is
P

1
m 1
m
n
·
·0+ 1−
2
n m
n
 

E[dTV (p, q)] =





ε
n



1 12
1 12
3 4
·
+ ·
+ ·
8 100 8 100 4 100



3
m
ε
1−
ε>
,
100
n
100


=

where
!

!

!


12
1 46 24
1 7 
7 3 ,
=
−
100
100 24 6
100 3

!


4
1 26 24
1 5 
5 5 .
=
−
100
100 24 26
100 5

40



h



Thus, E[dTV (p, q)] = Ω(ε), which by Lemma 2.2 implies that E dTV p, P{0,1},{0,1}|[n]

i

= Ω(ε).

The next claim shows that, for each z ∈ [n], the first three norms of the conditional distribution pz (i, j)
match, hence do not provide any information towards distinguishing between the yes- and no-cases. Therefore, we need to get at least 4 samples (X, Y, Z) with the same value of Z — that is a 4-collision with regard
to Z — in order to have useful information.
Notation. Given a 4-variable function R = R[X1 , X2 , X3 , X4 ] and a real 2 × 2 matrix M ∈ M2 (R), we
will denote R(M ) := R(M1,1 , M1,2 , M2,1 , M2,2 ).
We have the following:
Claim 6.1. Let Y1 , Y2 , N1 , N2 , N3 the probability matrices in the definition of the yes and no-instances. For
every 4-variable polynomial R ∈ R[X1 , X2 , X3 , X4 ] of degree at most 3, the following holds:
1
1
3
1
1
R(N1 ) + R(N2 ) + R(N3 ) = R(Y1 ) + R(Y2 ) .
8
8
4
2
2
Proof. The first crucial observation is that the associated matrices can be expressed in the form
(N1 , Y1 , N3 , Y2 , N2 ) = (A + kB)0≤k≤4 ,
where A, B are the following matrices:
!

!

1
6 24
,
A=
100 24 46

1 10
0
B=
100 0 −10

.

Therefore, for any function 4-variable function R (not necessarily a polynomial), we have


1
1
3
R(N1 ) + R(N2 ) + R(N3 ) −
8
8
4


1
1
R(Y1 ) + R(Y2 )
2
2
1
3
1
1
1
= R(N1 ) − R(Y1 ) + R(N3 ) − R(Y2 ) + R(N2 )
8
2
4
2
8
!
!
4
4
X
X
1
1
k 4
4−k 4
R(A + kB) =
R(A + kB) ,
=
(−1)
(−1)
k
k
8 k=0
8 k=0




which is the 4th -order forward difference of R at A (more precisely, the fourth finite difference of f (k) =
R(A + kB)). Using the fact that the (d + 1)th -order forward difference of a degree-d polynomial is zero,
we get that the above RHS is zero for every degree-3 polynomial R.
For the sake of simplicity and without loss of generality, we can use the Poissonization trick for the
analysis of our lower bound construction (cf. [DK16, CDKS17]). Specifically, instead of drawing m independent samples from p, we assume that our algorithm is provided with mz samples from the conditional
distribution pz (i.e., conditioned on Z = z), where the (mz )’s are independent Poisson random variables
with mz ∼ Poisson(mpZ (z)).
 
Consider the following process: we let U ∼ Bern 21 be a uniformly random bit, and choose p to be
selected as follows: (i) If U = 0, then p is drawn from the yes-instances, (ii) If U = 1, then p is drawn from
01 10 11
the no-instances. For every z ∈ [n], let az = (a00
z , az , az , az ) be the 4-tuple of counts of (i, j)i,j∈{0,1}
among the mz samples (X, Y ) ∼ pz . Accordingly, we will denote A = (az )z∈[n] .
41

Following
complexity
lower
 the mutual
 information method used in [DK16], to show the desired sample


6/7
8/7
6/7
8/7
bound of Ω n /ε
, it suffices to show that I (U ; A) = o(1), unless m = Ω n /ε
. Since the
Pn
(az )z∈[n] ’s are independent conditioned on U , we have that I (U ; A) ≤ z=1 I (U ; az ), and therefore it
suffices to bound from above separately I (U ; az ) for every z. We proceed to establish such a bound in the
following lemma:
Lemma 6.1. For any z ∈ [n], we have I (U ; az ) = O



ε8 m7
n7



.

Before proving the lemma, we show that it implies the desired lower bound. Indeed, assuming Lemma 6.1,
we get that
!
!
n
n
X
X
ε8 m7
ε8 m 7
I (U ; A) ≤
I (U ; az ) =
O
=O
,
n7
n6
z=1
z=1




which is o(1) unless m = Ω n6/7 /ε8/7 . It remains to prove the lemma.
Proof of Lemma 6.1. By symmetry, it is sufficient to show the claim for z = 1. To simplify the notation, let
a := a1 . We first bound I (U ; a) from above using [CDKS17, Fact 4.12] (see also [DK16]) as follows:
I (U ; a) ≤

Pr[ a = α | U = 1 ]
Pr[ a = α ] 1 −
Pr[ a = α | U = 0 ]
4

X



2

:= Φ(n, m, ε) .

α∈N

Our next step is to get a hold on the conditional probabilities Pr[ a = α | U = 0 ] and Pr[ a = α | U = 1 ].
For notational convenience, we set
16
24
36
6
26
46
p1 :=
, p2 :=
, p3 :=
, q1 :=
, q2 :=
, q3 :=
,
100
100
100
100
100
100
1
on it. Note that by
and let Ξ denote the event that the bin is “heavy”, i.e., that pZ puts probability mass m
m
construction this event happens with probability n and that Ξ is independent of U .
We start with the yes-case. Recall
 that with
 probability m/n, Ξ holds: the probability pZ (1) equals
1
1/m, in which case we draw Poisson m · m samples from Z = 1 and each sample is uniformly random


on {0, 1} × {0, 1}. That is, each of the four outcomes is an independent Poisson m ·
ε
n

1
m

·

1
4





random

variable. With probability 1 − m/n, Ξ̄ holds: pZ (1) equals ε/n and we draw Poisson m ·
samples from
Z = 1. With probability 1/2, all samples follow the first case, and with probability 1/2 all samples follow
the second case.
For any α = (α1 , α2 , α3 , α4 ) ∈ N4 , we can explicitly calculate the associated probabilities. Specifically,
we can write:
εm


h
i
e− n
1 α1 α2 α3 α4 1 α1 α2 α3 α4
Pr a = α U = 0, Ξ̄ =
p1 p2 p2 p3 + p3 p2 p2 p1
,
(37)
α1 !α2 !α3 !α4 ! 2
2
and
h

i

h

Pr[ a = α | U = 0 ] = Pr[ a = α | U = 0, Ξ ] · Pr[ Ξ ] + Pr a = α U = 0, Ξ̄ · Pr Ξ̄
P4

m −· 1 ·4 4− `=1 α`
m
=
·e 4
+ 1−
n
α1 !α2 !α3 !α4 !
n
− εm
1
e n
pα1 pα2 pα3 pα4 +
2 α1 !α2 !α3 !α4 ! 1 2 2 3
εm

1
e− n
pα3 1 pα2 2 pα2 3 pα1 4 .
2 α1 !α2 !α3 !α4 !


42

 

·

i

Similarly, for the no-case, we have
εm

e− n
Pr a = α U = 1, Ξ̄ =
α1 !α2 !α3 !α4 !
h

i



1 α1 α2 α3 α4 1 α1 α2 α3 α4 3 α1 α2 α3 α4
q p p q + q3 p2 p2 q1 + q2 p2 p2 q2
8 1 2 2 3
8
4



,
(38)

Pr[ a = α | U = 0, Ξ ] = Pr[ a = α | U = 1, Ξ ] ,
and
h

i

h

Pr[ a = α | U = 1 ] = Pr[ a = α | U = 1, Ξ ] · Pr[ Ξ ] + Pr a = α U = 1, Ξ̄ · Pr Ξ̄
P4

i

m − 1 ·4 4− `=1 α`
m
=
·
·e 4
+ 1−
n
α1 !α2 !α3 !α4 !
n
εm
εm
εm
1

e− n
e− n
e− n
1
3
q1α1 pα2 2 pα2 3 q3α4 +
q3α1 pα2 2 pα2 3 q1α4 +
q2α1 pα2 2 pα2 3 q2α4 .
8 α1 !α2 !α3 !α4 !
8 α1 !α2 !α3 !α4 !
4 α1 !α2 !α3 !α4 !




With these formulas in hand, we can write
Φ(n, m, ε) :=



X

Pr[ a = α ]

α∈N4

m
= 1−
n


2 X



≤

α∈N4

h

i

h

h

i

h

Pr a = α U = 0, Ξ̄ − Pr a = α U = 1, Ξ̄
Pr[ a = α | U = 0 ]
h

i 2


Pr[ a = α | U = 0 ]

α∈N4

Pr[ a = α ] 

2

Pr a = α U = 0, Ξ̄ − Pr a = α U = 1, Ξ̄

Pr[ a = α ] 


X

Pr[ a = α | U = 0 ] − Pr[ a = α | U = 1 ]
Pr[ a = α | U = 0 ]

i 2
 .
i

h

By Eqs. (37) and (38) and Claim 6.1, we observe that the difference Pr a = α U = 0, Ξ̄ −Pr a = α U = 1, Ξ̄

43

i

is zero for any |α| :=

P

i αi

≤ 3. We thus obtain


Φ(n, m, ε) ≤

=

X

h

i

Pr[ a = α ] 

α∈N4
|α|≥4
∞ X
X

h

Pr a = α U = 0, Ξ̄ − Pr a = α U = 1, Ξ̄

i 2


Pr[ a = α | U = 0 ]
h

Pr[ a = α ] · Pr |a| = k Ξ̄

i2

·

(†)

k=4 α∈N4
|α|=k



h

i



=



k!2

n

i

h

h

Pr a = α U = 0, Ξ̄, |a| = k − Pr a = α U = 1, Ξ̄, |a| = k

k!2

n



Pr[ a = α ] ·

α∈N4
|α|=k

i 2


Pr[ a = α | U = 0 ]


∞ − 2εm 
X
e n
εm 2k X
k=4

Pr[ a = α ] ·

α∈N4
|α|=k



≤

i 2


Pr[ a = α | U = 0 ]


∞ − 2εm 
X
e n
εm 2k X
k=4

h

Pr a = α U = 0, Ξ̄, |a| = k − Pr a = α U = 1, Ξ̄, |a| = k

2
Pr[ a = α | U = 0 ]

2

,

where for (†) we used the fact that, |a| is independent of U to write
εm

e− n
Pr |a| = k U = 1, Ξ̄ = Pr |a| = k U = 0, Ξ̄ = Pr |a| = k Ξ̄ =
k!
i

h

i

h

i

h



εm
n

k

.

To conclude the proof, we will handle the denominator using the bound
Pr[ a = α | U = 0 ] ≥ Pr[ a = α | Ξ, U = 0 ] · Pr[ Ξ | U = 0 ] =
and rewrite Pr[ a = α ] = Pr[ a = α | |a| = |α| ] · Pr[ |a| = |α| ]. Using x :=
sions for conciseness, we now get:

m −1 4−k
e
,
n
k!

εm
n

in the following expres-

∞
∞
2
X
X
n2 −2x X
2k
2 n
−2x
(4x) Pr[ |a| = k ]
Pr[ a = α | |a| = k ] = 4e 2 e
(4x)2k Pr[ |a| = k ]
Φ(n, m, ε) ≤ 4e 2 e
m
m
4
k=4
k=4
2

α∈N
|α|=k

∞
X
n2
= 4e 2 e−2x
(4x)2k
m
k=4
2

∞
n2 X
≤ 40 2
(4x)2k
m k=4

m −x xk
m e−1
+ 1−
e
n k!
n
k!


m1
xk
+
n k!
k!



!

= 40

!

∞
∞
n X
1
n2 X
1  2/3 3k
(4x)2k + 40 2
4 x
m k=4 k!
m k=4 k!

n (4x)8
n 8
n2 1  2/3 12
n2 12
= 40
+o
x + 40 2
4 x
+o
x
m 24
m
m 24
m2




ε8 m7
ε8 m7
=2
·
+
o
3
n7
n7
16 5

!

,

44

!

(‡)

where for (‡) we relied on the Taylor series expansion of exp, recalling that
proof of Lemma 6.1.

6.2





1
n3/8

Second Lower Bound Regime: Ω n7/8 /ε for

≤ε≤

εm
n

 1. This completes the

1
n1/8

Assume we are in the regime where
max(min(n6/7 /ε8/7 , n7/8 /ε),

√

n/ε2 ) = n7/8 /ε ,

1
1
≤ ε ≤ n1/8
. Suppose there is a testing algorithm for conditional independence using m ≤ cn7/8 /ε
i.e., n3/8
samples, for a sufficiently small universal constant c > 0. In this regime, we also have m  n.
Our construction of yes- and no- instances in this case is similar to those of the previous lower bound,
although some specifics about how pZ is generated will change.

The yes-instance. A pseudo-distribution p is drawn from the yes-instances as follows: Independently for
each value 1 ≤ z ≤ n − 1, we set:
1
• With probability 12 , pZ (z) = m
and pz (i, j) = 14 for all i, j ∈ {0, 1}. That is, we select uniform
marginals for pz (i, j).
• With probability 12 , pZ (z) = nε and pz (i, j) is defined by the same 2 × 2 matrices as in the previous
case:

– With probability 1/2, Y1 ,
– With probability 1/2, Y2 .
Furthermore, we set pZ (n) = 1, and pn (i, j) = 14 for all i, j ∈ {0, 1}. The last condition ensures that
kpZ k1 = Θ(1). It is clear that the resulting pseudo-distribution p satisfies kpZ k1 = Θ(1) and that p ∈
P{0,1},{0,1}|[n] .
The no-instance. A pseudo-distribution p is drawn from the no-instances as follows: Independently for
each value z ∈ [n], we set:
1
1
• With probability m
n , pZ (z) = m and pz (i, j) = 4 for all i, j ∈ {0, 1}, as before.
ε
• With probability 1 − m
n , pZ (z) = n , and pz (i, j) is defined by the same 2 × 2 matrix as before:

– With probability 1/8, N1 ,
– With probability 1/8, N2 ,
– With probability 3/4, N3 .
Furthermore,
we set as beforeipZ (n) = 1, and pn (i, j) =

h

that E dTV p, P{0,1},{0,1}|[n]

1
4

for all i, j ∈ {0, 1}. This construction ensures

= Ω(ε).





The bulk of the proof of the Ω n6/7 /ε8/7 remains the same. In particular, setting a := a1 , we can
bound from above as before the mutual information I (U ; a) by
I (U ; a) ≤


∞ − 2εm 
X
e n
εm 2k X
k=4

k!2

n



Pr[ a = α ] ·

α∈N4
|α|=k

45

2
Pr[ a = α | U = 0 ]

2

.

In the current setting, we use the fact that
1
4−k
Pr[ a = α | U = 0 ] ≥ Pr[ a = α | Ξ, U = 0 ] · Pr[ Ξ | U = 0 ] = e−1
,
2
k!
to obtain that
2 − 2εm
n

I (U ; a) ≤ 16e e


∞ 
X
4εm 2k X

n

k=4

2 − 2εm
n

Pr[ a = α ] = 16e e

α∈N4
|α|=k
εm

2k

ε8 m8
I (U ; a) = O
n8


8

8





8

ε m
Therefore, I (U ; A) ≤ n−1
= O ε nm7
z=1 O
n8
pletes the proof of this branch of the lower bound.

6.3

√ 

Third Lower Bound Regime: Ω

n

k=4

1 εm
1
Recalling that Pr[ |a| = k ] = 12 e−1 k!
and that
+ 12 e− n k!
n
expansion of the first term of the sum, we finally get that

P


∞ 
X
4εm 2k

n
ε2

8



εm
n



=Θ

1
n1/8



Pr[ |a| = k ] .

 1, with a Taylor series

!

.




, which is o(1) unless m = Ω n7/8 /ε . This com-

for ε < n−3/8

Finally, assume we are in the regime where
√
√
max(min(n6/7 /ε8/7 , n7/8 /ε), n/ε2 ) = n/ε2 ,
i.e., ε < n−3/8 . In this case, we can show the desired lower bound by a simple reduction from the
knownhard instances
for uniformity testing. Let N be an even positive integer. It is shown in [Pan08]

√
2
that Ω N /ε samples are required to distinguish between (a) the uniform distribution on [N ], and (b) a
distribution selected at random
by pairing

  consecutive elements 2i, 2i + 1 and setting the probability mass
1+2ε 1−2ε
1+2ε
of each pair to be either N , N
or 1−2ε
independently and uniformly at random.
N , N
We map these instances to our conditional independence setting as follows: Let N = 4n. We map
[N ] to the set {0, 1} × {0, 1} × [n] via the mapping Φ : [N ] → {0, 1} × {0, 1} × [n] defined as follows:
Φ(2i) = (0, 0, i), Φ(2i + 1) = (0, 1, i), Φ(2i + 2) = (1, 0, i), and Φ(2i + 3) = (1, 1, i).
For a distribution p ∈ ∆([N ]) adversarially selected as described above, the following conditions are
satisfied: (1) In case (a), Φ(p) is the uniform distribution on {0, 1} × {0, 1} × [n] and therefore Φ(p) ∈
P{0,1},{0,1}|[n] . (2) In case (b), it is easy to see that for each fixed value of the third coordinate, the conditional
distribution on the first two coordinates is one of the following:
!

1 1 + 2ε 1 − 2ε
,
4 1 + 2ε 1 − 2ε

!

1 1 + 2ε 1 − 2ε
,
4 1 − 2ε 1 + 2ε

!

1 1 − 2ε 1 + 2ε
, or
4 1 − 2ε 1 + 2ε

1 1 − 2ε 1 + 2ε
4 1 + 2ε 1 − 2ε

!

.

It is clear that each
is Ω(ε)-far from independent. Therefore, by Lemma 2.2 it
 of these four distributions

follows that dTV Φ(p), P{0,1},{0,1}|[n] = Ω(ε).
In conclusion, we have established that any testing algorithm for P{0,1},{0,1}|[n] can be used (with
parameter
ε)
uniformity over [N ] (with parameter ε0 = O(ε)). This implies a lower bound of
√
 to test

√
Ω N /ε0 2 = Ω n/ε2 on the sample complexity of ε-testing conditional independence.
46

7

Sample Complexity Lower Bound for X = Y = Z = [n]




In this section, we outline the proof of the (tight) sample complexity lower bound of Ω n7/4 for testing
conditional independence in the regime that X = Y = Z = [n], and ε = Ω(1).
Specifically, we will show that, when |X | = |Y| = |Z| = n and ε = 1/20, it is impossible to distinguish
between conditional independence and ε-far from conditional independence with s = o(n7/4 ) samples. To
do this, we begin by producing an adversarial ensemble of distributions. The adversarial distribution will be
designed to match the cases where the upper bound construction will be tight. In particular, each conditional
marginal distribution will have about n3/4 heavy bins and the rest of the bins light. The difference between
our distribution and the product of the marginals (if it exists) will be uniformly distributed about the light
bins.
First we will come up with an ensemble where X and Y are conditionally independent and then we
will tweak it slightly. For the first distribution, we let the distribution over Z be uniform. For each value of
Z = z, we pick random subsets Az , Bz ⊂ [n] of size n3/4 (which we assume to be an integer), which will
be the heavy bins of the conditional distribution. We then let the conditional probabilities be defined by
(

Pr[ X = j | Z = z ] =

n−3/4 /2
if j ∈ Az
3/4
1/(2(n − n )) else

and
(

Pr[ Y = j | Z = z ] =

n−3/4 /2
if j ∈ Bz
3/4
1/(2(n − n )) else.

We then let X and Y be conditionally independent on Z, defining the distribution over (X, Y, Z).
Finally, we introduce a new one bit variable W . In ensemble D0 , W is an independent random bit. In
ensemble D1 , the conditional distribution on W given (X, Y, Z) is uniform random if X ∈ AZ or Y ∈ BZ ,
but otherwise is given by a uniform random function f : [n] × [n] × [n] → {0, 1}. In particular, in this case,
W is determined by the values of X, Y, Z, though different elements of D1 will give different functions.
Note that elements of D0 have XW and Y conditionally independent on Z, whereas elements of D1 are
ε-far from any such distribution. We show that no algorithm that takes s samples from a random distribution
from one of these families can reliably distinguish which family the samples came from.
In particular, let F be a uniform random bit. Let S be a sequence of s quadruples (W, X, Y, Z) obtained
by picking a random element p from DF and taking s independent samples from p. It suffices to show that
one cannot reliably recover F from S. Note that with high probability S contains at most t := 2s/n +
log n = o(n3/4 ) samples for each value of Z. Therefore of we let Tz be a sequence of t independent
samples from p conditioned on Z = z for each z, it suffices to show that F cannot be reliably recovered
from T1 , . . . , Tn . For this it suffices to show that I (F ; T1 , . . . , Tn ) = o(1). Since the Tz are conditionally
Pn
independent on F , we have that I (F ; T1 , . . . , Tn ) ≤
z=1 I (F ; Tz ) and thus it suffices to show that
I (F ; Tz ) = o(1/n) for every z. Since this shared information is clearly the same for each z, we will
suppress the subscript.
We say that two distinct elements of T collide if they have the same values of X and Y . We note that if
we condition on the values of X and Y in the elements of T , that the values of W for the elements that do not
collide are uniform random bits independent of the other values of W and of F . Therefore, no information
can be gleaned from these W ’s. This means that all of the information comes from W ’s associated with
collisions. Unfortunately, most collisions (as we will see) come from heavy values of either X or Y (or
both), and these cases will also provide no extra information.
47

More formally, note that
"

Pr[ T = M | F = 0 ]
O min 1, 1 −
Pr[ T = M | F = 1 ]

I (F ; T ) = EM ∼T



2 !!#

.

We claim that if M has C pairs that collide, then


1−

Pr[ T = M | F = 0 ]
Pr[ T = M | F = 1 ]

2

= O(C 2 /n).

Our result will then follow from the observation that E C 2 = O(t2 /n3/2 ) = o(1).
Given M , call a value of X (resp. Y ) extraneous if it




• occurs as an X- (resp. Y -) coordinate of an element of M ; and
• does not occur as an X- (resp. Y -) coordinate of an element of M involved in a collision.
We claim that our bound


1−

Pr[ T = M | F = 0 ]
Pr[ T = M | F = 1 ]

2

= O(C 2 /n).

(39)

holds even after conditioning on which extraneous X are in A and which extraneous Y are in B. This will
be sufficient since which X and Y are extraneous, and which of them are in A or B, are both independent
of F . Let ME be the X and Y values of M along with the information on which extraneous X and Y are
in A or B. It is enough to bound


Pr[ T = M | F = 0 ]
1−
Pr[ T = M | F = 1 ]

2



=

Pr[ T = M | F = 1 ] − Pr[ T = M | F = 0 ]
Pr[ T = M | F = 1 ]

2

.

It thus suffices to bound
(dTV ((T | F = 0, ME ), (T | F = 1, ME )))2 .
Say that a pair of colliding elements of M is light if neither the corresponding X nor the corresponding Y
are in A or B. Observe that if we condition on no collision in M being light, the conditional distributions
of T on F = 0 and on F = 1 are the same. Therefore, the expression above is bounded by
Pr[ There exists a light collision in M ]2 .
Thus, it suffices to bound the probability that M contains a light collision given its values of X and Y . If M
(ignoring the values of W and Z) is ((x1 , y1 ), . . . , (xt , yt )) the probability of seeing this M conditioned on
the values of A and B is Θ(n−2t+|{ i∈[t] : xi ∈A }|/4+|{ i∈[t] : yi ∈B }|/4 ). Now given some set of a of values
of X appearing in M , the prior probability that those are the values of X in M appearing in A is n−a/4 · φ,
where φ is some quantity that changes by only a 1 + o(1) factor if a single element is added or removed from
the set of X’s. A similar relation holds for Y ’s. Given this, and conditioning on whether the extraneous X’s
and Y ’s are in A and B, we note that each non-extraneous X or Y that is in A or B contributes a factor of
roughly n−1/4 to the prior probability of having that configuration of elements in A or B, but contributes a
factor of at least n1/2 to the conditional probability of seeing the M that we saw given those values of A and
B. Therefore, even conditioned on the previously determined lightness of other collisions, each collision
has only a O(n−1/2 ) probability of being light. Therefore the probability that there is a light collision is
O(C/n1/2 ). This implies Eq. (39) and completes the proof of our lower bound.

48

References
[ADK15]

J. Acharya, C. Daskalakis, and G. Kamath. Optimal testing for properties of distributions.
In Proceedings of NIPS’15, 2015. 1.1, 1.2, 5.3

[Agr92]

A. Agresti. A survey of exact inference for contingency tables. Statist. Sci., 7(1):131–153,
02 1992. 1.1

[BFF+ 01]

T. Batu, E. Fischer, L. Fortnow, R. Kumar, R. Rubinfeld, and P. White. Testing random
variables for independence and identity. In Proc. 42nd IEEE Symposium on Foundations of
Computer Science, pages 442–451, 2001. 1.2, 5.3

[BFR+ 00]

T. Batu, L. Fortnow, R. Rubinfeld, W. D. Smith, and P. White. Testing that distributions are
close. In IEEE Symposium on Foundations of Computer Science, pages 259–269, 2000. 1.1

[BH07]

R. Blundell and J. L. Horowitz. A non-parametric test of exogeneity. The Review of Economic Studies, 74(4):1035–1058, 2007. 1.1

[BKR04]

T. Batu, R. Kumar, and R. Rubinfeld. Sublinear algorithms for testing monotone and unimodal distributions. In ACM Symposium on Theory of Computing, pages 381–390, 2004.
1.1

[BT14]

T. Bouezmarni and A. Taamouti. Nonparametric tests for conditional independence using
conditional distributions. Journal of Nonparametric Statistics, 26(4):697–719, 2014. 1.1

[Can15]

C. L. Canonne. A survey on distribution testing: Your data is big. but is it blue? Electronic
Colloquium on Computational Complexity (ECCC), 22:63, 2015. 1.1

[CDGR16]

C. L. Canonne, I. Diakonikolas, T. Gouleakis, and R. Rubinfeld. Testing shape restrictions
of discrete distributions. In 33rd Symposium on Theoretical Aspects of Computer Science,
STACS 2016, pages 25:1–25:14, 2016. See also [CDGR17] (full version). 1.1

[CDGR17]

C. L. Canonne, I. Diakonikolas, T. Gouleakis, and R. Rubinfeld. Testing shape restrictions
of discrete distributions. Theory of Computing Systems, pages 1–59, 2017. 7

[CDKS17]

C. L. Canonne, I. Diakonikolas, D. M. Kane, and A. Stewart. Testing Bayesian networks.
In Proceedings of the 30th Conference on Learning Theory, COLT 2017, pages 370–448,
2017. 1.1, 6.1, 6.1

[CDS17]

C. L. Canonne, I. Diakonikolas, and A. Stewart. Fourier-based testing for families of distributions. CoRR, abs/1706.05738, 2017. 1.1

[CDVV14]

S. Chan, I. Diakonikolas, P. Valiant, and G. Valiant. Optimal algorithms for testing closeness
of discrete distributions. In SODA, pages 1193–1203, 2014. (document), 1.1, 1.2

[Coc54]

W. G. Cochran. Some methods for strengthening the common χ2 tests. Biometrics,
10(4):417–451, 1954. 1.1

[Daw79]

A. P. Dawid. Conditional independence in statistical theory. Journal of the Royal Statistical
Society. Series B (Methodological), 41(1):1–31, 1979. 1.1

49

[DDK18]

C. Daskalakis, N. Dikkala, and G. Kamath. Testing Ising models. In SODA, 2018. To
appear. 1.1

[DDS+ 13]

C. Daskalakis, I. Diakonikolas, R. Servedio, G. Valiant, and P. Valiant. Testing k-modal
distributions: Optimal algorithms via reductions. In SODA, pages 1833–1852, 2013. 1.1

[DGPP16]

I. Diakonikolas, T. Gouleakis, J. Peebles, and E. Price. Collision-based testers are optimal for uniformity and closeness. Electronic Colloquium on Computational Complexity
(ECCC), 23:178, 2016. 1.1

[DGPP17]

I. Diakonikolas, T. Gouleakis, J. Peebles, and E. Price. Sample-optimal identity testing with
high probability. CoRR, abs/1708.02728, 2017. 1.1

[DK16]

I. Diakonikolas and D. M. Kane. A new approach for testing properties of discrete distributions. In FOCS, pages 685–694, 2016. Full version available at abs/1601.05557. (document), 1.1, 1.2, 1.2, 1.2, 1.3.1, 1.3.2, 1.3.3, 2.2, 2.2, 2.3, 5.3, 6.1, 6.1, B.3

[DKN15a]

I. Diakonikolas, D. M. Kane, and V. Nikishkin. Optimal algorithms and lower bounds for
testing closeness of structured distributions. In 56th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2015, 2015. 1.1

[DKN15b]

I. Diakonikolas, D. M. Kane, and V. Nikishkin. Testing identity of structured distributions.
In Proceedings of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms,
SODA 2015, San Diego, CA, USA, January 4-6, 2015, 2015. 1.1

[DKN17]

I. Diakonikolas, D. M. Kane, and V. Nikishkin. Near-optimal closeness testing of discrete
histogram distributions. In 44th International Colloquium on Automata, Languages, and
Programming, ICALP 2017, pages 8:1–8:15, 2017. 1.1

[DM01]

M. A. Delgado and W. G. Manteiga. Significance testing in nonparametric regression based
on the bootstrap. The Annals of Statistics, 29(5):1469–1507, 2001. 1.1

[dMASdBP14] P. de Morais Andrade, J. M. Stern, and C. A. de Braganca Pereira. Bayesian test of significance for conditional independence: The multinomial model. Entropy, 16(3):1376–1395,
2014. 1.1
[Dob59]

R. L. Dobrušin. A general formulation of the fundamental theorem of Shannon in the theory
of information. Uspehi Mat. Nauk, 14(6 (90)):3–104, 1959. 2.1

[DP17]

C. Daskalakis and Q. Pan. Square Hellinger subadditivity for Bayesian networks and its
applications to identity testing. In Proceedings of the 30th Conference on Learning Theory,
COLT 2017, pages 697–703, 2017. 1.1

[EO87]

D. Easley and M. O’Hara. Price, trade size, and information in securities markets. Journal
of Financial Economics, 19(1):69 – 90, 1987. 1.1

[Fis24]

R. A. Fisher. The distribution of the partial correlation coefficient. Metron, 3:329–332,
1924. 1.1

[Gol17]

O. Goldreich. Introduction to Property Testing. Cambridge University Press, 2017. 1.1

50

[Gra80]

C.W.J. Granger. Testing for causality: A personal viewpoint. Journal of Economic Dynamics and Control, 2(Supplement C):329 – 352, 1980. 1.1

[GS10]

G. Geenens and L. Simar. Nonparametric tests for conditional independence in two-way
contingency tables. Journal of Multivariate Analysis, 101(4):765–788, 2010. 1.1

[HPS16]

M. Hardt, E. Price, and N. Srebro. Equality of opportunity in supervised learning. In
Advances in Neural Information Processing Systems 29: Annual Conference on Neural
Information Processing Systems 2016, pages 3315–3323, 2016. 1.1

[Hua10]

T.-M. Huang. Testing conditional independence using maximal nonlinear conditional correlation. Ann. Statist., 38(4):2047–2091, 08 2010. 1.1

[LG96]

O. Linton and P. Gozalo. Conditional Independence Restrictions: Testing and Estimation. Cowles Foundation Discussion Papers 1140, Cowles Foundation for Research in Economics, Yale University, 1996. 1.1

[LRR11]

R. Levi, D. Ron, and R. Rubinfeld. Testing properties of collections of distributions. In
ICS, pages 179–194, 2011. 1.2, 5.3

[MH59]

N. Mantel and W. Haenszel. Statistical aspects of the analysis of data from retrospective
studies of disease. Journal of the National Cancer Institute, 22(4):719–748, April 1959.
PMID: 13655060. 1.1, 1.1

[Nea03]

R. E. Neapolitan. Learning Bayesian Networks. Prentice-Hall, Inc., 2003. 1.1

[NUU17]

K. Natori, M. Uto, and M. Ueno. Consistent learning Bayesian networks with thousands of
variables. In Proceedings of The 3rd International Workshop on Advanced Methodologies
for Bayesian Networks, volume 73 of Proceedings of Machine Learning Research, pages
57–68. PMLR, 20–22 Sep 2017. 1.1

[Pan08]

L. Paninski. A coincidence-based test for uniformity given very sparsely-sampled discrete
data. IEEE Transactions on Information Theory, 54:4750–4755, 2008. 1.1, 6.3

[Pea88]

J. Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference.
Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. 1.1

[Pin05]

M. S. Pinsker. On the estimation of information via variation. Problemy Peredachi Informatsii, 41(2):3–8, 2005. A

[Rub12]

R. Rubinfeld. Taming big probability distributions. XRDS, 19(1):24–28, 2012. 1.1

[SGS00]

P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT press,
2nd edition, 2000. 1.1, 1.1

[Son09]

K. Song. Testing conditional independence via Rosenblatt transforms.
37(6B):4011–4045, 12 2009. 1.1

[SW07]

L. Su and H. White. A consistent characteristic function-based test for conditional independence. Journal of Econometrics, 141(2):807 – 834, 2007. 1.1

51

Ann. Statist.,

[SW08]

L. Su and H. White. A nonparametric Hellinger metric test for conditional independence.
Econometric Theory, 24(4):829–864, 2008. 1.1

[SW14]

L. Su and H. White. Testing conditional independence via empirical likelihood. Journal
of Econometrics, 182(1):27 – 44, 2014. Causality, Prediction, and Specification Analysis:
Recent Advances and Future Directions. 1.1

[TBA06]

I. Tsamardinos, L. E. Brown, and C. F. Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine Learning, 65(1):31–78, Oct 2006. 1.1

[VV14]

G. Valiant and P. Valiant. An automatic inequality prover and instance optimal identity
testing. In FOCS, 2014. 1.1

[WH17]

X. Wang and Y. Hong. Characteristic function based testing for conditional independence:
a nonparametric regression approach. Econometric Theory, pages 1–35, 2017. 1.1, 1.1

[Wyn78]

A. D. Wyner. A definition of conditional mutual information for arbitrary ensembles. Inform. and Control, 38(1):51–59, 1978. 2.1

[ZPJS11]

K. Zhang, J. Peters, D. Janzing, and B. Schölkopf. Kernel-based conditional independence
test and application in causal discovery. In Proceedings of the Twenty-Seventh Conference
on Uncertainty in Artificial Intelligence, UAI’11, pages 804–813. AUAI Press, 2011. 1.1

A

Testing with Respect to Mutual Information

We conclude by considering a slightly different model from the one considered thus far. In particular, while
the total variation metric is a reasonable one to measure what it means for X and Y to be far from conditionally independent, there is another metric that is natural in this context: conditional mutual information.
Specifically, we modify the testing problem to distinguish between the cases where X and Y are conditionally independent on Z and the case where I(X; Y |Z) ≥ ε. Our picture here is somewhat less complete, but
we are still able to say something in the case where X, Y are binary.
Theorem A.1. If X and Y are binary random variables and Z has a support of size n, there exists a
sample-efficient algorithm that distinguishes between I(X; Y |Z) = 0 and I(X; Y |Z) ≥ ε with sample
complexity
√
O(max(min(n6/7 log8/7 (1/ε)/ε8/7 , n7/8 log(1/ε)/ε), n log2 (1/ε)/ε2 )) .
Proof. This follows immediately upon noting that by Lemma A.1 (stated and proven later), that if X and
Y are ε-close in total variation distance from being conditionally independent on Z, then I(X; Y |Z) ≤
O(ε log(1/ε)); or, by the contrapositive, that I(X; Y |Z) ≥ ε implies that X and Y are Ω(ε/ log(1/ε))far in total variation distance from being conditionally independent on Z. Therefore, it suffices to run our
existing conditional independence tester with parameter ε0 := Ω(ε/ log(1/ε)). The sample complexity of
this tester is as specified.
Remark A.2 (On the optimality of this bound). It is not difficult to modify the analysis slightly in order to
remove the logarithmic factors from the first two terms in the above expression. Intuitively, this is because
these terms arise only when at least half of the mutual information comes from “light” bins, with mass at
52

most 1/m. In this case, these bins contribute at least m4 z ε2z pZ (z)4  m4 z (pZ (z)εz log(1/εz ))4 
m4 ε4 /n3 to the expectation of Z, and the analysis proceeds from there as before.
It is also easy to show that in this regime our lower bounds still apply, as the hard instances also produced
distributions with mutual information Ω(ε).1 Therefore, we have matching upper and lower bounds as long
as ε  n−3/8 / log2 n.
However, it seems likely that the correct behavior in the small ε regime is substantially different when
testing with respect to mutual information. The difficult cases for total variation distance testing actually
end up with mutual information merely I(X; Y |Z) = O(ε2 ). It is quite possible that a better algorithm or a
better analysis of the existing algorithm could give substantially improved performance when ε < n−3/8 . In
fact, it is conceivable that the sample complexity of O(n7/8 /ε) could be maintained for a broad range of ε.
The only lower bound that we know preventing this is a lower bound of Ω(ε log(1/ε)) by noting that there
are distributions with I(X; Y |Z) ≥ ε, but where (X, Y, Z) is O(ε/ log(1/ε))-far in variation distance from
being conditionally independent.
P

P

Lemma A.1. Assume (X, Y, Z) ∼ p, where p ∈ ∆(X × Y × Z) with |X | = `1 , |Y| = `2 , and |Z| = n.
Then, for every ε ∈ (0, 1),








• If dTV p, PX ,Y|Z ≤ ε, then I(X; Y |Z) ≤ O(ε log(`1 `2 /ε));
• If dTV p, PX ,Y|Z ≥ ε, then I(X; Y |Z) ≥ 2ε2 .
Proof. The second item is simply an application of Pinsker’s inequality, recalling that
I(X; Y |Z) = dKL ((X, Y ) | Z || (X | Z) ⊗ (Y | Z)) .
i.e. the Kullback–Leibler divergence between the joint distribution of (X, Y | Z) and the product of
marginals (X | Z) and (Y | Z). As for the first, it follows from the relation between conditional mutual information and total variation distance obtained in [Pin05] (and Lemma 2.2).
1

I.e., the conditional mutual information of “no-distributions” is easily seen to actually be Ω(ε), while applying the relation
between total variation distance and conditional mutual information as a black-box to the ε distance in total variation distance would
incur a quadratic loss in ε.

53

B
B.1

Deferred Proofs from Section 2
Proof of Lemma 2.1

The proof follows from the following chain of (in-)equalities:
2dTV p, p0 =

p(i, j, z) − p0 (i, j, z)

X



(i,j,z)∈X ×Y×Z

pZ (z) · pz (i, j) − p0Z (z) · p0z (i, j)

X

=

(i,j,z)∈X ×Y×Z

pZ (z) · (pz (i, j) − p0z (i, j)) + (pZ (z) − p0Z (z)) · p0z (i, j)

X

=

(i,j,z)∈X ×Y×Z

pZ (z) · pz (i, j) − p0z (i, j) +

X

≤

(i,j,z)∈X ×Y×Z

(i,j,z)∈X ×Y×Z



=

X



=2



pz (i, j) − p0z (i, j)  +

X

pZ (z) ·

z∈Z

pZ (z) ·

X



 pZ (z) − p0 (z) ·
Z

z∈Z

(i,j)∈X ×Y

X

pZ (z) − p0Z (z) · p0z (i, j)

X


dTV pz , p0z

+


2dTV pZ , p0Z

p0z (i, j)

X
(i,j)∈X ×Y

,

z∈Z

where the fourth line used the triangle inequality and the last line used the fact that (i,j)∈X ×Y p0z (i, j) = 1.
This completes the proof of the first part of the lemma. For the second part, we note that the equality in (2)
holds if and only if the triangle inequality in the fourth line above holds with equality, i.e., when pZ = p0Z .
This completes the proof of Lemma 2.1.
P

B.2

Proof of Lemma 2.2

Let p0 ∈ PX ,Y|Z be such that dTV (p, p0 ) ≤ ε and q = z∈Z pZ (z)qz . Since dTV (p, q) ≤ dTV (p, p0 ) +
dTV (p0 , q) ≤ ε + dTV (p0 , q), it suffices to show that dTV (p0 , q) ≤ 3ε. By Lemma 2.1, we have that
P

dTV q, p0 ≤


X

qZ (z) · dTV qz , p0z + dTV qZ , p0Z




z∈Z

=

X

pZ (z) · dTV qz , p0z + dTV pZ , p0Z




z∈Z

=

X





pZ (z) · dTV pz,X ⊗ pz,Y , p0z,X ⊗ p0z,Y + dTV pZ , p0Z



z∈Z

≤

X









pZ (z) · dTV pz,X , p0z,X + dTV pz,Y , p0z,Y



+ dTV pZ , p0Z



z∈Z

=

X
z∈Z





pZ (z)dTV pz,X , p0z,X +

X





pZ (z)dTV pz,Y , p0z,Y + dTV pZ , p0Z



z∈Z

≤ 3ε ,
where the second line uses the fact that qZ = pZ , the third line uses the fact that qz = pz,X ⊗ pz,Y
(Definition 2.1) and that p0z = p0z,X ⊗ p0z,Y (since p0 ∈ PX ,Y|Z ), the fourth line uses the sub-additivity of
total variation distance for product distributions, and the last line uses the fact that each of the three terms in
the fifth line is bounded from above by dTV (p, p0 ). This completes the proof of Lemma 2.2.
54

B.3

Proof of Lemma 2.3

This lemma is essentially shown in [DK16]. The only difference is that we require a proof for (ii) when S
is a set of m independent samples (as opposed to Poi(m) samples from p in [DK16]). We show this by an
explicit calculation below.
P
Let ai equal one plus the number of copies of i in S, i.e. ai := 1 + j∈S 1{i=j} . We note that the
hP

i

i
expected squared `2 -norm of pS is E ni=1 aj=1
p2i /a2i =
1 + X where X is a Bin(m, pi ) random variable. Therefore,

P

Pn

2
i=1 pi E[1/ai ] .

Further, ai is distributed as


m
m
X
X
1
1
m k
m + 1 k+1
1
=
E
pi (1 − pi )m−k =
pi (1 − pi )(m+1)−(k+1)
1+X
k
+
1
(m
+
1)p
k
k
+
1
i k=0
k=0
!



!

m+1
X m+1
1
1
1 − (1 − pi )m+1
≤
.
=
p`i (1 − pi )(m+1)−` =
(m + 1)pi `=1
(m + 1)pi
(m + 1)pi
`

!

h

i

This implies E [kpS k22 ≤

B.4

Pn

2
i=1 pi /(mpi )

= (1/m)

Pn

i=1 pi

= 1/m. which completes the proof.

Proof of Claim 2.1

Recalling that E[N ] = λ and E N 2 = λ + λ2 , we get


h

i



E N 1{N ≥4} = e

∞
X
λk

−λ

k

k=4

k!

=λ−e

−λ



1
λ + λ + λ3 := f (λ) ,
2


2

and


2

Var N 1{N ≥4} = λ + λ − e

−λ

3
λ + 2λ + λ3
2



2





−λ

− λ−e



1
λ + λ + λ3
2
2

2

:= g(λ) .

(x)
(x)
From these expressions, it is easy to check that (i) limx→0 fg(x)
= 14 , and (ii) limx→∞ fg(x)
= 1. From the
definition as a variance of a non-constant random variable, it follows that g(x) > 0 for all x > 0, from which
we get that (iii) fg is continuous and positive on [0, ∞). Combining these three statements, we get that fg
achieves a minimum c on [0, ∞), and that this minimum is positive. This implies the result with C := 1/c.
The value 4.22 comes from studying numerically this ratio, whose minimum is achieved for x ' 1.1457.

B.5

Proof of Claim 2.2

Let a, b ≥ 0, λ > 0, and assume X ∼ Poisson(λ).
Without loss of generality, suppose 0 < a ≤ b (the case
p
a = 0 being trivial). We can rewrite X min(X, a) min(X, b)1{X≥4} as Y with
√
√
Y := X 2 1{X≤a} + aX 3/2 1{a<X≤b} + abX1{X>b}
which implies
Y 2 = X 4 1{X≤a} + aX 3 1{a<X≤b} + abX 2 1{X>b} .
By linearity of expectation, the original claim boils down to proving there exists C > 0 such that
h

i

h

i

h

i

E X 4 1{4≤X≤a} + E aX 3 1{a<X≤b} 1{X≥4} + E abX 2 1{X>b} 1{X≥4}
 h
i
h√
i
h√
i
≤ C E X 2 1{4≤X≤a} + E aX 3/2 1{a<X≤b} 1{X≥4} + E abX1{X>b} 1{X≥4}
 h
i
h√
i
h√
i2
+ E X 2 1{4≤X≤a} + E aX 3/2 1{a<X≤b} 1{X≥4} + E abX1{X>b} 1{X≥4}
55

and since (x + y + z)2 ≥ x2 + y 2 + z 2 for x, y, z ≥ 0, it is enough to show
h

i

h

i

h

E β 2 X 2α 1{X∈S} ≤ CβE X α 1{X∈S} + β 2 E X α 1{X∈S}

i2

for α, β > 0, and S ⊆ R+ an interval. This in turn follows from arguments similar to that of the proof
of Claim 2.1.

B.6

Proof of Claim 2.3
h p

i

For λ < 8, we can take bound the expectation by the contribution of X = 4 as E X min(X, a) min(X, b)1{X≥4} ≥
4 min(4, a) min(4, p
b)λ4 /4! ≥ λ4 /3. Then λ min(λ, a) min(λ, b) ≥ λ min(λ, 2) min(λ, 2) ≥ λ2 /4 ≥
4
λ /256. Thus min(λ min(λ, a) min(λ, b), λ4 ) ≤ 256λ4 . Putting this together we have that for λ < 8,
p

p

 q

p



q

E X min(X, a) min(X, b)1{X≥4} ≥ (1/768) min(λ min(λ, a) min(λ, b), λ4 ) .
To deal with the λ ≥ 8 case, we claim that Pr[ X ≥ bλ/2c ] ≥ 1/2 in this case. To see this, we just need to
P
k
expand 1 = exp(−λ) ∞
k=0 f (k)λ /k! and note that for 1 ≤ k ≤ λ/2, the ratio of the k term to the k − 1
term is at least λ/k ≥ 2. Thus the sum of the first bλ/2c terms is smaller than the k = bλ/2c term and so
bλ/2c−1

exp(−λ)

X

k

λ /k! ≤ exp(−λ)

k=0

∞
X

λk /k .

k=bλ/2c

The RHS is Pr[ X ≥ bλ/2c ] and the LHS is Pr[ X < bλ/2c ] = 1 − Pr[ X ≥ bλ/2c ] and so we have that
Pr[ X ≥ bλ/2c ] ≥ 1/2 as claimed.
For 8 ≤ λ ≤ 2 min a, b, we have
 q



h

E X min(X, a) min(X, b)1{X≥4} ≥ E X 2 1{X≥bλ/2c}

i

≥ (1/2)(bλ/2c)2 ≥ λ2 /3
q

≥ min(λ min(λ, a) min(λ, b)/6, λ4 /3)
q

≥ (1/6) min(λ min(λ, a) min(λ, b), λ4 )
For λ ≥ 2 max a, b, 4, noting that for X ≥ λ/2,
 q

p

min(X, a) min(X, b) =



√

ab, we have

 q



E X min(X, a) min(X, b)1{X≥4} ≥ E X min(X, a) min(X, b)1{X≥bλ/2c}
√
≥ bλ/2c ab
√
≥ abλ/3
q

≥ (1/6) min(λ min(λ, a) min(λ, b), λ4 )
The final case we need to consider is when λ is between 2a and 2b and the maximum of those is over 4

56

Supposing without loss of generality that a ≤ b, for max 2a, 4 ≤ λ ≤ 2b, we have
 q



 q



E X min(X, a) min(X, b)1{X≥4} ≥ E X min(X, a) min(X, b)1{X≥bλ/2c}
√
≥ bλ/2c3/2 a
√
≥ abλ3/2 /4
q

≥ (1/8) min(λ min(λ, a) min(λ, b), λ4 ) .

57

