Model-Powered Conditional Independence Test
Rajat Sen1,* , Ananda Theertha Suresh2,* , Karthikeyan Shanmugam3,* , Alexandros G.
Dimakis1 , and Sanjay Shakkottai1

arXiv:1709.06138v1 [stat.ML] 18 Sep 2017

1

The University of Texas at Austin
2
Google, New York
3
IBM Research, Thomas J. Watson Center
September 20, 2017
Abstract
We consider the problem of non-parametric Conditional Independence testing (CI testing) for continuous random variables. Given i.i.d samples from the joint distribution f (x, y, z) of continuous random
vectors X, Y and Z, we determine whether X ⊥
⊥ Y |Z. We approach this by converting the conditional
independence test into a classification problem. This allows us to harness very powerful classifiers
like gradient-boosted trees and deep neural networks. These models can handle complex probability
distributions and allow us to perform significantly better compared to the prior state of the art, for
high-dimensional CI testing. The main technical challenge in the classification problem is the need for
samples from the conditional product distribution f CI (x, y, z) = f (x|z)f (y|z)f (z) – the joint distribution
if and only if X ⊥
⊥ Y |Z. – when given access only to i.i.d. samples from the true joint distribution f (x, y, z).
To tackle this problem we propose a novel nearest neighbor bootstrap procedure and theoretically show
that our generated samples are indeed close to f CI in terms of total variational distance. We then develop
theoretical results regarding the generalization bounds for classification for our problem, which translate
into error bounds for CI testing. We provide a novel analysis of Rademacher type classification bounds
in the presence of non-i.i.d near-independent samples. We empirically validate the performance of our
algorithm on simulated and real datasets and show performance gains over previous methods.

1

Introduction

Testing datasets for Conditional Independence (CI) have significant applications in several statistical/learning
problems; among others, examples include discovering/testing for edges in Bayesian networks [15, 27, 7, 9],
causal inference [23, 14, 29, 5] and feature selection through Markov Blankets [16, 31]. Given a triplet of
random variables/vectors (X, Y, Z), we say that X is conditionally independent of Y given Z (denoted by
X ⊥
⊥ Y |Z), if the joint distribution fX,Y,Z (x, y, z) factorizes as fX,Y,Z (x, y, z) = fX|Z (x|z)fY |Z (y|z)fZ (z).
The problem of Conditional Independence Testing (CI Testing) can be defined as follows: Given n i.i.d samples
from fX,Y,Z (x, y, z), distinguish between the two hypothesis H0 : X ⊥
⊥ Y |Z and H1 : X ⊥
6 ⊥ Y |Z.
In this paper we propose a data-driven Model-Powered CI test. The central idea in a model-driven approach
is to convert a statistical testing or estimation problem into a pipeline that utilizes the power of supervised
learning models like classifiers and regressors; such pipelines can then leverage recent advances in classification/regression in high-dimensional settings. In this paper, we take such a model-powered approach
(illustrated in Fig. 1), which reduces the problem of CI testing to Binary Classification. Specifically, the key
steps of our procedure are as follows:
(i) Suppose we are provided 3n i.i.d samples from fX,Y,Z (x, y, z). We keep aside n of these original samples
in a set U1 (refer to Fig. 1). The remaining 2n of the original samples are processed through our first module,
the nearest-neighbor bootstrap (Algorithm 1 in our paper), which produces n simulated samples stored in
* Equal Contribution

1

3n Original
Samples

XY Z

x 1 y 1 z1

U1

n
Original Samples

.
.
.

..
.
1

Training
Set
Dr
G

ĝ (Trained Classifier)

y3n z3n

Shu✏e

U20

XY Z

`

XY Z

.
.
.

2n Original
Samples

1
+

..
.

x3n

`

XY Z

.
.
.

Nearest
Neighbor
Bootstrap

0

+

..
.

0
n samples
close to f CI

Test
Set De

L(ĝ, De ) (Test Error)

Figure 1: Illustration of our methodology. A part of the original samples are kept aside in U1 . The rest of
the samples are used in our nearest neighbor boot-strap to generate a data-set U20 which is close to f CI in
distribution. The samples are labeled as shown and a classifier is trained on a training set. The test error is
measured on a test set there-after. If the test-error is close to 0.5, then H0 is not rejected, however if the test
error is low then H0 is rejected.
U20 . In Section 3, we show that these generated samples in U20 are in fact close in total variational distance
(defined in Section 3) to the conditionally independent distribution f CI (x, y, z) , fX|Z (x|z)fY |Z (y|z)fZ (z).
(Note that only under H0 does the equality f CI (.) = fX,Y,Z (.) hold; our method generates samples close to
f CI (x, y, z) under both hypotheses).
(ii) Subsequently, the original samples kept aside in U1 are labeled 1 while the new samples simulated from
the nearest-neighbor bootstrap (in U20 ) are labeled 0. The labeled samples (U1 with label 1 and U20 labeled
0) are aggregated into a data-set D. This set D is then broken into training and test sets Dr and De each
containing n samples each.
(iii) Given the labeled training data-set (from step (ii)), we train powerful classifiers such as gradient boosted
trees [6] or deep neural networks [17] which attempt to learn the classes of the samples. If the trained
classifier has good accuracy over the test set, then intuitively it means that the joint distribution fX,Y,Z (.)
is distinguishable from f CI (note that the generated samples labeled 0 are close in distribution to f CI ).
Therefore, we reject H0 . On the other hand, if the classifier has accuracy close to random guessing, then
fX,Y,Z (.) is in fact close to f CI , and we fail to reject H0 .
For independence testing (i.e whether X ⊥
⊥ Y ), classifiers were recently used in [19]. Their key observation
was that given i.i.d samples (X, Y ) from fX,Y (x, y), if the Y coordinates are randomly permuted then the
resulting samples exactly emulate the distribution fX (x)fY (y). Thus the problem can be converted to a
two sample test between a subset of the original samples and the other subset which is permuted - Binary
classifiers were then harnessed for this two-sample testing; for details see [19]. However, in the case of CI
testing we need to emulate samples from f CI . This is harder because the permutation of the samples needs
to be Z dependent (which can be high-dimensional). One of our key technical contributions is in proving
that our nearest-neighbor bootstrap in step (i) achieves this task.
The advantage of this modular approach is that we can harness the power of classifiers (in step (iii) above),
which have good accuracies in high-dimensions. Thus, any improvements in the field of binary classification
imply an advancement in our CI test. Moreover, there is added flexibility in choosing the best classifier based
on domain knowledge about the data-generation process. Finally, our bootstrap is also efficient owing to fast
algorithms for identifying nearest-neighbors [24].

1.1

Main Contributions

(i) (Classification based CI testing) We reduce the problem of CI testing to Binary Classification as
detailed in steps (i)-(iii) above and in Fig. 1. We simulate samples that are close to f CI through a novel
nearest-neighbor bootstrap (Algorithm 1) given access to i.i.d samples from the joint distribution. The
2

problem of CI testing then reduces to a two-sample test between the original samples in U1 and U20 , which
can be effectively done by binary classifiers.
(ii) (Guarantees on Bootstrapped Samples) As mentioned in steps (i)-(iii), if the samples generated by
the bootstrap (in U20 ) are close to f CI , then the CI testing problem reduces to testing whether the data-sets U1
and U20 are distinguishable from each other. We theoretically justify that this is indeed true. Let φX,Y,Z (x, y, z)
denote the distribution of a sample produced by Algorithm 1, when it is supplied with 2n i.i.d samples
from fX,Y,Z (.). In Theorem 1, we prove that dT V (φ, f CI ) = O(1/n1/dz ) under appropriate smoothness
assumptions. Here dz is the dimension of Z and dT V denotes total variational distance (Def. 1).
(iii) (Generalization Bounds for Classification under near-independence) The samples generated
from the nearest-neighbor bootstrap do not remain i.i.d but they are close to i.i.d. We quantify this property
and go on to show generalization risk bounds for the classifier. Let us denote the class of function encoded
by the classifier as G. Let R̂ denote the probability of error of the optimal classifier ĝ ∈ G trained on the
training set (Fig. 1). We prove that under appropriate assumptions, we have



q
√
r0 − O(1/n1/dz ) ≤ R̂ ≤ r0 + O(1/n1/dz ) + O
V n−1/3 + 2dz /n
with high probability, upto log factors. Here r0 = 0.5(1 − dT V (f, f CI )), V is the VC dimension [30] of the
class G. Thus when f is equivalent to f CI (H0 holds) then the error rate of the classifier is close to 0.5. But
when H1 holds the loss is much lower. We provide a novel analysis of Rademacher complexity bounds [4]
under near-independence which is of independent interest.
(iv) (Empirical Evaluation) We perform extensive numerical experiments where our algorithm outperforms
the state of the art [32, 28]. We also apply our algorithm for analyzing CI relations in the protein signaling
network data from the flow cytometry data-set [26]. In practice we observe that the performance with respect
to dimension of Z scales much better than expected from our worst case theoretical analysis. This is because
powerful binary classifiers perform well in high-dimensions.

1.2

Related Work

In this paper we address the problem of non-parametric CI testing when the underlying random variables are
continuous. The literature on non-parametric CI testing is vast. We will review some of the recent work in
this field that is most relevant to our paper.
Most of the recent work in CI testing are kernel based [28, 32, 10]. Many of these works build on the study
in [11], where non-parametric CI relations are characterized using covariance operators for Reproducing
Kernel Hilbert Spaces (RKHS) [11]. KCIT [32] uses the partial association of regression functions relating X,
Y , and Z. RCIT [28] is an approximate version of KCIT that attempts to improve running times when the
number of samples are large. KCIPT [10] is perhaps most relevant to our work. In [10], a specific permutation
of the samples is used to simulate data from f CI . An expensive linear program needs to be solved in order to
calculate the permutation. On the other hand, we use a simple nearest-neighbor bootstrap and further we
provide theoretical guarantees about the closeness of the samples to f CI in terms of total variational distance.
Finally the two-sample test in [10] is based on a kernel method [3], while we use binary classifiers for the same
purpose. There has also been recent work on entropy estimation [13] using nearest neighbor techniques (used
for density estimation); this can subsequently be used for CI testing by estimating the conditional mutual
information I(X; Y |Z).
Binary classification has been recently used for two-sample testing, in particular for independence testing [19].
Our analysis of generalization guarantees of classification are aimed at recovering guarantees similar to [4],
but in a non-i.i.d setting. In this regard (non-i.i.d generalization guarantees), there has been recent work in
proving Rademacher complexity bounds for β-mixing stationary processes [21]. This work also falls in the
category of machine learning reductions, where the general philosophy is to reduce various machine learning
settings like multi-class regression [2], ranking [1], reinforcement learning [18], structured prediction [8] to
that of binary classification.
3

2

Problem Setting and Algorithms

In this section we describe the algorithmic details of our CI testing procedure. We first formally define our
problem. Then we describe our bootstrap algorithm for generating the data-set that mimics samples from
f CI . We give a detailed pseudo-code for our CI testing process which reduces the problem to that of binary
classification. Finally, we suggest further improvements to our algorithm.
Problem Setting: The problem setting is that of non-parametric Conditional Independence (CI) testing
given i.i.d samples from the joint distributions of random variables/vectors [32, 10, 28]. We are given 3n i.i.d
samples from a continuous joint distribution fX,Y,Z (x, y, z) where x ∈ Rdx , y ∈ Rdy and z ∈ Rdz . The goal is to
test whether X ⊥
⊥ Y |Z i.e whether fX,Y,Z (x, y, z) factorizes as, fX,Y,Z (x, y, z) = fX|Z (x|z)fY |Z (y|z)fZ (z) ,
f CI (x, y, z)
This is essentially a hypothesis testing problem where: H0 : X ⊥
⊥ Y |Z and H1 : X ⊥
6 ⊥ Y |Z.
Note: For notational convenience, we will drop the subscripts when the context is evident. For instance we
may use f (x|z) in place of fX|Z (x|z).
Nearest-Neighbor Bootstrap: Algorithm 1 is a procedure to generate a data-set U 0 consisting of n
samples given a data-set U of 2n i.i.d samples from the distribution fX,Y,Z (x, y, z). The data-set U is broken
into two equally sized partitions U1 and U2 . Then for each sample in U1 , we find the nearest neighbor in U2 in
terms of the Z coordinates. The Y -coordinates of the sample from U1 are exchanged with the Y -coordinates
of its nearest neighbor (in U2 ); the modified sample is added to U 0 .
Algorithm 1 DataGen - Given data-set U = U1 ∪ U2 of 2n i.i.d samples from f (x, y, z) (|U1 | = |U2 | = n ),
returns a new data-set U 0 having n samples.
function DataGen(U1 , U2 , 2n)
U0 = ∅
for u in U1 do
Let v = (x0 , y 0 , z 0 ) ∈ U2 be the sample such that z 0 is the 1-Nearest Neighbor (1-NN)
of z (in `2 norm) in the whole data-set U2 , where u = (x, y, z)
5:
Let u0 = (x, y 0 , z) and U 0 = U 0 ∪ {u0 }.
6:
end for
7: end function
1:
2:
3:
4:

One of our main results is that the samples in U 0 , generated in Algorithm 1 mimic samples coming from the
distribution f CI . Suppose u = (x, y, z) ∈ U1 be a sample such that fZ (z) is not too small. In this case z 0
(the 1-NN sample from U2 ) will not be far from z. Therefore given a fixed z, under appropriate smoothness
assumptions, y 0 will be close to an independent sample coming from fY |Z (y|z 0 ) ∼ fY |Z (y|z). On the other
hand if fZ (z) is small, then z is a rare occurrence and will not contribute adversely.
CI Testing Algorithm: Now we introduce our CI testing algorithm, which uses Algorithm 1 along with
binary classifiers. The psuedo-code is in Algorithm 2 (Classifier CI Test -CCIT).
Algorithm 2 CCITv1 - Given data-set U of 3n i.i.d samples from f (x, y, z), returns if X ⊥
⊥ Y |Z.

function CCIT(U, 3n, τ, G)
Partition U into three disjoint partitions U1 , U2 and U3 of size n each, randomly.
Let U20 = DataGen(U2 , U3 , 2n) (Algorithm 1). Note that |U20 | = n.
Create Labeled data-set D := {(u, ` = 1)}u∈U1 ∪ {(u0 , `0 = 0)}u0 ∈U20
Divide data-set D into train and test
Pset Dr and De respectively. Note that |Dr | = |De | = n.
Let ĝ = argming∈G L̂(g, Dr ) := |D1r | (u,`)∈Dr 1{g(u) 6= l}. This is Empirical Risk Minimization for
training the classifier (finding the best function in the class G).
7:
If L̂(ĝ, De ) > 0.5 − τ , then conclude X ⊥
⊥ Y |Z, otherwise, conclude X ⊥
6 ⊥ Y |Z.
8: end function
1:

2:
3:
4:
5:
6:

4

In Algorithm 2, the original samples in U1 and the nearest-neighbor bootstrapped samples in U20 should be
almost indistinguishable if H0 holds. However, if H1 holds, then the classifier trained in Line 6 should be
able to easily distinguish between the samples corresponding to different labels. In Line 6, G denotes the
space of functions over which risk minimization is performed in the classifier.
We will show (in Theorem 1) that the variational distance between the distribution of one of the samples
in U20 and f CI (x, y, z) is very small for large n. However, the samples in U20 are not exactly i.i.d but close
to i.i.d. Therefore, in practice for finite n, there is a small bias b > 0 i.e. L̂(ĝ, De ) ∼ 0.5 − b, even when H0
holds. The threshold τ needs to be greater than b in order for Algorithm 2 to function. In the next section,
we present an algorithm where this bias is corrected.
Algorithm with Bias Correction: We present an improved bias-corrected version of our algorithm as
Algorithm 3. As mentioned in the previous section, in Algorithm 2, the optimal classifier may be able to
achieve a loss slightly less that 0.5 in the case of finite n, even when H0 is true. However, the classifier is
expected to distinguish between the two data-sets only based on the Y, Z coordinates, as the joint distribution
of X and Z remains the same in the nearest-neighbor bootstrap. The key idea in Algorithm 3 is to train a
classifier only using the Y and Z coordinates, denoted by ĝ 0 . As before we also train another classier using
all the coordinates, which is denoted by ĝ. The test loss of ĝ 0 is expected to be roughly 0.5 − b, where b is the
bias mentioned in the previous section. Therefore, we can just subtract this bias. Thus, when H0 is true
L̂(ĝ 0 , De0 ) − L̂(ĝ, De ) will be close to 0. However, when H1 holds, then L̂(ĝ, De ) will be much lower, as the
classifier ĝ has been trained leveraging the information encoded in all the coordinates.
Algorithm 3 CCITv2 - Given data-set U of 3n i.i.d samples, returns whether X ⊥
⊥ Y |Z.

function CCIT(U, 3n, τ, G)
Perform Steps 1-5 as in Algorithm 2.
Let Dr0 = {((y, z), `)}(u=(x,y,z),`)∈Dr . Similarly, let De0 = {((y, z), `)}(u=(x,y,z),`)∈De . These are the
training and test sets without the X-coordinates.
P
4:
Let ĝ = argming∈G L̂(g, Dr ) := |D1r | (u,`)∈Dr 1{g(u) 6= l}. Compute test loss: L̂(ĝ, De ).
P
5:
Let ĝ 0 = argming∈G L̂(g, Dr0 ) := |D10 | (u,`)∈Dr0 1{g(u) 6= l}. Compute test loss: L̂(ĝ 0 , De0 ).
1:

2:
3:

r

6:
7:

3

6 ⊥ Y |Z, otherwise, conclude X ⊥
⊥ Y |Z.
If L̂(ĝ, De ) < L̂(ĝ 0 , De0 ) − τ , then conclude X ⊥
end function

Theoretical Results

In this section, we provide our main theoretical results. We first show that the distribution of any one of the
samples generated in Algorithm 1 closely resemble that of a sample coming from f CI . This result holds for a
broad class of distributions fX,Y,Z (x, y, z) which satisfy some smoothness assumptions. However, the samples
generated by Algorithm 1 (U2 in the algorithm) are not exactly i.i.d but close to i.i.d. We quantify this and
go on to show that empirical risk minimization over a class of classifier functions generalizes well using these
samples. Before, we formally state our results we provide some useful definitions.
Definition 1. The total variational distance between two continuous probability distributions f (.) and
g(.) defined over a domain X is, dT V (f, g) = supp∈B |Ef [p(X)] − Eg [p(X)]| where B is the set of all measurable
functions from X → [0, 1]. Here, Ef [.] denotes expectation under distribution f .
We first prove that the distribution of any one of the samples generated in Algorithm 1 is close to f CI in
terms of total variational distance. We make the following assumptions on the joint distribution of the original
samples i.e. fX,Y,Z (x, y, z):
Smoothness assumption on f (y|z): We assume a smoothness condition on f (y|z), that is a generalization
of boundedness of the max. eigenvalue of Fisher Information matrix of y w.r.t z.
5

Assumption 1. For z ∈ Rdz , a such that ka − zk2 ≤ 1 , the generalized curvature matrix Ia (z) is,
#
"
!
Z
f (y|z)
∂2
δ 2 log f (y|z 0 )
Z=z
log
Ia (z)ij =
=E −
f (y|z)dy
∂zi0 ∂zj0
f (y|z 0 )
δzi0 δzj0
z 0 =a
0

(1)

z =a

dz

We require that for all z ∈ R and all a such that ka − zk2 ≤ 1 , λmax (Ia (z)) ≤ β. Analogous assumptions
have been made on the Hessian of the density in the context of entropy estimation [12].
Smoothness assumptions on f (z): We assume some smoothness properties of the probability density
function f (z). The smoothness assumptions (in Assumption 2) is a subset of the assumptions made in [13]
(Assumption 1, Page 5) for entropy estimation.
Definition 2. For any δ > 0, we define G(δ) = P (f (Z) ≤ δ). This is the probability mass of the distribution
of Z in the areas where the p.d.f is less than δ.
Definition 3. (Hessian Matrix) Let Hf (z) denote the Hessian Matrix of the p.d.f f (z) with respect to z i.e
Hf (z)ij = ∂ 2 f (z)/∂zi ∂zj , provided it is twice continuously differentiable at z.
Assumption 2. The probability density function f (z) satisfies the following:
(1) f (z) is twice continuously differentiable and the Hessian matrix Hf satisfies kHf (z)k2 ≤ cdz almost
everywhere, where cdz is only dependent on the dimension.
R
(2) f (z)1−1/d dz ≤ c3 , ∀d ≥ 2 where c3 is a constant.
Theorem 1. Let (X, Y 0 , Z) denote a sample in U20 produced by Algorithm 1 by modifying the original sample
(X, Y, Z) in U1 , when supplied with 2n i.i.d samples from the original joint distribution fX,Y,Z (x, y, z). Let
φX,Y,Z (x, y, z) be the distribution of (X, Y 0 , Z). Under smoothness assumptions (1) and (2), for any  < 1 ,
n large enough, we have:
dT V (φ, f CI ) ≤ b(n)
s



1
1 β c3 ∗ 21/dz Γ(1/dz ) βG (2cdz 2 )
dz +2
2
+
exp
−
nγ
,
+
c

+
G
2c

.
d
d
d
z
z
z
2 4 (nγdz )1/dz dz
4
2
Here, γd is the volume of the unit radius `2 ball in Rd .
Theorem 1 characterizes the variational distance of the distribution of a sample generated in Algorithm 1
with that of the conditionally independent distribution f CI . We defer the proof of Theorem 1 to Appendix A.
Now, our goal is to characterize the misclassification error of the trained classifier in Algorithm 2 under both
H0 and H1 . Consider the distribution of the samples in the data-set Dr used for classification in Algorithm 2.
Let q(x, y, z|` = 1) be the marginal distribution of each sample with label 1. Similarly, let q(x, y, z|` = 0)
denote the marginal distribution of the label 0 samples. Note that under our construction,
 CI
f (x, y, z)
if H0 holds
q(x, y, z|` = 1) = fX,Y,Z (x, y, z) =
6= f CI (x, y, z) if H1 holds
q(x, y, z|` = 0) = φX,Y,Z (x, y, z)

(2)

where φX,Y,Z (x, y, z) is as defined in Theorem 1.
Note that even though the marginal of each sample with label 0 is φX,Y,Z (x, y, z) (Equation (2)), they are
not exactly i.i.d owing to the nearest neighbor bootstrap. We will go on to show that they are actually close
to i.i.d and therefore classification risk minimization generalizes similar to the i.i.d results for classification [4].
First, we review standard definitions and results from classification theory [4].
Ideal Classification Setting: We consider an ideal classification scenario for CI testing and in the process
define standard quantities in learning theory. Recall that G is the set of classifiers under consideration. Let
CI
q̃ be our ideal distribution for q given by q̃(x, y, z|` = 1) = fX,Y,Z (x, y, z), q̃(x, y, z|` = 0) = fX,Y,Z
(x, y, z)
and q̃(` = 1) = q̃(` = 0) = 0.5. In other words this is the ideal classification scenario for testing CI. Let
L(g(u), `) be our loss function for a classifying function g ∈ G, for a sample u , (x, y, z) with true label `.
In our algorithms the loss function is the 0 − 1 loss, but our results hold for any bounded loss function s.t.
|L(g(u), `)| ≤ |L|. For a distribution q̃ and a classifier g let Rq̃ (g) , Eu,`∼q̃ [L(g(u), `)] be the expected risk
6

of the function g. The risk optimal classifier gq̃∗ under q̃ is given by gq̃∗ , arg ming∈G Rq̃ (g). Similarly for
P
1
a set of samples S and a classifier g, let RS (g) , |S|
u,`∈S L(g(u), `) be the empirical risk on the set of
samples. We define gS as the classifier that minimizes the empirical loss on the observed set of samples
S that is, gS , arg ming∈G RS (g).
If the samples in S are generated independently from q̃, then standard results from the learning theory states
that with probability ≥ 1 − δ,
r
r
V
2 log(1/δ)
∗
Rq̃ (gS ) ≤ Rq̃ (gq̃ ) + C
+
,
(3)
n
n
where V is the VC dimension [30] of the classification model, C is an universal constant and n = |S|.
Guarantees under near-independent samples: Our goal is to prove a result like (3), for the classification
problem in Algorithm 2. However, in this case we do not have access to i.i.d samples because the samples in
U20 do not remain independent. We will see that they are close to independent in some sense. This brings us
to one of our main results in Theorem 2.
Theorem 2. Assume that the joint distribution f (x, y, z) satisfies the conditions in Theorem 1. Further
assume that f (z) has a bounded Lipschitz constant. Consider the classifier ĝ in Algorithm 2 trained on the
set Dr . Let S = Dr . Then according to our definition gS = ĝ. For  > 0 we have:
(i) Rq (gS ) − Rq (gq∗ ) ≤ γn
! 
!
!
r
1/3 r d
√
log(n/δ)
1
4 z log(n/δ) + on (1/)
, C|L|
+ G() ,
V + log
+
δ
n
n
with probability at least 1 − 8δ. Here V is the V.C. dimension of the classification function class, G is as
defined in Def. 2, C is an universal constant and |L| is the bound on the absolute value of the loss.
(ii) Suppose the loss is L(g(u), `) = 1g(u)6=` (s.t |L| ≤ 1). Further suppose the class of classifying functions is
such that Rq (gq∗ ) ≤ r0 + η. Here, r0 , 0.5(1 − dT V (q(x, y, z|1), q(x, y, z|0))) is the risk of the Bayes optimal
classifier when q(` = 1) = q(` = 0). This is the best loss that any classifier can achieve for this classification
problem [4]. Under this setting, w.p at least 1 − 8δ we have:
 b(n)
 b(n)
1
1
1 − dT V (f, f CI ) −
≤ Rq (gS ) ≤
1 − dT V (f, f CI ) +
+ η + γn
2
2
2
2
where b(n) is as defined in Theorem 1.
We prove Theorem 2 as Theorem 3 and Theorem 4 in the appendix. In part (i) of the theorem we prove
that generalization bounds hold even when the samples are not exactly i.i.d. Intuitively, consider two sample
inputs ui , uj ∈ U1 , such that corresponding Z coordinates zi and zj are far away. Then we expect the
resulting samples u0i and u0j (in U20 ) to be nearly-independent. By carefully capturing this notion of spatial
near-independence, we prove generalization errors in Theorem 3. Part (ii) of the theorem essentially implies
that the error of the trained classifier will be close to 0.5 (l.h.s) when f ∼ f CI (under H0 ). On the other
hand under H1 if dT V (f, f CI ) > 1 − γ, the error will be less than 0.5(γ + b(n)) + γn which is small.

4

Empirical Results

In this section we provide empirical results comparing our proposed algorithm and other state of the art
algorithms. The algorithms under comparison are: (i) CCIT - Algorithm 3 in our paper where we use
XGBoost [6] as the classifier. In our experiments, for each data-set we boot-strap the samples and run our
algorithm B times. The results are averaged over B bootstrap runs1 . (ii) KCIT - Kernel CI test from [32].
We use the Matlab code available online. (iii) RCIT - Randomized CI Test from [28]. We use the R package
that is publicly available.
1 The

python package for our implementation can be found here (https://github.com/rajatsen91/CCIT).

7

4.1

Synthetic Experiments

We perform the synthetic experiments in the regime of post-nonlinear noise similar to [32]. In our experiments
X and Y are dimension 1, and the dimension of Z scales (motivated by causal settings and also used in [32, 28]).
X and Y are generated according to the relation G(F (Z) + η) where η is a noise term and G is a non-linear
function, when the H0 holds. In our experiments, the data is generated as follows: (i) when X ⊥
⊥ Y |Z, then
each coordinate of Z is a Gaussian with unit mean and variance, X = cos(aT Z + η1 ) and Y = cos(bT Z + η2 ).
Here, a, b ∈ Rdz and kak = kbk = 1. a,b are fixed while generating a single dataset. η1 and η2 are zero-mean
Gaussian noise variables, which are independent of everything else. We set V ar(η1 ) = V ar(η2 ) = 0.25. (ii)
when X ⊥
6 ⊥ Y |Z, then everything is identical to (i) except that Y = cos(bT Z + cX + η2 ) for a randomly
chosen constant c ∈ [0, 2].
In Fig. 2a, we plot the performance of the algorithms when the dimension of Z scales. For generating each
point in the plot, 300 data-sets were generated with the appropriate dimensions. Half of them are according
to H0 and the other half are from H1 Then each of the algorithms are run on these data-sets, and the ROC
AUC (Area Under the Receiver Operating Characteristic curve) score is calculated from the true labels (CI
or not CI) for each data-set and the predicted scores. We observe that the accuracy of CCIT is close to 1 for
dimensions upto 70, while all the other algorithms do not scale as well. In these experiments the number
√
of bootstraps per data-set for CCIT was set to B = 50. We set the threshold in Algorithm 3 to τ = 1/ n,
which is an upper-bound on the expected variance of the test-statistic when H0 holds.

4.2

Flow-Cytometry Dataset

We use our CI testing algorithm to verify CI relations in the protein network data from the flow-cytometry
dataset [26], which gives expression levels of 11 proteins under various experimental conditions. The ground
truth causal graph is not known with absolute certainty in this data-set, however this dataset has been widely
used in the causal structure learning literature. We take three popular learned causal structures that are
recovered by causal discovery algorithms, and we verify CI relations assuming these graphs to be the ground
truth. The three graph are: (i) consensus graph from [26] (Fig. 1(a) in [22]) (ii) reconstructed graph by
Sachs et al. [26] (Fig. 1(b) in [22]) (iii) reconstructed graph in [22] (Fig. 1(c) in [22]).
For each graph we generate CI relations as follows: for each node X in the graph, identify the set Z
consisting of its parents, children and parents of children in the causal graph. Conditioned on this set Z, X
is independent of every other node Y in the graph (apart from the ones in Z). We use this to create all CI
conditions of these types from each of the three graphs. In this process we generate over 60 CI relations for
each of the graphs. In order to evaluate false positives of our algorithms, we also need relations such that
X⊥
6 ⊥ Y |Z. For, this we observe that if there is an edge between two nodes, they are never CI given any other
conditioning set. For each graph we generate 50 such non-CI relations, where an edge X ↔ Y is selected at
random and a conditioning set of size 3 is randomly selected from the remaining nodes. We construct 50 such
negative examples for each graph. In Fig. 2, we display the performance of all three algorithms based on
considering each of the three graphs as ground-truth. The algorithms are given access to observational data
for verifying CI and non-CI relations. In Fig. 2b we display the ROC plot for all three algorithms for the
data-set generated by considering graph (ii). In Table 2c we display the ROC AUC score for the algorithms
for the three graphs. It can be seen that our algorithm outperforms the others in all three cases, even when
the dimensionality of Z is fairly low (less than 10 in all cases). An interesting thing to note is that the edges
(pkc-raf), (pkc-mek) and (pka-p38) are there in all the three graphs. However, all three CI testers CCIT,
KCIT and RCIT are fairly confident that these edges should be absent. These edges may be discrepancies in
the ground-truth graphs and therefore the ROC AUC of the algorithms are lower than expected.

8

1.0
1.0

CCIT
RCIT
KCIT

0.8
True Positives

ROC AUC

0.9

0.8

0.7

0.6
0.4
CCIT
RCIT
KCIT

0.2
0.6
05

20

50
70
100
Dimension of Z

(a)

150

0.00.0

0.2

0.4
0.6
False Positives
(b)

0.8

1.0

Algo.

Graph (i)

Graph (ii)

Graph (iii)

CCIT
RCIT
KCIT

0.6848
0.6448
0.6528

0.7778
0.7168
0.7416

0.7156
0.6928
0.6610

(c)

Figure 2: In (a) we plot the performance of CCIT, KCIT and RCIT in the post-nonlinear noise synthetic data. In
generating each point in the plots, 300 data-sets are generated where half of them are according to H0 while the rest
are according to H1 . The algorithms are run on each of them, and the ROC AUC score is plotted. In (a) the number
of samples n = 1000, while the dimension of Z varies. In (b) we plot the ROC curve for all three algorithms based on
the data from Graph (ii) for the flow-cytometry dataset. The ROC AUC score for each of the algorithms are provided
in (c), considering each of the three graphs as ground-truth.

5

Conclusion

In this paper we present a model-powered approach for CI tests by converting it into binary classification,
thus empowering CI testing with powerful supervised learning tools like gradient boosted trees. We provide
an efficient nearest-neighbor bootstrap which makes the reduction to classification possible. We provide
theoretical guarantees on the bootstrapped samples, and also risk generalization bounds for our classification
problem, under non-i.i.d near independent samples. In conclusion we believe that model-driven data dependent
approaches can be extremely useful in general statistical testing and estimation problems as they enable us
to use powerful supervised learning tools.

References
[1] Maria-Florina Balcan, Nikhil Bansal, Alina Beygelzimer, Don Coppersmith, John Langford, and Gregory
Sorkin. Robust reductions from ranking to classification. Learning Theory, pages 604–619, 2007.
[2] Alina Beygelzimer, John Langford, Yuri Lifshits, Gregory Sorkin, and Alex Strehl. Conditional probability
tree estimation analysis and algorithms. In Proceedings of the Twenty-Fifth Conference on Uncertainty
in Artificial Intelligence, pages 51–58. AUAI Press, 2009.
[3] Karsten M Borgwardt, Arthur Gretton, Malte J Rasch, Hans-Peter Kriegel, Bernhard Schölkopf,
and Alex J Smola. Integrating structured biological data by kernel maximum mean discrepancy.
Bioinformatics, 22(14):e49–e57, 2006.
[4] Stéphane Boucheron, Olivier Bousquet, and Gábor Lugosi. Theory of classification: A survey of some
recent advances. ESAIM: probability and statistics, 9:323–375, 2005.
[5] Eliot Brenner and David Sontag. Sparsityboost: A new scoring function for learning bayesian network
structure. arXiv preprint arXiv:1309.6820, 2013.
[6] Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22Nd
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 785–794.
ACM, 2016.
9

[7] Jie Cheng, David Bell, and Weiru Liu. Learning bayesian networks from data: An efficient approach
based on information theory. On World Wide Web at http://www. cs. ualberta. ca/˜ jcheng/bnpc. htm,
1998.
[8] Hal Daumé, John Langford, and Daniel Marcu. Search-based structured prediction. Machine learning,
75(3):297–325, 2009.
[9] Luis M De Campos and Juan F Huete. A new approach for learning belief networks using independence
criteria. International Journal of Approximate Reasoning, 24(1):11–37, 2000.
[10] Gary Doran, Krikamol Muandet, Kun Zhang, and Bernhard Schölkopf. A permutation-based kernel
conditional independence test. In UAI, pages 132–141, 2014.
[11] Kenji Fukumizu, Francis R Bach, and Michael I Jordan. Dimensionality reduction for supervised learning
with reproducing kernel hilbert spaces. Journal of Machine Learning Research, 5(Jan):73–99, 2004.
[12] Weihao Gao, Sewoong Oh, and Pramod Viswanath. Breaking the bandwidth barrier: Geometrical
adaptive entropy estimation. In Advances in Neural Information Processing Systems, pages 2460–2468,
2016.
[13] Weihao Gao, Sewoong Oh, and Pramod Viswanath. Demystifying fixed k-nearest neighbor information
estimators. arXiv preprint arXiv:1604.03006, 2016.
[14] Markus Kalisch and Peter Bühlmann. Estimating high-dimensional directed acyclic graphs with the
pc-algorithm. Journal of Machine Learning Research, 8(Mar):613–636, 2007.
[15] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press,
2009.
[16] Daphne Koller and Mehran Sahami. Toward optimal feature selection. Technical report, Stanford
InfoLab, 1996.
[17] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional
neural networks. In Advances in neural information processing systems, pages 1097–1105, 2012.
[18] John Langford and Bianca Zadrozny. Reducing t-step reinforcement learning to classification. In Proc.
of the Machine Learning Reductions Workshop, 2003.
[19] David Lopez-Paz and Maxime Oquab.
arXiv:1610.06545, 2016.

Revisiting classifier two-sample tests.

arXiv preprint

[20] Colin McDiarmid. On the method of bounded differences. Surveys in combinatorics, 141(1):148–188,
1989.
[21] Mehryar Mohri and Afshin Rostamizadeh. Rademacher complexity bounds for non-iid processes. In
Advances in Neural Information Processing Systems, pages 1097–1104, 2009.
[22] Joris Mooij and Tom Heskes. Cyclic causal discovery from continuous equilibrium data. arXiv preprint
arXiv:1309.6849, 2013.
[23] Judea Pearl. Causality. Cambridge university press, 2009.
[24] V Ramasubramanian and Kuldip K Paliwal. Fast k-dimensional tree algorithms for nearest neighbor
search with application to vector quantization encoding. IEEE Transactions on Signal Processing,
40(3):518–531, 1992.
[25] Bero Roos. On the rate of multivariate poisson convergence. Journal of Multivariate Analysis, 69(1):120–
134, 1999.
[26] Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan. Causal proteinsignaling networks derived from multiparameter single-cell data. Science, 308(5721):523–529, 2005.
[27] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press,
2000.
10

[28] Eric V Strobl, Kun Zhang, and Shyam Visweswaran. Approximate kernel-based conditional independence
tests for fast non-parametric causal discovery. arXiv preprint arXiv:1702.03877, 2017.
[29] Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian
network structure learning algorithm. Machine learning, 65(1):31–78, 2006.
[30] Vladimir N Vapnik and A Ya Chervonenkis. On the uniform convergence of relative frequencies of events
to their probabilities. In Measures of Complexity, pages 11–30. Springer, 2015.
[31] Eric P Xing, Michael I Jordan, Richard M Karp, et al. Feature selection for high-dimensional genomic
microarray data. In ICML, volume 1, pages 601–608. Citeseer, 2001.
[32] Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. arXiv preprint arXiv:1202.3775, 2012.

11

A

Guarantees on Bootstrapped Samples

In this section we prove that the samples generated in Algorithm 1, through the nearest neighbor bootstrap,
are close to samples generated from f CI (x, y, z) = fX|Z (x|z)fY |Z (y|z)fZ (z). The closeness is characterized
in terms of total variational distance as in Theorem 1. Suppose 2n i.i.d samples from distribution f (x, y, z)
are supplied to Algorithm 1. Consider a typical sample (X, Y, Z) ∼ f (x, y, z), which is modified to produce a
typical sample in U20 (refer to Algorithm 1) denoted by (X, Y 0 , Z). Here, Y 0 are the Y -coordinates of a sample
(X 0 , Y 0 , Z 0 ) in U2 such that Z 0 is the nearest neighbor of Z. Let us denote the marginal distribution of a
typical sample in U20 by φX,Y,Z (x, y, z), i.e (X, Y 0 , Z) ∼ φX,Y,Z (x, y, z). Now we are at a position to prove
Theorem 1.
Proof of Theorem 1. Let fZ 0 |z (z 0 ) denote the conditional p.d.f of the variable Z 0 (that is the nearest neighbor
of sample Z in U2 ), conditioned on Z = z. Therefore, the distribution of the new-sample is given by,
Z
φX,Y,Z (x, y, z) = fX|Z (x|z)fZ (z) fY |Z (y|z 0 )fZ 0 |z (z 0 )dz 0 .
(4)
CI
We want to bound the total variational distance between φX,Y,Z (x, y, z) and fX,Y,Z
(x, y, z). We have the
following chain:

2 ∗ dT V (φ, f CI ) =

Z

fX|Z (x|z)fY |Z (y|z)fZ (z) − fX|Z (x|z)fZ (z)

x,y,z

Z 

Z
=
x,y,z

≤

fX|Z (x|z)fY |Z (y|z)fZ (z)

Z

Z

x,y,z

Z

0

=

Z

≤

x,z,kz 0 −zk2 ≤

Z

fX|Z (x|z)fZ (z)fZ 0 |z (z 0 )

Z



fY |Z (y|z 0 )fZ 0 |z (z 0 )dz 0 dxdydz
fZ 0 |z (z 0 )dz 0 dxdydz

fY |Z (y|z 0 )
fZ 0 |z (z 0 )dz 0 dxdydz
fY |Z (y|z)
!

fY |Z (y|z) − fY |Z (y|z 0 ) dy dz 0 dxdz

!
fY |Z (y|z) − fY |Z (y|z 0 ) dy dz 0 dxdz+

fX|Z (x|z)fZ (z)fZ 0 |z (z 0 )dz 0 dxdz

2
x,z,kz 0 −zk2 >

≤

Z

fX|Z (x|z)fZ (z)fZ 0 |z (z )

x,z,z 0

fY |Z (y|z 0 )
1−
fY |Z (y|z)

1−

fX|Z (x|z)fY |Z (y|z)fZ (z)

Z

Z
x,z,kz 0 −zk2 ≤

fX|Z (x|z)fZ (z)fZ 0 |z (z 0 )

Z

!
fY |Z (y|z) − fY |Z (y|z 0 ) dy dz 0 dxdz+

2 ∗ P (kz 0 − zk2 > )

(5)

By Pinsker’s inequality, we have:
Z
y

s Z
1
f (y|z)
fY |Z (y|z) − fY |Z (y|z ) dy ≤
log
f (y|z)dy
2 y
f (y|z 0 )
0

By Taylor’s expansion with second-order residual, we have:
Z
1
f (y|z)
log
f (y|z)dy = (z 0 − z)T Ia (z)(z 0 − z)
0
f (y|z )
2
12

(6)

(7)

for some a = λz + (1 − λ)z 0 where 0 ≤ λ ≤ 1.
Under Assumption 1 and  < 1 we have,
(z 0 − z)T Ia (z)(z 0 − z) ≤ βkz 0 − zk22 .

(8)

Then, (8), (7), (6) and (5) imply:
r
2 ∗ dT V (φ, f

CI

)≤

β
E[kz 0 − zk2 1kz0 −zk2 ≤ ] + 2P (kz 0 − zk2 > )
4

(9)

We now bound both terms separately. Let Z1 , Z2 ..., Zn be distributed i.i.d according to f (z). Then, fZ 0 |z (·)
is the pdf of the nearest neighbor of z among Z1 , . . . , Zn .
A.0.1

Bounding the first term

In this section we will use d in place of dz for notational simplicity. Let γd be the volume of the unit `2 ball
in dimension d. Let S = {z : f (z) ≥ 2 ∗ cd 2 }. This implies, that for z ∈ S:
f (z) − cd 2 ≥ f (z)/2

(10)

Let Z 0 = argminZ1 ,Z2 ...,Zn kZi − zk2 be the random variable which is the nearest neighbor to a point z
among n i.i.d samples Zi drawn from the distribution whose pdf is f (z) that satisfies assumption 2. Let
r(z) = ||z − z 0 ||2 . Let F (r) be the CDF of the random variable R. Since R is a non-negative random variable,
Z
Z
Z

ER [r(z)1r≤ ] =
rdF (r) = [rF (r)]0 −
F (r)dr ≤
P (R > r)dr
(11)
r≤

r≤

r≤

For any r ≤ , observe that
Pr(R > r) = Pr(@i : zi ∈ B(z, r))

n

= (1 − Pr(Z ∈ B(z, r)))

≤ exp(−nPr(Z ∈ B(z, r)))

(12)

We have the following chain to bound Pr(Z ∈ B(z, r)). Let a = λz + (1 − λ)t.
Z
d
|Pr(Z ∈ B(z, r)) − f (z)γd r | ≤
(f (t) − f (z))dt
t∈B(z,r)

Z
=
t∈B(z,r)

(∇f (z)T .(t − z) + (t − z)T Hf (a)(t − z))dt

≤ max kHf (t)k2
t∈B(z,r)

Z
t∈B(z,r)

kt − zk22 dt

≤ cd γd rd+2

(13)

By putting together (11),(12) and(13), we have:
Z

−nγd r d (f (z)−cd r 2 )
ER [r(z)1r≤ ]≤1z∈S
e
dr + 1z∈S c
r≤
Z

d
1
≤ 1z∈S
e− 2 nγd f (z)r dr + 1z∈S c
r≤

!
e−t
−1+1/d
≤ 1z∈S
t
dt + 1z∈S c
1
1/d
t≤ 12 nγd f (z)d d[ 2 nγd f (z)]


21/d
Γ(1/d)
+ 1z∈S c
≤ 1z∈S
d(nγd f (z))1/d
Z

13

Therefore, the first term in bounded by:
E[kz 0 − zk2 1kz0 −zk2 ≤ ] ≤ EZ [ER [r(z)1r≤ ]]


21/d
c
Γ(1/d) + 1z∈S
≤ EZ
d(nγd f (z))1/d

c3 ∗ 21/d
≤
Γ(1/d) +  ∗ G 2cd 2
1/d
(nγd ) d

(14)
(15)
(16)
(17)

A.0.2

Bounding the second term

We now bound the second term as follows:
Pr(||z − z 0 ||2 > ) ≤ EZ [Pr(R > )]



nγd f (z)d
≤ EZ 1z∈S exp −
+ Pr(z ∈ S c )
2


≤ exp −nγd cd d+2 /2 + G 2cd 2
Substituting in (9), we have:
s
2 ∗ dT V (g, f CI ) ≤



β c3 ∗ 21/d Γ(1/d) βG (2cd 2 )
d+2
2
+
+
2
exp
−nγ
c

/2
+
2G
2c

d
d
d
4 (nγd )1/d d
4

Substitute dz in place of d to recover Theorem 1.

B

Generalization Error Bounds on Classification

In this section, we will prove generalization error bounds for our classification problem in Algorithm 2.
Note the samples in U20 are not i.i.d, so standard risk bounds do not hold. We will leverage a spatial near
independence property to provide generalization bounds under non-i.i.d samples. In what follows, we will
prove the results for any bounded loss function L(g(u), `) ≤ |L|. Let S , Dr i.e., the set of training samples.
For 1 ≤ i ≤ 3, let Zi , {z : (x, y, z) ∈ Ui }. Let Z , Z1 ∪ Z2 . Observe that
Rq (gS ) ≤ Rq (gq∗ ) + 2 sup(RS (g) − Rq (g)),

(18)

g∈G

and hence in the rest of the section we upper bound supg∈G (RS (g) − Rq (g)). To this end, we define conditional
risk RS (q|Z) as
1 X
RS (g|Z) ,
E[L(g(u), `)|Z].
n
(u,`)∈S

By triangle inequality,
sup(RS (g) − Rq (g)) ≤ sup(RS (g) − RS (g|Z)) + sup(RS (g|Z) − Rq (g)).
g∈G

g∈G

g∈G

We first bound the second term in the right hand side of Equation (19) in the next lemma.
Lemma 1. With probability at least 1 − δ,
r
r
V
2 log(1/δ)
sup(RS (g|Z) − Rq (g)) ≤ |L|C
+ |L|
,
n
n
g∈G
where V is the VC dimension of the classification model.
14

(19)

Proof. For a sample u = (x, y, z), observe that Z → z → u forms a Markov chain. Hence,
E[L(g(u), `)|Z] = E[L(g(u), `)|Z].
Let
h(Z) , E[L(g(u), `)|Z]

(20)

Hence,
RS (g|Z) − Rq (g) =

1 X
h(Z) − Eq [h(Z)].
n
z∈Z

The above term is the average of n independent random variables h(Z) and hence we can apply standard
tools from learning theory [4] to obtain
r
r
Vh
2 log(1/δ)
sup(RS (g|Z) − Rq (g)) ≤ |L|C
+ |L|
,
n
n
g∈G
where Vh is the VC dimension of the class of models of h. The lemma follows from the fact that VC dimension
of h is smaller than the VC dimension of the underlying classification model.
We next bound the first term in the RHS of Equation (19). Proof is given in Appendix B.1.
Lemma 2. Let  > 0. If the Hessian of the density f (z)and the Lipscitz constant of the same is bounded,
then with probability at least 1 − 7δ,
! 
!
r
1/3 r d
√
log(n/δ)
1
4 log(n/δ) + on (1/)
sup(RS (g) − Rn (g|Z)) ≤ |L|
V + log
+
δ
n
n
g∈G
+ |L|G().

(21)

Lemmas 1 and 2, together with Equations (18) and (19) yield the following theorem.
Theorem 3. Let  > 0. If the Hessian and the Lipscitz constant of f (z) is bounded, then with probability at
least 1 − 8δ,
! 
!
!
r
1/3 r d
√
1
log(n/δ)
4 log(n/δ) + on (1/)
∗
V + log
+
+ G() , (22)
Rq (ĝ) ≤ Rq (gq ) + c|L|
δ
n
n
where c is a universal constant and ĝ is the minimizer in Step 6 of Algorithm 2.

B.1

Proof of Lemma 2

We need few definitions to prove Lemma 2. For a point z, let Bn (z) be a ball around it such that
2

Pr (Z ∈ Bn (z)) =

Z∼f (z)

log nδ
, αn .
n

Intuitively, with high probability the nearest neighbor of each sample z lies in Bn (z). We formalize it in the
next lemma.
Lemma 3. With probability > 1 − δ, the nearest neighbor of each sample z ∈ Z2 in Z3 lie in B(z).
Proof. The probability that none of Z3 appears in B(z) is 1 − (1 − αn )n ≤ δ/n. The lemma follows by the
union bound.
We now bound the probability that the the nearest neighbor balls Bn () intersect for two samples.
15

Lemma 4. Let  > 0. If the Hessian of the density (f (z)) is bounded by c and the Lipschitz constant is
bounded by γ, then for any given z1 such that f (z1 ) ≥  and a sample z2 ∼ f ,
Pr (Bn (z1 ) ∩ Bn (z2 ) 6= ∅) ≤ βn , 4d αn (1 + on (1/)).

z2 ∼f

Proof. Let rn (z) denote the radius of Bn (z). Let B(z, r) denote the ball of radius r around z and V (z, r) be
its volume. We can rewrite βn as
= Pr(Bn (z1 ) ∩ Bn (z2 ) 6= ∅)

= Pr(Bn (z1 ) ∩ Bn (z2 ) 6= ∅, 3rn (z1 ) ≥ rn (z2 )) + Pr(Bn (z1 ) ∩ Bn (z2 ) 6= ∅, 3rn (z1 ) < rn (z2 )).
We first bound the first term. Note that
Z

f (z 0 )dz 0 = αn ,

z 0 ∈Bn (z)

Hence, by Taylor’s series expansion and the bound on Hessian yields,
Z

αn =
f (z 0 )dz 0 = V (z, rn (z)) f (z) + O(rn2 (z)c) .
z 0 ∈Bn (z)

Similarly,

Pr(z 0 ∈ V (z, 4rn (z))) = V (z, 4rn (z)) f (z) + O(9rn2 (z)c)
= 4d αn (1 + on (1/)),

where the last equality follows from the fact that V (z, 4rn (z))/V (z, rn (z)) = 4d in d dimensions.
Then the first term can be bounded as
Pr(Bn (z1 ) ∩ Bn (z2 ) 6= ∅, 3rn (z1 ) ≥ rn (z2 )) = Pr(z2 ∈ B(z1 , rn (z1 ) + rn (z2 )), 3rn (z1 ) ≥ rn (z2 ))
≤ Pr(z2 ∈ B(z1 , 4rn (z1 )), 3rn (z1 ) ≥ rn (z2 ))
≤ Pr(z2 ∈ B(z1 , 4rn (z1 )))

≤ 4d αn (1 + on (1/)).

To bound the second term, observe that if Bn (z1 ) ∩ Bn (z2 ) 6= ∅ and 3rn (z1 ) < rn (z2 ). There exists a point z 0
on the line joining z1 and z2 at distance 3rn (z1 ) from z1 such that
Pr(z 00 ∈ B(z 0 , 3rn (z1 ))) < αn .
As before bound on the Hessian yields,
αn > V (z 0 , 3rn (z1 ))(f (z 0 ) − O(9rn2 (z1 )c)).
Hence,
f (z 0 ) < 3−d (f (z) + O(rn2 (z)c)) + O(9rn2 (z1 )c).
However, f (z) >  and f (z 0 ) ≥ f (z) − 3rn (z1 )γ and rn (z1 ) → 0. Hence, a contradiction. Thus
Pr(Bn (z1 ) ∩ Bn (z2 ) 6= ∅, 3rn (z1 ) < rn (z2 )) = 0.

Consider the graph on indices [n], such that two indices are connected if and only if Bn (zi ) ∩ Bn (zj ) 6= ∅,
f (z1 ) ≥ , f (z2 ) ≥  . Let ∆(Z1n ) be the maximum degree of the resulting graph. We first show that the
maximum degree of this graph is small.
16

Lemma 5. With probability ≥ 1 − δ,

∆(Z1n ) ≤ 4nβn .

Proof. For index 1, by Lemma 4 that probability of j points intersect is at most
j  
X
n i
βn (1 − βn )n−i .
i
i=0
Hence, the degree of vertex 1 is dominated by a binomial distribution with parameters n and βn . The lemma
follows from the union bound and the Chernoff bound.
Let k > 2∆(Z1n ) and S0 , S1 , S2 . . . Sk be k independent sets of the above graph such that maxt≥1 |St | ≤ 2n/k.
Note that such independent sets exists by Lemma 10. We set the exact value of k later. Let S0 contains all
indices such that f (zi ) < .
Lemma 6. With probability > 1 − δ,
r
1
|S0 | ≤ nG() + n log .
δ
Proof. Observe that |S0 | is the sum of n independent random variables and changing any of them changes S0
by at most 1. The lemma follows by McDiarmid’s inequality.
We can upper bound the LHS in Equation (21) as
sup(RS (g) − Rn (g|Z)) ≤
g∈G

k
X
|St |

sup

t=1

g∈G

n

1 X
|S0 |
|L|.
(L(g(ui ), `i ) − h(Zi )) +
|St |
n
i∈St

Let N (Zi ) denote the number of elements of Z3 that are in B(Zi ) and Let Ai be always true if Zi ∈ Z1 and
otherwise Ai be the event such that nearest neighbor of samples in N (Zi ) > 0. We first show the following
inequality.
Lemma 7. With probability ≥ 1 − δ, for all sets St .
1 X
1 X
(L(g(ui ), `i ) − h(Zi )) = sup
1Ai (L(g(ui ), `i ) − h(Zi )) .
sup
g∈G |St |
g∈G |St |
i∈St

i∈St

Proof. Let Xi = (L(g(ui ), `i ) − h(Zi )) and Observe that LHS can be written as
1 X
1 X
1 X
sup
Xi = sup
Xi 1Ai + sup
(Xi − Xi 1Ai ).
g∈G |St |
g∈G |St |
g∈G |St |
i∈St

i∈St

i∈St

If the conditions of Lemma 3 hold, the nearest sample of Zi ’s lie within Bn (Zi ). Hence, with probability
≥ 1 − δ, the second term is 0. Hence the lemma.
Let r be defined as follows.
r(Zi , Ni ) , E[1Ai (L(g(ui ), `i ) − h(Zi )) |N (Zi ) = Ni ].
Observe that
"
E

#
X

i∈St

=

X
i∈St

=

X
i∈St

=

X

1Ai (L(g(ui ), `i ) − h(Zi )) |N (Z1 ), . . . N (Zn )
E [1Ai (L(g(ui ), `i ) − h(Zi )) |N (Z1 ), . . . N (Zn )]
E [1Ai (L(g(ui ), `i ) − h(Zi )) |N (Zi )]
r(Zi , Ni ).

(23)

i∈St

17

Hence, we can split the term as
1 X
1Ai (L(g(ui ), `i ) − h(Zi ))
sup
g∈G |St |
i∈St
1 X
1 X
≤ sup
1Ai (L(g(ui ), `i ) − r(Zi , N (Zi ))) + sup
1Ai (r(Zi , N (Zi )) − h(Zi ))
g∈G |St |
g∈G |St |
i∈St

(24)

i∈St

Given {Zi , N (Zi )}, the first term in the RHS of the Equation (24) is a function of |St | independent random
variables as 1Ai ∗ L(g(ui ), `i ) are mutually independent given {Zi , N (Zi )}. Thus we can use standard tools
from VC dimension theory and state that with probability ≥ 1 − δ, the first term in the RHS of Equation (24)
can be upper bounded as
s
s
V
log(1/δ)
1 X
1Ai (L(g(ui ), `i ) − r(Zi , N (Zi ))) ≤ |L|C
sup
+ |L|
.
|S
|
|S
|
|St |
t
t
g∈G
i∈St

conditioned on {Zi , N (Zi )}
To bound the second term in the RHS of Equation (24), observe that unlike the first term, the N (Zi )s are
dependent on each other. However note that N (Z1 ), . . . N (Z|St | ) are distributed according to multinomial
distribution with parameters n and αn . However, if we replace them by independent Poisson distributed
N (Zi )s we expect the value not to change. We formalize it by total variation distance. By Lemma 9, the
total variation distance between a multinomial distribution and product of Poisson distributions is
O(|St |αn ),
and hence any bound holds in the second distribution holds in the first one with an additional penalty
of
O (|St |αn |L|) .
Under the new independent sampling distribution, again the samples are independent and we can use standard
tools from VC dimension and hence, with probability ≥ 1 − δ, the term is upper bounded by
s
s
V
log(1/δ)
|L|C
+ |L|
.
|St |
|St |
Hence, summing over all the bounds, we get

s

s
k
1
X
|S0 |
V
|St |  log δ
sup(RS (g) − Rn (g|Z)) ≤ |L|O 
+
+c
+ αn |St |
n
n
|S
|
|S
t
t|
g∈G
t=1


s
!
r
r
1
√
log δ
k
1

≤ |L|O G() +
log + c V + αn max |St | +
t
n
δ
n
r

s
!
r
1
√
log δ
k
1
n
.
≤ |L|O 
log + c V + αn +
n
δ
k
n
2/3

conditioned on Zi , N (Zi ) Choose k = nαn + 8nβn , and note that the conditioning can be removed as the
term on the r.h.s are constants.This yields the result. The error probability follows by the union bound.
Theorem 4. Assume the conditions for Theorem 3. Suppose the loss is L(g(u), `) = 1g(u)6=` (s.t |L| ≤
1). Further suppose the class of classifying function is such that Rq (gq∗ ) ≤ r0 + η. Here, r0 , 0.5(1 −
dT V (q(x, y, z|1), q(x, y, z|0))) is the risk of the Bayes optimal classifier when P(` = 1) = P(` = 0). This is the
best loss that any classifier can achieve for this classification problem [4]. Under this setting, w.p at least
1 − 8δ we have:
 b(n)
 b(n)
1
1
1 − dT V (f, f CI ) −
≤ Rq (gS ) ≤
1 − dT V (f, f CI ) +
+ η + γn
2
2
2
2
18

Proof. Assume the bounds of Theorem 3 holds which happens w.p at least 1 − 8δ. From Theorem 3 we have
that
Rq (gS ) ≤ Rq (gq∗ ) + γn .
(25)
Also, note that from Theorem 1 we have the following:
dT V (q(x, y, z|1), q(x, y, z|0)) = dT V (φ, f )
≤ dT V (φ, f CI ) + dT V (f CI , f )
≤ b(n) + dT V (f CI , f )

(26)

Under our assumption we have Rq (gq∗ ) ≤ r0 + η. Combining this with (25) and (26) we get the r.h.s. For, the
l.hs note that Rq (gq∗ ) ≥ r0 as the bayes optimal classifier has the lowest risk. We can now use (26) to prove
the l.h.s.

C

Tools from probability and graph theory

Lemma 8 (McDiarmid’s inequality [20]). Let X1 , X2 , . . . Xm be m independent random variables and f be a
function from xn1 → R such that changing any one of the Xi s changes the function f at most by ci , then


−22
Pr(f − E[f ] ≥ ) ≤ exp Pm 2 .
i=1 ci
Lemma 9 (Special case of Theorem 1 in [25]). Let fm be the multinomial distribution with parameters n
Pk
and p1 , p2 , . . . pk , 1 − i=1 pi , and fp be the product of Poisson distributions with mean npi for i ≤ 1 ≤ k,
then
k
X
dT V (fm , fs ) ≤ 8.8
pi .
i=1

Lemma 10. For a graph with maximum degree ∆, there exists a set of independent sets S1 , S2 , . . . Sk such
that k ≥ 2∆ and
max |Si | ≤ 2n/k.
1≤i≤k

Proof. We show that the following algorithm yields a coloring (and hence independent sets) with the required
property.
Let 1, 2, . . . k be k colors, where k > 2∆. We arbitrarily order the nodes, and sequentially color nodes with a
currently least used color from among the ones not used by its neighbors. Consider the point in time when i
nodes have been colored, and we evaluate the options for the (i + 1)th node. The number of possible choices
of color for that node is c ≥ k − ∆. Out these c colors, the average number of nodes belonging to each color at
this point is at-most i/c. Therefore by pigeonholing, the minimum is less than the average; thus the number
of nodes belonging to chosen color is no larger than i/c ≤ i/(k − ∆).
Hence at the end when all n nodes are colored, each color has been used no more than (n − 1)/(k − ∆) + 1 <
2n/k.

19

