K-Nearest-Neighbor Local Sampling Based Conditional Independence Testing
Authors: Shuai Li
Yingjie Zhang
Hongtu Zhu
Christina Dan Wang
Hai Shu
Ziqi Chen
Zhuoran Sun
Yanfeng Yang
Venue: NeurIPS 2023 poster (2023)
Keywords: Conditional Independence testing
causal inference
conditional mutual information
k-nearest neighbor
conditional randomization test
conditional permutation test
TL;DR: Local Sampling Based CIT

Abstract: Conditional independence (CI) testing is a fundamental task in statistics and machine learning, but its effectiveness is hindered by the challenges posed by high-dimensional conditioning variables and limited data samples. This article introduces a novel testing approach to address these challenges and enhance control of the type I error while achieving high power under alternative hypotheses. The proposed approach incorporates a computationally efficient classifier-based conditional mutual information (CMI) estimator, capable of capturing intricate dependence structures among variables. To approximate a distribution encoding the null hypothesis, a $k$-nearest-neighbor local sampling strategy is employed. An important advantage of this approach is its ability to operate without assumptions about distribution forms or feature dependencies. Furthermore, it eliminates the need to derive asymptotic null distributions for the estimated CMI and avoids dataset splitting, making it particularly suitable for small datasets. The method presented in this article demonstrates asymptotic control of the type I error and consistency against all alternative hypotheses. Extensive analyses using both synthetic and real data highlight the computational efficiency of the proposed test. Moreover, it outperforms existing state-of-the-art methods in terms of type I and II errors, even in scenarios with high-dimensional conditioning sets. Additionally, the proposed approach exhibits robustness in the presence of heavy-tailed data.

K K-Nearest-Neighhor Local Sampling Based Conditional Independence Testing
Shuai Li, Yingjie Zhang
Hongtu Zhu Departments of Biostatistics, Statistics, Computer Science, and Geneties The University of North Carolina at Chapel Hill
School of Statistics KLATASDS-MOE
East China Normal University
Christina Dan Wang
Hai Shu Department of Biostatistics School of Global Public Health New York University
Business Division
New York University Shanghai
Ziqi Chen; Zhuoran Sun, Yanfeng Yang
School of Statistics
KLATASDS-MO
East China Normal University
Abstract
Conditional independence (CI) testing is fundamental task in statisties and machine learning, but its effectiveness hindered by the challenges posed by high-dimensional conditioning variables and limited data samples. This article introduces a novel testing approach to address these challenges and enhance con- trol of the type I eгтor while achieving high power under alternative hypotheses. The proposed approach incorporates a computationally efficient classifier-basee conditional mutual information (CMI) estimator, capable of capturing intricate dependence structures among variables. To approximate distribution encoding the null hypothesis, k-nearest-neighhorr local sampling strategy employed. An important advantage of this approach is its ability to operate without assumptions about distribution forms or feature dependencies. Furthermore, eliminates the need to derive asymptotic null distributions for the estimated CMI and avoids dataset splitting, making particularly suitable for small datasets. The method presented in this article demonstrates asymptotic control of the type error and consistency against all alternative hypotheses. Extensive analyses using both syn- thetic and real data highlight the computational efficiency of the proposed test. Moreover, it outperforms existing state-of of-the-art methods in terms of type and II error,, even in scenarios with nigh-dimensional conditioning sets. Additionally, the proposed approach exhibits robustness the presence of heavy-tailed data.
Introduction
Testing for conditional independence (CI) is crucial and challenging task in statistics and machine learning, with wide-ranging applications graphical models (26] [7], causal inference 44 37/ 22). and variable selection [I]] 24]. The objective is to determine whether two random variables, X and are independent given set of conditioning variables Z, based on observations of the joint
"Corresponding author (zgchen@fem.ecmm @ww..ccu.emmc..cc....
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

distribution (x,,, Specifically, the hypothesis to be tested is:
Ho: X TT Z versus H1 X 1o Z
(I)
However, CI testing becomes challenging due to the high-dimensionality of the conditioning variables Z, and the limited availability of data samples [5] 35 40, []. Existing tests may struggle to control the type I error, particularly when handling digennimensional conditioning variable sets with complex dependency structures 40). Moreover, even when a test is valid, the limited data availability can make it exceedingly challenging to distinguish between null and alternative hypotheses, resulting in low testing power [40].
In this article, present novel conditional independence (CI) testing method based k-nearest- neighbor local sampling scheme. Our approach incorporates two essential components: a conditional mutual information (CMI) estimator utilizing classification techniques, and k-nearest-neighhor local sampling strategy to approximate the conditional distribution Pxiz (z). This approximation enables us to simulate data sets from distribution that represents the null hypothesis, allowing us to estimate the CMI using these simulated data sets. By comparing the CMI estimator based on real data with those computed from the simulated data sets, we make informed decisions regarding the hypothesis testing. Theoretical analysis demonstrates that our proposed test achieves valid control of the type error asymptotically and exhibits consistency against all alternatives in H1. Synthetic and real data analyses showcase that our test outperforms previous methods in terms of both the type error and power under H1, even when dealing with high-diig high-ddmensional conditioning sets and/or small datasets. Furthermore, our approach remains computationally efficient as the dimension of the conditioning set and/or sample size increase. Additionally, our method is robust even the presence of heavy-tailed data.
2 Related Work
2.1 CMI estimation
In our work, we utilize CMI as a measure of conditional dependence. CMI offers a strong theoretical guarantee for CI testing, where IZ) 0 & X IL Y [4]. has the ability to capture complex dependence structures among variables, even in non-linear scenarios. Estimators based on k-nearest neighbors and kernel methods have been widely employed for CMI estimation
27 38 23 19. However, these methods may encounter efficiency issues when dealing with high-
dimensional conditioning variables, known as the curse of dimensionality [34 32. 32 approached the CMI estimation problem by formulating it minimax optimization problem and proposed the use of generative adversarial networks (GANs) to optimize and obtain the CMI estimator, which can handle -igh-dimmnsional scenarios. However, the training of GANs is often challenging, with the risk of collapse if hyperparameters and regularizers are not carefully chosen [15). Recently, [34 30 proposed classifier- CMI estimator capable of handling high dimensions. They initially developed an estimator for the Bullback-Leeibed (KI.) divergence using a classifier and then derived mutual information (MT) estimators based on it. CMI is defined as the difference between two MI values, and the CMI estimator is obtained by computing the difference between two MI estimators. In contrast, our method directly utilizes the classifier-based KL-divergence estimator obtain the CMI estimator. This approach also handles high dimensions effectively and is computationally more efficient compared the method employed by [34]] 30
2.2 CI Testing
In recent years, a considerable body of literature on CI testing has emerged. Here, we provide a brief overview of some existing methods, and for more comprehensive review, we refer readers to (28
One important category of CI testing methods involves using different measures for CI S 46
(18) proposed a CI testing method based on the empirical Hilbert-Schmidt norm of the conditional cross-covariance operator. Another approach, KCIT, was introduced by 47 which utilizes the partial association of reproducing kernel Hilbert spaces to measure conditional independence. (40] proposed a novel kernel-based CI testing method using the I, distance between two well-chosen analytic kernel mean embeddings evaluated at finite set of locations. 42) proposed the generalized covariance measure for testing based on a regression method. However, obtaining the exact distribution of the test statistic derived from the conditional independence measure under
For example.

---

Hu can be challenging. Instead, researchers have employed alternative methods to derive asymptotic distributions of the test statistics under the null hypothesis 47 effectiveness of asymptotic distributions can be compromised when the sample size is small or the dimension of Z is high [T6 38 As a result, tests based on asymptotic distributions may exhibit inflated type-I errors or inadequate power when dealing with small sample sizes or high-dimensional Z [4 In our work, we choose to utilize CMI as a metric for conditional dependence.
46
45
40
Nevertheless, the
To overcome the reliance on asymptotic null distributions, M proposed the model-X framework, assuming knowledge of the conditional distribution px (z) Under this assumption, a set of test statistics can be computed that are exchangeable under the null hypothesis, either through direct resampling [TI] or permutation methods [8]. However, in practice. the distribution of X Z is rarely available. Therefore, accurately approximating the distribution of X Z becomes crucial for maintaining the type error control of these tests. To address this challenge, [5] developed the generative conditional independence test (GCIT) using Wasserstein GANs to approximate the distribution of X IZ, while [4] proposed using Sinkhorn GANs for the same purpose. However, as highlighted in [5 30), limited data and noise may lead to inaccurate learning of conditional distributions using neural networks, resulting in inflated type errors for GCIT and the double GANs- based CI test (DGCIT). To mitigate this issue, [30] proposed using the 1 l-nearest neighbor method to generate samples from the approximated conditional distribution of X given Z. However, their method requires splitting the dataset into two parts, with the testing dataset used for calculating the test statistics comprising only one-third of the total samples. This will reduce the statistical power of the test, particularly when working with small datasets. [88 utilized Kozachenko- Lorzchenko-Leonenko estimator for CMI as the test statistic and proposed a permutation scheme based on k-nearest neighbors to generate samples under the null distribution. However, the curse of dimensionality adversely affects CMI estimator, leading to poor performance when the conditioning variable set Z has high dimensions. Additionally, no theoretical guarantee is provided. 31 and (25) employed a binning strategy, discretizing Z into finite number of bins based on the proximity of conditioning variables Z, followed by the ppermute-within -ithitggg---g strategy. However, selecting bins in high- sigh-ddmmnsionall settings presents a significant challenge [8. and their methods are limited handling conditioning variables with very few dimensions.
the performance of the
groups"
In our work, we propose the utilization of k-nearest- neiuhstr local sampling strategy as an appealing alternative the binning strategy. Specifically, our strategy generates samples locally based on the k-nearest neighbors of the conditioning variables Z. We demonstrate that the distribution of samples generated from this n-neaorest-neighbor local sampling scheme closely approximates the true conditional distribution 2x|2 (z]=) (x]2(x)22) in terms of total variation distance. One significant advantage of our proposed method is that allows the entire dataset to be used for computing the testing statistics, eliminating the need for dataset splitting. This feature makes our method more effective compared to 30). which requires dataset splitting, particularly when dealing with small datasets. Moreover, our method does not require the derivation of asymptotic null distributions for the estimated CMI, and it can easily handle nigh-dimensional conditioning variables Z. We provide theoretical and empirical evidence that our test achieves a valid control of the type I error and attains high power under the alternative hypothesis H1.
3
K K--eearstt-Nighhorr Local Sampling Based CI Testing
3.1
Classifier-based CMI estimator
The CMI for triplet of random variables/vectors (X, Y. Z) is defined as:
1(x; = JIʃ Px. log Px,Y.2 y,2)
dadydz,
px,z (x, 2)pr12(y)=)
where Px,Y. y. is the joint density of (X, Y, 2). Px,z(x, 2) is the joint density of (X, Z). and py12 (2)2(1)2 22) is the conditional density of Y given Z = 2. One approach to estimate CMI is by directly estimating the joint and conditional densities from the available data and plugging them
into (2) However, accurately estimating the density functions can be challenging, especially in
high- nigh-dimensional settings where it is more difficult than directly estimating CMI 2 special case of the Kullback-Leinn (KL) divergence, and thus we have:
[34] CMI is
(()) = DRI(px.Y.zZz, y, 2) (px. 2)priz (v)2)),
3

---

Algorithm 1 I-Nearest-Neighhor sampling (I-NN(V,, V2, Datasets Vi and V2, both with sample size and V Vi U V2 consisting of 2n independently and identically distributed (i.i.d.) samples from px,Y 2 (I. y, 2)
Input:
Output: A new data set V" consists of n samples.
Let V 0.
2: for (Xx, Y. 2) in V/ do
3: Go to V1 find the sample (X', YY,," Z') such that Z' the 1-nearest neighbor of Z in terms
of the 12 norm
4: V V" ((X,Y",
5: end for
6:
return V'
where ()lgg denotes the KI. divergence between two distribution functions F and G, with den-
sity functions f((x) and g(r), respectively. The Varadhan (DV) representation of Dx (fl given by:
sup (sw) - log Ew-g (w))),
(4)
where the function class S includes all functions with finite expectations. The optimal function in given by s (a:) = (x)))(())) 4. which leads to:
(fll))))) = E log (f(w)/9(w)) - log E (f(w)/o(w))
(5)
Building upon the classifier-based CMI estimation method proposed by [3 we propose CMI estimation method using the 1-nearest-neighee docode for estimating CMI is outlined in Algorithm 2 The main objective to empirically es- timate 5 with f Px,y, y, 2) and g px (x, px122 which requires samples from both (x, 2) and Px,z (2r 2) 2) (yl z). The available data only consists of samples from Px,Y.z y, =). However, generating samples from Px.z(2 2)Pr12 (vlz) requires knowledge of the unknown conditional distribution py)z ()) To address this challenge, we propose utilizing the 1-NN sampling algorithm [4] 30), which is outlined in Algorithm Empirical and theoretical results presented in (30 indicate that the I-NN sampling algorithm can accurately approximate the conditional distribution.
(1-NN) sampling algorithm 41 30
The pseu-
Next, we formalize the classifier-based CMI estimator. We consider a data set V consisting of 2n i.idd samples {W; := (X,, Y, 2, with (X,, Y,, Zz) Px,Y (zz y, z). The data set V is divided into two equally sized parts Vi and V2, where |V = |V/2 = n. For data sets Vi and V2, we use the I-NN data set V with samples. We assign labels / = 1 for all samples in V22 and 0 for all samples in V" In this supervised classification task, a binary classifier can be trained using advanced binary classification model, such XGBoost [41 [1] or deep neural networks [21] The classifier produces predicted probability am = P(( 1|,m) for a given
sampling algorithm D to generate
sample W.. leading to an estimator of the likelihood ratio on Wm. given by Ż(W... = Qm/ (1 am)
follows from 33 and (5) that an estimator of I (X; Z) given by
T(X;YIZ) D& Dww.px,, y. px px, 20)
= d log i(w:) - log [d Σi(w))
(6)
where d = [n/3] with [t] being the largest integer not greater than t, W is sample in vlest and W is a sample in test where v/lst and yyest are defined in Algorithm According to Theorem
1 in 34.. Y|2; is a consistent estimator of (Z).
In contrast to the approach taken by 34, 30 which estimate CMI by using the difference between two mutual informations (i.e., I(X; Y|) I(X; Y. Z) ) I(X; Z)) and therefore require training two binary classifiers, our method achieves CMI estimation with just single binary classifier. This leads to significantly improved computational efficiency. As depicted in Algorithm calculating a single p-value requires the computation of (B + 1) CMIs. Hence, in practical applications involving real data analysis, especially when dealing with large sample sizes or high-dimensional conditional variables, efficient computation of the CMI becomes crucial.

---

Algorithm 2 Classifier-based CMI Estimator
Input: Dataset V containing 2n i.idd. samples drawn from px,Y (2, y, 2) Output: An estimator of CMI.
1: Equally split V into two parts V1 and V2, each containing samples.
2: Apply Algorithm to generate a new dataset V with -=
3: Form the labeled datasets V = (WW/ 1): w E V23 and V = {(W), wy E V*).
4: Divide V, into training and testing subsets vyrain and vjest, at a ratio of 2:1.
5: Similarly, split Vg into training and testing subsets v-train and vttst at a ratio of 2:1.
6: Merge the datasets to form vtrain = vtrain ytrain and ytest = vjest u ylest
7: Train the classifier C using Vtrain
8: For each E Vifest, where Visatt includes all features in ylest, compute the classifier- based predicted probability P(l 1/w).
9: Calculate [(X; 2 as per formula 6).
10: return ((X;YIZ)
3.2 The Proposed CI Testing Procedure
While we have proposed a consistent estimator for CMI, accurately estimating as zero under the null hypothesis, Ho. is practically unattainable due to random errors in the sample data. These errors lead to deviations the estimator from the actual value. To enhance the effectiveness of CI testing based on the CMI estimator, we propose a test that takes into account the statistical variation inherent in the CMI estimator.
Under Ho, we can express the following representation:
Px (zly, 2) = px|z(z]2).
(7)
The null model in (7) preserves the dependence between X and Z, while breaking any dependence between X and Y. Therefore, if direct causal link exists between X and Y, replacing X with a null sample X ~ (x]z(>]2) would disrupt this relationship. Thus, we can conclude that X and Y are conditionally independent given Z, i.e., T(X"; =0
n) Consider i.i.d. copies A = { Mi, 2) : :i 1, n. of (X,Y, Z). When the distribution of X is known, conditional on Z = (21...., )m, we can independently draw a pseudo sample (b) Px (=]22) for each i across yT)TT independent of Y (V1,..., " and also X = (#) DT, where B is the number of
(b) 2 y are I..... B such that all X (z) .....
z
repetitions. Under Ho. (x((), Y, 2) (X, Y, Z) for all 6, where denotes equality in distribution. We denote the CMI estimator of I [(X; Z) based on (x(1), Y. Z) as CMI and denote the estimator based on (X, Y, Z) as CMI. We can approximate the p-value by
+ Σ ((CM 2 CMI)
p
(8)
1 B
where 1(-) is the indicator function. To begin with, we demonstrate that the test based on (8) can effectively control the type error. Specifically, under the null hypothesis Ho. the (B + 1) triples ( Y, Z) are exchangeable, and thus the p-value is valid and satisfies P(p S a Ho) I for any given a (0, ( Furthermore, it intuitively suggests that our test can achieve high power under the alternative hypothesis H1. As approaches infinity, (b) converges to zero in probability. But under H1, we know that I(X; Y|Z)= > 0. As result, CMI should be positive with high probability. Consequently, the p-value calculated using (8) is very small with high probability, indicating the consistency of our test against all alternatives stated in Hi.
(X, Y, Z). (x Y. Z)..... (x(
We propose using the k-nearest-neighbor local sampling strategy to approximate px (2 The approximated distribution is denoted as Pxiz (x) 2). To generate pseudo samples from pxdz (x) 22), we follow the steps below:
However, Px]z(x') 2)) is rarely known in practice 5 43
5

---

Algorithm 3 K T-Neearest-eiihhoor local sampling based CI testing
Input: Dataset (X, Y, Z) consisting of n i..d. samples from px.y, (x, Parameter: The number of repetitions B: the neighbor order k; the significance level a. Output: Accept Ho TT Z or H1 X ff Z
Use Algorithm 2 obtain CMI based on
Z).
2: for E {1, 21, ..... do
3:
Obtain the set of indices of the k-nearest neighbor of 2; M, {j (1.2...,,1) n..,.
|| - - ) where dţ denotes the distance of 2 its k-nearest neighbor.
4: end for
5: b=
6: while < B do
7: Initialize empty arra X of length n.
8: for E (1.2....,)) n. do
9: Shuffle M₁
j= [1]
11:
X X
12:
end for
13: Use Algorithm to obtain (M) based on 2).
6 b
15:
end while
Compute p-value: p:= [1+ E 1/CMI() 2 CMI))/( /11 + B).
17: if 2 then
Accept Ho X Z
19: else
20:
Accept H1 : X M4 Y
21: end if
1. Obtain the set of indices of the k-nearest neighbors of based the la norm. Denote this set as M, = {j E [I,...,.] ||| zill2 S di where dl is the distance of z to k-nearest neighbor;
2. Shuffle M, and let j be the first element of M; Then, define the pseudo sample Te as æyj
Repeat the above two steps for all E 21.2...., n to obtain the resultant pseudo sample vector X ",)"".
We establish in Section 4 that the total variation distance between the true distribution of X X and the distribution of samples generated by the neeerrsttneighbor local sampling strategy tends to zero in probability as n goes infinity. This evidence suggests that the latter close approximation of the former. Thus, (X, Y, Z) approximates directly a distribution that encodes the null hypothesis. We then compute the CMI estimator based on it. This process is repeated B times, resulting in B realizations of the CMI estimator under the null hypothesis denoted by (CM ..... (M) We can determine whether to reject the null hypothesis by comparing these estimators to the one of the original sample set A. Specifically, by (t) If the resulting p-value is smaller than the prespecified significance level, reject the null hypothesis. The pseudo code is presented in Algorithm 3 In Section 4 we will demonstrate through theoretical analysis that our test asymptotically achieves a valid control of type error and is consistent against all the alternatives in H,
calculate the p-value using (8) but with (M) replaced
30 employed the 1-NN sampling strategy approximate the distribution of X Z. To ensure substantial dissimilarity between the pseudo-sample datasets generated across repetitions, their method requires dividing the dataset into two parts. As consequence, the dataset used for calculating the test statistics consists of only one-third of the total samples. This limitation can result in reduced statistical power of the test, especially when working with small datasets. In contrast, our proposed procedure eliminates the need for dataset splitting and allows the entire dataset to be used in computing CMI. thereby avoiding the loss of testing power.

---

4 Theoretical Results
In this section, present our main theoretical results. All the detailed proofs for these results can be found in the Supplementary Materials. Let P1 and P2 be two probability distributions defined on the same probability space. The total variation distance between P1 and P3 defined as drv (Pr, = supAco (A() (2)),, where the supremum is taken over all measurable subsets A of the sample space f. We will show in Theorem 2 that the distribution of X generated by the d-nearest-neighhbor local sampling strategy is very close to the true conditional distribution in terms of the total variation distance.
demonstrates that the k-nearest neighbor of Z has similar property to the work of [3), where it has been shown that the nearest neighbor of Z converges almost surely (a.s.) to Z as the sample size n approaches infinity.
Lemma plays crucial role in proving Theorem 2
Lemma 1. Let Z and Z....., Zn be i.i.d. random vectors from p(z). For given positive integer k. let h be the k-nearest neighbor of 7. from the set ( 7
Then, 7 7 as n + 0c.
We make the following regularity conditions, which have been introduced in 19 20 and [11] To simplify the notation, we may drop the subscripts when referring to probability density functions.
Assumption 1 (Smoothness on p(æ[z)). For any 2 E id and any such that la 2 =11 S we assume S ( S B. where d2 the dimension of the random vector Z, |l-112 denotes the 12 norm, 0 < 3 < 00, and Ia (2), called the generalized curvature matrix, is defined with entries
Ta(=);; = (- 82 log p( Z = 0z,0z,
82
p(x12) p(x13) p(alz)do
7 z=
log
E=o
02,02,
12=a
Assumption 2 (Smoothness on p(z)). Assume that the probability density function p(z) is twice continuously differentiable. Let (2) denote the Hessian matrix of p(z) with respect to 2. We assume that || H>)2)]l2 S Cd. holds almost everywhere, where ca, depends only on the dimension d2 of the random vector Z.
Note that Assumption I can be viewed as an extension of the requirement that the maximum eigenvalue of the Fisher information matrix for with respect to z bounded. Assumptions I and 2 can be validated when (X, 2) follows multivariate Gaussian distribution (MVD) and when Z follows MVD, respectively.
For Z, we define Z as its l-nearest neighbor for I= == 1,...,.. k. According to Algorithm each z is selected with probability 1/k. Let & be a random variable with probability mass function 1/k for 1 1, ..., k. and is independent of both Y and Z. Given Z. we denote by X the sample generated by the k-nearest-neighhor local sampling mechanism, which follows the distribution p(z|Z)- p(z| )
P(c
x... ( The sketch proof of Theorem proceeds as follows. We first apply Pinsker's ss inequality, which relates the total variation distance between p(x 2) and f(x 2) to their KL divergence. We then establish a connection between the KL divergence and the discrepancy between Z and z(!) Finally, by Lemma | we derive Theorem 2 Theorem 2. Under Assumptions and 2, we have drv {p(x) Z), f(x Z)) = op(1) as n cc.
Remark 1. The constant | in Assumption I is used to establish an upper bound on the KL divergence between p(z) and p(æ ( which can be used to bound drv {p(x) (((al)), through Pinsker's inequality.
In Theorem we bound the excess type I error conditionally on Y and Z by the total variation distance between P(- [Z) and p(-1Z).
Theorem 3. Assume Hu X IL YIZ is true. Under Assumptions 1 and 2, for any signif- icance level a E (0, 1), the p-value obtained from Algorithm satisfies P(p S a a Z) s drv (p( T (()
The type error rate can unconditionally controlled based on Theorems 2 and which implies that P(p A a: Hn) I a + o1) as approaches infinity. Now we turn to analyze the power of our test in asymptotic scenarios in Theorem 4 Some intuitive explanations on both type error rate and
power can be found in the third paragraph of Section 3.2
7

---

Theorem 4. For any E (0, 1) and the number of repetitions B, satisfying B, + co, let p be the p-value calculated in Algorithm ] Under Assumptions 1 and and the Assumptions stated in the Supplementary Materials, P(p S H)H + as - co.
The condition Bn - oo is mild and has been made in 7 43y. The consistency of our test heavily relies on the consistency of the CMI estimator, which hinges on whether the joint density of (X, X,Y', Z) generated by the I-NN sampling (Algorithm I), denoted as o approximates p(z, z)p(ylz) (denoted as g) well, in terms of TV distance. B and ca_ in Assumptions 1 and 2 are used to bound (Tv(6.9).
Remark 2.
Remark 3. According to Theorem f our test can achieve consistency against all the alternatives stated in Hi In contrast, the method proposed by 143) only achieves consistency against subset of the alternatives in H1, and [30] does not provide any theoretical results on the testing power.
Empirical Results
In this section, we present a comparative evaluation of our proposed method against the state-of-the-art (SOTA) methods on synthetic datasets. Our code is publicly available at: https://github.cco Shaaai kenwitch/NNi. SCIT. We specifically compare our method with six commonly used competitive CI
testing methods: KCIT 47 LPC CCIT [4]). CMIknn (38), GCIT (5). and NNSCIT [30]. We set the number of repetitions B = 200 and the neighbor order k = 7 for our tests. The XGBoost classifier was used in all of our experiments. Further elaboration on the choice of is given in Figure 3 in the Supplementary Materials. We set the significance level to a 0.05 and report the type I error rate and the testing power under H₁ for all methods evaluated in our experiments. All the results are presented as an average over 200 independent trials. We provide additional simulation studies and real data analysis in the Supplementary Materials.
Scenario I: the post-nonlinear model. The datasets used in our experiments are generated using the post-nonlinear model similar to those in 47 16 5 40). Specifically, define YX,Y,2 under Ho and H1 as follows:
Ho X = f1(Z 0.5c.), Y = f2(7 (0.kg)))cg)).
H1 X = fi(Z 0.5c..) + 0.5c,, Y == =f2( + 0.56g) + 0.5cb,
where Z represents the sample mean of Z = (=1.... ). all 2/ in Z, Gz2 cy and c are i.i.d. samples generated from the standard Gaussian the standard Laplace distribution, and fi and fa are randomly sampled from the set {x, x,, x², tanh(x), (x)
We conduct a comparative analysis of various tests from two perspectives. The results are shown in Figure [ First, we the sample size at n 500 and vary the dimension of Z from 30 to 100. In both the Gaussian and Laplace cases, our test shows good and stable performance in controlling type I error and achieving high power under H1 as the dimension of Z increases. LPCIT and NNSCIT have satisfactory performance in controlling type error, but LPCIT loses power under Hi, when the dimension exceeds 70, and NNSCIT has inadequate power for all dimensions. Although GCIT, CCTT, and KCIT have adequate power under H1, they inflate the type I error in almost all scenarios. CMIkan shows weak performance on both type error and testing power Figure 6 in the Supplementary Materials contains additional results on the performance of various tests in low dimensions of Z ranging from 5 to 30, further demonstrating the superiority of our test. Second, we vary the sample size from 300 2000 while fixing the dimension of Z 80. Our test maintains good control of the type I error and achieves high power, with the power approaching 1 when the sample size exceeds 500. However, LPCIT, NNSCIT, and CMIknn lose power under H1, while GCIT, CCIT, and KCIT either always or sometimes fail to control the type error well.
In the Supplementary Materials, Figure 7 shows the timing performance of all methods for single test. Our test is found to be highly computationally efficient even when dealing with large sample sizes and high- nigh-dimensional conditioning sets. In contrast, CMIknn and CCIT for sample sizes exceeding 1000, and LPCIT for dimension of higher than 50 are impractical due to their prohibitively long running time.
Scenario II: the heavy tailed model. We compare the performance of our test with the SOTA tests under the heavy tailed error distributions, as described in [00]. The data (X, Y, Z) is generated
8

---

Gauss, 500
Gauss, dz- 80
0.35
0.30
0C11
0.30
Die lerror
0.20
H42CTT
0.20
0.15
Oun
0.15
0.10
0.05
0.05
0.00
0.00
30
40
50
60
70
80
90
100
300
500
1000
1500
2000
Dimension of
Sa Sample size
Gauss, samples-500 500
Gauss, dz-80
1.0
1.0
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
30
40
50
60
70
80
90
100
300
500
1000
1500
2000
Dimension of
Sample
Laplace, 500
Laplace, dz=80
0.45
0.45
CCT
0.35
K211
0.35
d8 Sror
0.30
0.30
0.25
0.20
Dun
0.20
0.15
00.55
0.10
0.05
0.00
0.00
30
40
50
60
70
80
90
100
300
500
1000
1500
2000
Dimension of
Sample
Laplace, samples=500 500
Laplace, dz=80
1.0
1.0
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.4
0.4
0.3
0.3
0.2
30
40
50
60
70
80
90
100
300
500
1000
1500
2000
Dimension of
Sample
Figure 1: Comparison of the type eтгor (lower is better) and power (higher is better) of our method with six SOTA methods on the post-nonlinear model under Gaussian Laplace distributions in Scenario I. Left: The results when varying the dimension of Z. Right: The results when varying the sample size.
according the following model:
Ho X = 7 + Ei, Y ==Z+e E2,
H1 X = E1 = Z +1 E2,
(10)
where Z is the sample mean of Z = (2ı,..,
zis are i.i.d. samples generated from the standard Gaussian distribution, and 6 and é2 are independently generated from the standard Cauchy distribution. We keep the sample size fixed at 500.

---

Cauchy, samples
Gauss, samples
0.40
1.0
COT
0.9
0.8
Ce
ANSOT
0.6
Ourt
0.5
Tyee
0.3
0.2
0.05
0.00
0.0
5
10
15
20
25
30
40
500
5
10
15
20
25
30
40
50
Dimension of
Dimension
Cauchy, samples=500
Gauss, samples=500
1.0
1.0
0.9
0.8
0.7
0.6
Power
0.7
0.4
0.6
0.2
5
10
15
20
25
30
40
50
5
10
15
20
25
30
40
50
Dimension of
Dimension
Figure 2: Comparison of the type error (lower is better) and power (higher is better) of our method with six SOTA methods Scenarios II and III. Left: The results in Scenario II. Right: The results in Scenario III.
the chain structure. We further the chain strueture Y + Z X in 30 to validate the proposed test. The detailed procedure for generating the simulated data can be found in the Supplementary Materials. The sample size is fixed
Scenario III:
500.
Figure 2 demonstrates that our test effectively controls the type error and obtains adequate power in both Scenarios II and III under various dimensions of Z.
6 Conclusion
In this paper, we introduce a novel method for conducting the conditional independence testing. Theoretical analysis shows that our asymptotically valid and consistent against all alternatives in H,. Extensive experiments both synthetic and real datasets demonstrate that our method consis- tently outperforms commonly used SOTA methods. Moreover, our approach maintains computational efficiency even when the sample size and/or dimension of the conditioning set increase. Our method has the potential to enhance the applicability of causal discovery to real-world problems, such as gene regulatory networks or complex social networks, and to facilitate the discovery of relationships and patterns in complex systems. Ethically, we believe that our rather fundamental work has minimal potential for misuse.
Acknowledgments
Dr. Ziqi Chen's work was partially supported by the National Key R&D Program of China 2021YFA1000100 and 2021YFA1000101) National Natural Science Foundation of China (NSFC) (12271167, 22331005 and 11871477), Natural Science Foundation of Shanghai (21ZR 1418800) and Basic Research Project of Shanghai Science and Technology Commission (22IC1400800), Dr. Christina Dan Wang's 's work was partially supported by the National Natural Science Foundation of China (NSFC) (12271363 and 11901395). We thank the anonymous reviewers for their helpful comments.
10

---

References
[1] Chunrong Ai, Li-Hsien Sun, Zheng Zhang, and Liping Zhu. Testing unconditional and condi- tional independence via mutual information. Journal of Econometrics, 2022.
[2] Mona Azadkia and Sourav Chatterjee. A simple measure of conditional dependence. The Annals of Statisties, 3000)3300--312,, 2021.
[3] Kunihiro Baba, Ritei Shibata, and Masaaki Sibuya. Partial correlation and conditional correlation measures of conditional independence. Australian & New Zealand Journal of Statistics, 46(4):657-664, 2004.
[4] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm. Mutual information neural estimation. In International Conference on Machine Learning, pages 531-540. PMLR, 2018.
[5] Alexis Bellot and Mihaela van der Schaar. Conditional independence testing using generative adversarial networks. Advances in Neural Information Processing Systems, 32, 2019.
[6] Wicher Pieter Bergsma. Testing conditional independence for continuous random variables. Eurandom, 2004.
[7] Thomas B Berrett and Richard J Samworth. Nonparametric independence testing via mutual information. Biometrika, 106(3):547-56) 2019.
[8] Thomas B Berrett, Yi Wang, Rina Foygel Barber, and Richard Samworth. The conditional permutation test for independence while controlling for confounders. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):175 197, 2020.
[9] Catherine Blake and Christopher J Merz. Uci repository of machine learning databases, 1998.
[10] Zhanrui Cai, Runze Li, and Yaowu Zhang. A distribution free conditional independence test with applications to causal discovery. Journal of Machine Learning Research, 23(55): 2022.
[11] Emmanuel Candès, Yingying Fan, Lucas Janson, and Jinchi Lv. Panning for gold: mmoeelxx knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3):551-577, 2018.
[12] 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.
[13] Thomas Cover and Peter Hart. Nearest neighbor pattern classification. IEEE transactions information theory, 1311:221-27, 1967.
[14] Thomas M Cover and Joy Thomas. Elements of Information Theory. John Wiley & Sons, 2012.
[15] Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans image synthesis. In Advances in Neural Information Processing Systems, 34, 2021.
[16] Gary Doran, Krikamol Muandet, Kun Zhang, and Bernhard Schölkopf. A permutation- based kernel conditional independence test. In Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence, pages 132-141, 2014.
[17] Sebastian Engelke and Adrien S Hitz. Graphical models for extremes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(4))871 1922, 2020.
[18] Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf. Kernel measures of conditional dependence. In Advances in Neural Information Processing Systems, 20, 2007.
[19] Weihao Gao, Sewoong Oh, and Pramod Viswanath. Breaking the bandwidth barrier: Geometri- cal adaptive entropy estimation. In Advances in Neural Information Processing Systems, 29, 2016.
11

---

[20] Weihao Gao, Sewoong Oh, and Pramod Viswanath. Demystifying fixed k-nearest neighbor information estimators. In Information Theory (ISIT), 2017 IEEE International Symposium on, pages 1267-1271,201
[21] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep learning. MIT Press, 2016.
[22] Guido W Imbens and Donald B Rubin. Causal inference in statistics, social, and biomedical sciences. Cambridge University Press, 2015.
[23] Kirthevasan Kandasamy, Akshay Krishnamurthy, Barnabas Poczos, Larry Wasserman, and James Robins. Nonparametric von mises estimators for entropies, divergences and mutual informations. In Advances in Neural Information Processing Systems, 28, 2015.
[24] Eugene Katsevich and Chiara Sabatti. Multilayer knockoff filter: Controlled variable selection at multiple resolutions. The Annals of Applied Statistics, 13(1):1-33, 2019.
[25] Immun Kim, Matey Neykov, Sivaraman Balakrishnan, and Larry Wasserman. Local permutation tests for conditional independence. The Annals of Statistics, 50(6))3388-344 2022.
[26] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT Press, 2009.
[27] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger. Estimating mutual information. Physical Review E, 6666066,,,,, 2004.
[28] Chun Li and Xiaodan Fan. nonparametric conditional independence tests for continuous variables. Wiley Interdisciplinary Reviews: Computational Statistics, 12(3):e 1489, 2020.
[29] Honghao Li, Vincent Cabeli, Nadir Sella, and Hervé Isambert. Constraint- causal structure learning with consistent separating sets. In Advances in Neural Information Processing Systems, 32, 2019.
[30] Shuai Li, Ziqi Chen, Hongtu Zhu, Christina Dan Wang, and Wang Wen. Nearest-neighbor sampling based conditional independence testing. arXiv preprint urXiv:2304.04160 04183, 2023.
[31] Dimitris Margaritis. Distribution-free learning of bayesian network structure in continuous domains. In Proceedings of the 20th National Conference Artificial Intelligence, volume pages 925-830, 2005.
[32] Arnab Mondal, Amnb Bhattacharjee, Sudipto Mukherjee, Himanshu Asnani, Sreeram Kannan, and AP Prathosh. C-mi-gan: Estimation of conditional mutual information using minmax formulation. In Proceedings of the 36th Conference Uncertainty in Artificial Intelligence, pages 849-88 858, 2020.
[33] Joris M Mooij and Tom Heskes. Cyclic causal discovery from continuous equilibrium data. In Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence, pages 431-439, 2013.
[34] Sudipto Mukherjee, Himanshu Asnani, and Sreeram Kannan. Cemi: Classifier based conditional mutual information estimation. In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence, pages 1083-1093, 2020.
[35] Matey Neykov, Sivaraman Balakrishnan, and Larry Wasserman. Minimax optimal conditional independence testing. The Annals of Statistics, 49(4) 2021.
[36] Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags. In Advances in Neural Information Processing Systems, 33, 2020.
[37] Judea Pearl. Causality. Cambridge University Press, 2009.
[38] Jakob Runge. Conditional independence testing based on nearest-neighhor estimator of conditional mutual information. In Proceedings of the 21th International Conference on Artificial Intelligence and Statistics, pages 938-947, 2018.
12

---

[39] Karen Sachs, Omar Perez. Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein signaling protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721):523-527 529, 2005.
[40] Meyer Scetbon, Laurent Meunier, and Yaniv Romano. An asymptotic test for conditional independence using analytic kernel embeddings. In International Conference on Machine Learning, pages
PMLR, 2022.
[41] Rajat Sen, Ananda Theertha Suresh, Karthikeyan Shanmugam, Alexandros G Dimakis, and Sanjay Shakkottai. Model-powered conditional independence test. In Advances in Neural Information Processing Systems, 30, 2017.
[42] Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. The Annals of Statistics, (8(3):1514-1553) 2020.
[43] Chengchun Shi, Tianlin Xu, Wicher Bergsma, and Lexin Li. Double generative adversar- ial networks for conditional independence testing. Journal of Machine Learning Research, 22(285):1-32, 2021.
[44] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT Press, 2000.
[45] Eric V Strobl, Kun Zhang, and Shyam Visweswaran. Approximate kernel-based conditional independence tests for fast non-parametric causal discovery. Journal of Causal Inference, 7(1), 2019.
[46] Xueqin Wang, Wenliang Pan, Wenhao Hu, Yuan Tian, and Heping Zhang. Conditional distance correlation. Journal of the American Statistical Association, 1106512) 1734, 2015.
[47] Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Proceedings of the 27th Conference on Uncertainty in Artificial Intelligence, pages 804-813, 2011.
[48] Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference Learning Representations, 2020.
13

---

Supplementary Materials for Neeeehbrstt- Local Sampling Based Conditional Independence Testing"
A Theoretical results
A.1 Proof of Lemma 1
Proof. Recall that is a random vector taking values in Euclidean space [R= (12) where d, is the dimension of Z and || ||2 is Euclidean distance. Z1, Zz. .... Zn are ii.. random vectors according to p(z). For a fixed x E IR we denote by zí"1 (z)...., z (z) a reordering of Z1, Zz,.... Zn
according to the increasing values of ||Z; 2112. that is,
1
=||2 S ... - 1 (2) 21|2-
Define the set G = {z E Rdd vo > 0, P([ : (ww) E Sz(0))) > 0}, where Sz(6) = for E Rd1 || - z1I2 < 8). For convenience, we omit w in probability in the following paper. For example,
write P(Z E (6)) instead of P([ Z(w) E (o))). By definition, for E G. V6 > 0,
P(Z E (6)) > 0. Let G be the complement of G. Then, for 2 E G., 3rz > 0, s.t. Vr < Tz+
P(Z (r)) O. Note that p( (z) z||2 > 6) P(z (k) (z) & (t))).
In order to prove (z) Zl12
0 a.s., it is sufficient to prove Vg >
lim, P(supm2n |/th (Z) (Z) Z > 5) O. We can obtain
P(sup (2/1/11(2))))) (Z) Z112 > 8) S P({supm 21 (Z) 2112 > 83 (Z
2) GJ) + P(ZEG)
1(s 2|12 > 0)p(z)dz + P(Z 3 6*)
JG
- Σ m2n { P(Z 69). (1I)
First, consider the first term of (11 We have
P(Z (2) & S.())) = P(Z (1) (z), z(2) (z),.... Z (m)) (2) $ S5)
m
+ P( (=) E Sz))), ( (z)...., Z(m) (z) $ (0))
P(Z (z), (2) (6). 7 (z) .))........................... Z(m) (z) & (o))
+ + ((
7 (2) E ( (z()), 2 (2), (z)...., 7(m) S.))
m
By setting P(Z Sz())) 7 we have
P( (z) & (f)) (1 7 - -7 m C 2 - -)
+ C k-1 - -) m-b
(12)
Consider the j-th term of 12 Let C1 := 70 /j1, C2 Cie' and C3 == Czel (1 73- By using Stirling's approximation, we have
m!
lim C T (1 -) lim m=000
7)
m- (m j]う!!
/2mm m m
lim C m-
22mmm - j)( m mj (1 7) m-
lim C2V (m - -j) (m 5) m- (1 - -) n-j
m
mm
m-000
lim Oze (m + (1)(1-1 7mm
lim Cze'm (1 (1 m-
m--
m-+00
lim Cym' (1 - mmm m=000
14

---

(12 (z) & Sz())) S (1m when m is large Thus, there exists C₁ > 0 such that P(Zm enough. holds that
Cam C-m 100-11111100 - 7) Cam++1(1 at+11( 7) +0, as m + +00.
m
(k) (2) & Sz()) o(m and Thus, for G and large enough, V8 > O, we have P P(Zm))
Σ P(Z (z) $ Sz())) S E m2n m
1
which shows 1im Σmzn PlZm (K) (2) (2)) for G. So by Lebesgue dominated
convergence theorem, we obtain lim Jo Emzn P(Zm) (z) & (z(0))((z)dz
Second, consider the second term of [1p. To prove that P(Z E G 0, aim to construct countable open cover of G and show that the probability of the random vector falling into each of these open balls is zero. By the property of Go, for every 2 E G there exists r= > 0 such that for all r < T=. P(Z E (r)) = 0. Furthermore, using the separability of Euclidean space and the density of the rational number set, we can approximate z using points from Qd with Q being the rational number set. Therefore, for every G there exist Q n (Tx) and Q 200
0% (2g 2r:
such that & E (r(r) S (r(r).. Because P(Z E (T())) = 0, we conclude that P(Z E S2(r)) 0. Define
F F= (r) 3z G", such that = E Sr(r)) S (r(rr) with a Q and r E =Q).
Note that the elements in set F are mutually distinct. By the construction of S,(r), F forms countable open cover of G", and the probability of Z falling into each open ball in F zero. Using the monotonicity and countable additivity properties of probability, we have P(Z E G*) < P(Z Sr) S (n) E (r) = O. Thus, we conclude that P(Z E GO) = O.
We therefore conclude that, V8 0, lim, P(sup molupmznllaa 11 (Z) - Z1I2 8) 0. This finish the proof.
A.2 Proof of Theorem
Proof. By Pinsker's inequality, have
Drv s (1)
Note that T(E= = + TfE = k) 1. By the definition of p(x) Z), we obtain
p(x Z)
DrLtp(a|/2,P(a1l1) = plz|Z) log
dr
p(x) -11==1111111111
p(z)
= | plzlZ) log
p(z] Z)) ((E-1)
1=1 |(ar| Z ) dr
p(x Z) dr p(x z(²)
Σ T/E 11 pxzl log
1=1
Il5 = (xx==)])xx(((
Then, by Taylor's expansion, we have
DkL10(=12)11(<122")
मे a Dxxtpcaa22111()t
=Dri(p(s +
(2 Z)
1==7
+ + (2 -
1
8²
D Dx.kolaa)))pplallllll)
(2(0) (2(0-2)) -2),
0'020
15

---

where = + 2) with 0 S 1.
Note that DEL (2x2))2)222 [ log p(a] 2) dr 0. By Lemma and Assumptions p(u]2)
and 2, we have
0 Oz kk 1xx.ooal22)llll
-= =- f p(z)
a
log p(z(22)
dx
0z
x'=2
Oz / p((ol')do
0
-0
1===
and
82 Oz D Dwl(p(a D =- =-f -fr(t)))
0²
0z' log logp(z|)')
dr Lo(Z).
d'=a
z'=0
This means
2 (2 - Lz) -Z)
Note that - + Z a.s. implies that - converges to Z in probability. Then 4A >> 0, for defined
in Assumption we have
- 2112 S {
P(IZ.." - ())
- P $11Z - 21 > 6) + +(11Z - Z112 >
o(1)
Thus, DEL (u(x) (xL(u(x\Z)\p(</// = (1).
Because 115 = 1} - 1 for = 1, ..., k, and is is finite, we obtain
[ = (S==19xx.oo<<(/3)bba21277 (
=
Finally, we conclude that drv (tt) ( (1))
A.3 Proof of Theorem
To prove Theorem the following lemmas needed.
Lemma 5.
Let A be drawn from Z(-1Z), independently Yf Y. x ( are i.i.d. samples
drawn from the neiibost-neihhor local sampling mechanism based (X, Y. Z). For any statistic T, the + statistics
T Z), T(X N,,,Y,Y V)Z).... (B) Y
are exchangeable conditionally on Y and Z.
Proof. We have that the B+1 triples (X, Y, Z), (1) (1) Y, Z).... (.....,((())) (B) Y, T) are i.i.d. samples
drawn from the same mechanism after conditionally on X Y and Z, where Xo is the order statistic of X. Note that. T is measurable. Thus, T(X, Y. Z). T(X (1) Y. Z). ( Y. Z) are i.i.d. after conditionally on X0, and Z. Conditionally on Xo: Y and Z, denote their cumulative conditional distribution F( X e Y. Z). Denote X as X (0), Then, for any .o. R and
16

---

any permutation π (TT()). T(m) of the indices {0, (0.1..... B), we have
P(T(X (0) Y, 2) lo, T(X Y, 2) S 41 ,..., T Y 2) S LB 2)
2 (P Y, 2) - to T(X (1) 2) S (..... T(X (B) ,Y, S (), Y, 2]}
=Ex,IX x (PTT (0) Y, 2) s tolx0)
P...,.,,(( (B) Y, Z) s <t (.2))
B
- =Explr.z
II X (1 Y. {
=Ex,,x,2 ( ((o) Y. 2) S to..... (( 2, S ttl (t Y, z))
(T())) = =P(T(X m()) Y, 2) 2 to. (
Y, Y 2) - tg Y Z).
Thus, the desired result follows.
Let denotes equality in distribution. We present the following Lemma:
Lemma 6. For any two bi-tuples (U, V) and (U', V'), if Vu, (VU = =u) = (V'IU" =u), we have drv {(U, V), (U' V") = dTv(U.U)
Proof. Denote the joint density functions of (U, V) and (U', V`) by PU.
pu.v (u, v) and pu", v (u', v').
respectively. According to the equivalent definition of the TV distance, we obtain
1 // 2 (pr.v (1, duudo
dTv (UU (U, VO)
2 (DVIo(olu)po(u)) - pv (w) (o) (u)
1
//
1
//
2 ( - (ur(()))d duudv
1
/ // (oler Ipo(u) pur du
2
1
2 /
(w)))du
dTv(U,U)
Now we present the proof of Theorem
(B) be i.i.d. drawn from the d-mearest-neighbor local sampling mechanism, see Algorithm 3 Now let X be an additional sample drawn from P(- Z) independently of Y. Let (1)
Proof. Let (1))
X (B) be i.i.d. drawn from the k-nearest-neighho local sampling mechanism after we observe X instead of X. Because (L)) ..... (B)), conditionally on X, Y and Z, is generated from the same mechanism (1(1), X(B),, conditionally X, Y and Z, for all æ E R", have
(1)
X x = x; Y, 2) (1)(1)
X x x. Y, 2).
Then, by applying Lemma b we obtain
drv drv[(XxXX0) (1) Y. V)Z), (X, x()) x 2i}
-rv( Z), ( Z)) = drv (p( ()))(
Define x := { (z. ()) x (1))
|
2 1(T(x (b) Y. Z) > T(x, Y, 2)} (1+
where 1(.) is the indicator function. Note that in case, the statistic T is selected to be CML. Then,
17

---

it follows that
P(p S a|Y, P((
X(8)) Ex Xa Y .2 X E Xa Y. Z) + P((XX ,..., X(O1 E 2))
= =P(XX, X(1)
- P((X, x(1)
(1)) E Xa YZ))
P( X
x(11) E EXol (Y,2)
ddvvl(X, X(",,
X X Y. (X, X(1,, (0,,,,
X
P((X,
( E Xa Y. 2) + drv (2(1))) Z)).
Applying Lemma 5 and the property of rank test, we obtain P((X (1), (1)
X (B)) E Xa Y Z) S
a. Finally, we have P(p S YYY, Z) S drv (p(-)))(,,))). p(-
Because the TV distance is bounded by 1, marginalizing the above inequality over Y and and applying Theorem 2 and Lebesgue dominated convergence theorem lead to
P(p S aoo) s
())))) o())
A.4 Proof of Theorem
Here present the assumptions given in [34] that ensure the consistency of CMI estimator. We denote the point (z, y, z) as w E IRd x Rd, x Rd, Let f(w)) = p(a, y, x) be the joint density function of (X, Y, Z), g(w) = 2) be the joint density function of (X. Y. Z) under Ho. and o(w) = 0(x, y. 2) be the joint density of (X,Y² Z)) produced by Algorithm For a given point w, ((w) = P( = 1(w) represents the predicted positive label probability generated by a classifier. When the prediction is from a classifier with parameter 0, it is denoted as (w(w). ). For convenience, the parameter 0 in 7 is dropped when it is understood from the context. According
Algorithm 2. P( 1) = P(1 = 0) = 1/2. Define the population binary-cross entropy loss over the joint distribution of data and label BCE(7) - (Fwi (1 log (WW) + (1 1) log(1 (W))).
Assumption (A1): f()) and (6-) admit densities in compact subset W C IR R& x Rd.
Assumption (A2): For some constant a, 3 0, a S fww), o(w) S 8, Vw.
Assumption (A3): Clip predictions such that (ww) E [T, - r]. Vw, with < S a/(a + 3).
Assumption (A4): The classifier class Co is parametrized by 0 within a compact domain Q Rt. There exists a constant such that ||0||2 S K, and the classifier's 's output is I-Lipschitz with respect to A.
Assumption (A5): f((2)) S C.. Vd 2 2, where Cg a constant.
We denote the CMI estimator CMI based on Algorithm
2 by Dig (fll)). The true CMI of (X,Y, Z)
CMI T(X; YlZ) D (fllg). Then we have the following Lemma:
Lemma 7. Under Assumptions I and 2 and (AI)-A5), we have DRL((110) P. DKL(fllg)
Proof. By the definition of convergence in probability, it is sufficient to prove vS 0, Vn 0, 3N,
when n > N, P(IDRL(S1))) (fllg) > 5) <n
Note that
P(IDR - = = - (kL(fll0) + DKI(fll)) DKI(fllg) > 6).
Applying Theorem 1 in 34), we have Dk)) -xx - 0, which means for 8/2 > 0 and 17 0, 3N1 when n No. (() KL (/1lo) (KI()]l0)
(n)
<.
Now consider the term ( - DKL( Let be the point-wise minimizer of binary-cross entropy loss based on fw) and ((w), and 7 be the point-wiss, minimizer of binary-cross entropy
18

---

loss based on f(w) and ((w). Applying Lemma in [34.. we have have>'(w))) (1w)/11-1'(00 = f(w)/(())
and [ f(w)/g( (w)). Then, by the definition of KL divergence, it follows that
f(w) Ō(w)
Bs(w) log f(w)
( - Dxl((1l))
Brw log
g(w)
7((w) -log (w)
B&
log
- 7 (ww)
(w)
1 y (w) loog
7" (w)
Es(w) log
7 (ww)
(w)
1 Em(w)
['ww 1- 7" (w)
((w)
(w)
(b(w) g(wp)
Bl(w) T
f(w)
1 1-T /I/
|o(I, v, g( y, 2) ddddddd
2(1 7)
drv(6.9).
The second inequality follows from Lagrange's mean value theorem and Assumption (A3).
Applying Theorem in [1 Ve1 S with e being defined in Assumption 1, we have drv (6.g) S
b(n), where
b(n) V 4 (nrja,) ) )/dd= + BeIG(2ca, + + exp - nid Cd, c (:t2) G(eca, f
3 C,2l/d: r(1/dd)
Here, 39 is defined in Assumption 1, Cs is defined in Assumption (A5), d. is the dimension of Z, I'()) is the gamma function, Yd. is the volume of the unit radius l2 ball in Rds, ca, is defined in Assumption 2, and V8 > 0, G(Ø) P(p(Z)+ I d).
Because E1 can be arbitrary small, we conelude that 1om b(n) 0. So arrive at lim DxL(fll)) - (fllg) = 0, which means for 8/2 0. NN2, when > No. [Dxz(fll) (fllg) < 6/2.
Then for 8 > 0 and v > 0, take N = max max(N,, N2), when n > N.
P(IDR20111 - +
Dxi((ll)) - DKL(/llg) 5) P(IDRL DRI(SIlO) (Dk((fl0) Dkr(fllo) 5)
< P(ID)) (fI1O) - DRI(fll0) 6/2) <n
holds. This finish the proof.
We therefore conclude that CMI is a consistent estimator of CMI. When considering CMI based on the sample (x()) Y, 2) (b 1,..., B) drawn from the k-nearest-neighhoor local sampling mechanism as depicted in Algorithm 3 we can state: Under Assumptions 1, 2, (A1)-(A2) with
(w) and ((w) replaced by densities of the distribution of (X, Z) and the corresponding 1-NN
distribution, respectively, and Assumptions (A3)-(A5), vb 1,..., B, CMI (b)
is consistent
estimator of (*) where (M) [(X Y|Z). Now let's present the proof of Theorem 4
19

---

Proof. Write 1)) as PH. (-). By Markov inequality, it follows that
1+ R.
=== 1(CMI
CMI)
Ph. (p a===== = PH. PH
1 Bn
(b)
EH. [1+
1(CMI
CMI
a(1 Bn)
0=1
Bn
PHr (CMI* (1)
2 CMI)
a(1 +B))
a(1
Ba)
1
1
(1)
+ -Hi (CMI
2 CMI).
a(1 + B,)
a
Because ((1 11 Z12, CMI ((
O. Then, vo >
PHY (CMI (1)
2 CMI) - Pm. ((CMI S CMI (1)
n {|CMI -
CMI s $0]
+ PH, (ICMI (1) ( >}) 6)
< Pn. (CMI S 5) + Prr, (ICMI (1)
CMI > 8).
Next, we have
PR: (CMI S 8) S Pn, ([CMI 5 n {|CMI - CMI| $ + PH. (ICMI CMI > 6)
S Pn, (CMI- 8 < CMI S 8) + PH, (ICMI CMI| 8).
Thus, we conclude that
PH. S aa)
a(1 + Br) PH,, (CMI 8 CMI S 8)
+ PH, (ICMI CMI| > 5) + PH, (ICMI (1) CMI > 8)].
Under H1, CMI > 0. Take 8 = CMI/4> > 0, we obtain
lim Pr. (p a) 1.
B Additional Empirical Results
B.1 The choice of the neighbor order k
To investigate the impact of the parameter is on our proposed approach, we employ linear uniform model. To accomplish this, we generate synthetic data in the following manner:
Ho X Eg. Y Ey. and Z
(aiform(-1,115
x = Ens Y = ax + 0.5eg, and Z Uniform(-1,1),
(13)
where e, and Ey are generated independently from the uniform distribution over the interval [-1, 1]. The parameter is randomly generated within the range of [0. 2]. As is shown in Figure method achieves effective control of type I error and exhibits the highest power under H₃ across all dimensions when A 7. Therefore, consistently set A in all experiments.
B.2 Empirical results for Scenario
We demonstrate the effectiveness of our approach and compare it with alternative methods in Scenario (13). The results are shown in Figure 4 which pertains to high-dimens -ine-niinnii Z, and Figure 5 which tocuses on low-dimensional Z. The results consistently demonstrate that our test achieves favorable performance in terms of the type error and power under H1. Although LPCIT, CMIknn, and NNSCIT effectively control the type error, they exhibit noticeably lower power compared to our method, when the dimension exceeds 60, often by a substantial margin. Furthermore, KCIT, GCIT, and CCTT all yield high power under H,, but they either always or sometimes suffer from inflated type errors.
20

---

0.10
Uniform, 5amples=500
1.00
Uniform, samples- 500
6412
1ype entor
0.00
0.04
0.93
0.00
100
0.85
30
50
Dimension of
Dimension of
Figure 3: Comparison of the type I error (lower is better) and power under Hi (higher is better) for our test in Scenario 13)
across different values of k.
Uniform, samples=500
Uniform, dcm80
0.35
0.30
0.25
0.25
LP3T
0.20
0.20
0.15
0.15
0.10
0.10
0.05
0.00
40
70
300
1000
1500
Dimension
Sample
Uniform, samples-500
Uniform dz-80
1.0
0.9
Powel
0.8
0.6
0.6
0.4
500
Dimension
Sample size
Figure 4: Comparison the type I error (lower is better) and power under H1 (higher is better) of our test with six SOTA tests in Scenario (13). Left: The results when varying the dimension of Z. Right: The results when varying the sample size.
B.3 Additional empirical results for Scenario
In Figure we present the type error and power under Hi in low dimensions of Z ranging from 5 to 30 for Scenario (Eq. (9) with Gaussian Laplace noises. It can be observed that our test and LPCIT consistently achieve good and stable performance in terms of type error and power under Hj, when the dimensionality of Z is lower than 30. On the other hand, GCIT, CCIT, and KCIT exhibit high power under H₃ but fail to control the type error. NNSCIT and CMIknn demonstrate relatively good control of type I errors but lack sufficient power under Hi.
B.4 Computational efficiency analysis
Figure shows the timing performance all methods for a single test under Scenario with Laplace noises. Our test is found to be highly computationally efficient even when dealing with large
21

---

samples-500
Unifortn, samples-500
V01
0.3
Tipe Lener
0.2
0.0
10
15
20
25
0.8
5
10
15
20
25
Dimension of
ot Z
Figure 5: Comparison the type error (lower is better) and power under H1 (higher is better) of our test with six SOTA tests in Scenario (13).
5a0ples=500
Laplace, samples-500
0.7
CCIT
0.6
OMirre
0.3
0.2
0.0
0.0
15
20
30
10
15
20
Dimension of Z
Dimension of Z
Gauss, samplcs-500
Laplace, somplos-500
1.0
1.0
0.9
b 0.8
0.6
0.6
0.5
10
15
20
23
30
5
10
20 Dimonsion of
15
25
30
Dimansion of
Figure 6: Comparison of the type error (lower is better) and power under Hi (higher is better) of our method with six SOTA methods on the post-nonlinear model under Gaussian or Laplace distributions in Scenario Left: The results under Gaussian distribution. Right: The results under Laplace distribution.
sample sizes and nigh-dimensional conditioning sets. contrast, CMIknn and CCIT for sample sizes exceeding 1000, and LPCIT for dimension of Z higher than 50 are impractical due to their prohibitively long running time.
B.5
The detailed experimental setup for Scenario III
For the chain structure Y + Z > X, we generate synthetic data as follows:
Ho Y N(L,)), 1), Z = Ya E X = 6T6 E2,
H Y ~ N(,, 1), = Yater, X = 2T6 + Y E2
22

---

samples=500
dz=80
300
330
2- onT
270
CMirn
UndI
240
210
210
180
Run time
180
Run
150
120
120
90
30
30
100
300
500
1000
1500
2000
Dimension
Sample
Figure 7: Running times in seconds function of sample size or dimension of Z on the post- nonlinear model under Laplace distribution Scenario Left: The results when varying the dimension of Z. Right: The results when varying the sample size.
where and b are both dimeennsii, the entries of and b are both randomly and uniformly sampled from [0, 0.3). is generated from dz-dimensional standard multivariate Gaussian distribution, and (2 is sampled from standard Gaussian distribution.
C Real Data Analysis
In order showcase the superior performance of test, we conduct a comparative evaluation against other state-of-the-art (SOTA) CI tests using real datasets. We assess the effectiveness of our method along with six SOTA approaches specific datasets: the ABALONE dataset [9] and the
Flow-Cytometry dataset 39
C.1 Real ABALONE dataset
The ABALONE dataset comprises measurements obtained from a study conducted to predict the age of abalones based on their physical characteristics. The dataset is pub- licly available the UCI Machine Learning Repository and can be downloaded from https//aarchiveiis..cem
In our evaluation, we consider the graph structure recovered by [31 the ground truth, as depicted in Figure of their paper. This graph represents the causal relationships among the 8 variables in the dataset. We specifically select 35 CI relations and non-CI relations from this graph. The philosophy used is that a node X independent of all other nodes Y in the graph when conditioned on its parents, children, and parents of children [26] 41). Additionally, if there exists a direct edge between node and node in the graph, they are never conditionally independent given any other set of variables. As result, the conditioning set Z can be arbitrarily selected from the remaining nodes. The dataset consists 4177 samples, and d, varies from to 6.
In order evaluate the performance various tests, we utilize precision, recall, and F-score as evaluation metrics. Precision is calculated as TP((TP+FP), where TP represents the number of true CI instances correctly identified, FP represents number of non-CI instances incorrectly identified as CI. Recall is calculated TP/(TP+FN), where FN represents the number of CI instances not identified. The F-score then computed as the harmonic mean of precision and recall, given by x precision recall / (precision + recall) 22). TN represents the number of correctly identified true
23

---

non-CI instances. Table ! presents the results for all methods. It should be noted that we do not record the results for GCIT does correctly identify any CI relations. Our approach successfully identifies 31 CI relations and 32 non-CI relations, achieving the highest F-score among the testing methods, while maintaining high precision and recall.
Table 1: The TP, TN, precision (pre), recall (rec) and F-score of test and six SOTA methods for the real ABALONE dataset.
Method
TP
TN
Pre
Rec
F-score
KCIT
35
0.1429
0.2501
CCTT
12
34
0.9231
0.3429
0.5
CMIknn
22
35
0.6286
0.7720
LPCIT
35
0.1429
0.2501
NNSCIT
33
6
0.5323
0.9429
0.6805
Ours 31
32
0.9118
0.8857
0.8986
C.2 Real Flow-Cytometry dataset
The Plow-Cytomedd dataset is a widely used benchmark in the field of causal structure learning
This dataset captures the expression levels of proteins and phospholipids in human cells [39]. The data can be obtained from the website https://www.ccieet wwt.nnce..r///o//..... 1.. 126/science. 105809. In our evaluation, we consider the consensus graph proposed in [33 as the ground truth, which has also been adopted by 41 for verifying CI relations. Figure 5(a) in [33 illustrates the causal relationships among the 11 proteins in the dataset. Following the philosophy outlined in Section C.T we select 50 CI relations and 40 non-Cl CI relations from this graph. The number of samples is 1755 and d, varies from to
36 48..
Table 2 presents the results for all tests. Our method outperforms other approaches by correctly adentifying 47 CI relations and achieving the highest recall and F-score.
Table 2: The TP, TN, precision (pre), recall (rec) and F-score of our test and six SOTA methods for the real Flow-Cytometry dataset.
Method
TP
TN
Pre
Rec
F-score
KCTT
32
30
0.7619
0.64
0.6957
CCrT
33
29
0.75
0.66
0.7021
CMIknn
41
26
0.7455
0.82
0.7810
GCIT
40
24
0.7143
0.8
0.7547
LPCIT
38
25
0.7170
0.76
0.7379
NNSCIT
33
26
0.7021
0.66
0.6804
Ours
47
0.7344
0.94
0.8246
24


Reviews:
[{"review_id": "4gfqjkb5Cb", "reviewer": "Reviewer_BxQV", "summary": "The paper proposes a novel method to test conditional independence (CI) by estimating conditional mutual information (CMI). A kNN local sampling strategy is employed to empirically sample from an unknown conditional distribution. The method is proved to have type-I error controlled, and power converging to 1.\n\nThe main contribution of the novel method consists of three points. First, it addresses challenges posed by high-dimensional conditioning variables and limited data samples. Second, it operates without assumptions about distribution forms or feature dependencies. Third, it avoids dataset splitting. Empirically, it is computationally efficient.", "questions": "(1) The paper claims novelty in avoiding sample splitting. It seems to me that Algorithm 1 & 2 are still splitting data into $V_1$ and $V_2$. Does this splitting not weaken test power?\n\n(2) More runtime analysis would be appreciated. The paper's claim on enhanced computational complexity sounds legit, but requires more evidence to support.\n\n(3) How does $\\beta$ and $c_{d_Z}$ in Assumptions 1 and 2 affect the convergence rates in Theorems 2 and 4?\n\n(4) Line 180: \"Furthermore, it intuitively suggests that our test can achieve high power under $H_1$\". This claim was based on $\\hat{\\text{CMI}}$ should be positive whp. This quantity, however, is replaced by $\\tilde{\\text{CMI}}^{(b)}$ in practice. It the latter still positive whp?", "limitations": "The paper claims that their novel method operates without assumptions about distribution forms or feature dependencies. However, the two assumptions on smoothness do not sound mild.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "(1) Good originality. Using kNN as a sampling method to circumvent prior knowledge on conditional distribution $p_{X|Z}(x|z)$ is a novel idea.\n\n(2) The algorithm is a high-quality one, in terms of enhanced computational efficiency.\n\n(3) The simulation settings is a nice touch to showcase that the novel method operates without assumptions about distribution forms.", "weaknesses": "(1) The paper claims an improvement upon the current literature under high-dimensional settings. But this contribution is not reflected in either theoretical results or simulations. Specifically, in simulation the sample size is 500, while dimension of $Z$ is at most 100. The results are still in low-dimensional regime.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the helpful review!\n\n**Weakness**: Please refer to the global response for details regarding the claim “ high-dimensional conditioning variables\".  Further, we conduct additional experiments to assess the performance of our tests in higher-dimensional scenarios. We generate simulation dataset based on `Scenario I: The post-nonlinear model'  in the paper. We keep the distribution as the standard Gaussian, maintain the sample size of $n=500$, and set the dimension of $Z$ to be 200 and 500. We present the type I error rate and testing power in Table 1 in the pdf file attached in the global response, where the results are averaged over 100 independent trials. It is shown that our test maintains good control of the type I error and achieves high testing power. Notably, LPCIT is  very time-consuming when the dimension of $Z$ is high, we thus do not report the results of this test. \n\n**Question 1**.  Please refer to the global response for details regarding the claim “our method avoids sample splitting\".\n\n**Question 2**. First, we provide the total runtime and the time taken for a single test run for all methods on the real ABALONE dataset ($n=4177$) and the Flow-Cytometry dataset ($n=1755$).  Table 2 in the pdf file attached in the global response shows that our test maintains a high level of computational efficiency even when handling large sample sizes (e.g., $n=4177$). Second, please refer to Lines 152-158 in the paper for more details. Third, we find that our method has a computational complexity of $O(B\\cdot dz \\cdot n\\log n)$. Further comparisons of computational complexity with other methods will be conducted in our future research.\n\n**Question 3**. In Theorem 2, $\\beta$ is employed to establish an upper bound on the KL divergence between $p(x|Z)$ and $p(x|Z\\_n^{(l)})$: $$D_{KL}\\\\{p(x|Z)||p(x|Z\\_{n}^{(l)})\\\\}=\\frac{1}{2}(Z_{n}^{(l)}-Z)^{T}I_{a}(Z)(Z_{n}^{(l)}-Z),$$\nwhere $0\\leq \\lambda_{\\max}(I_{a}(z))\\leq \\beta$. Therefore, given that $\\beta$ is finite, the convergence rate of $d\\_{TV}\\\\{p(x|Z),\\widehat{p}(x|Z)\\\\}$ is determined solely by the convergence rate of $\\Vert Z\\_{n}^{(l)}-Z\\Vert\\_{2}^{2}$ and remains unaffected by $\\beta$, as indicated by the inequalities:\n\\begin{align*}\n d\\_{TV}\\\\{p(x|Z),\\widehat{p}(x|Z)\\\\}  \\leq \\sqrt{\\frac{D_{KL}\\\\{p(x|Z),\\widehat{p}(x|Z)\\\\}}{2}},\\\\\nD\\_{KL}\\\\{p(x|Z),\\widehat{p}(x|Z)\\\\} = \\sum\\_{l=1}^{k} I\\\\{\\xi = l\\\\} D\\_{KL}\\\\{p(x|Z)||p(x|Z\\_{n}^{(l)})\\\\}.\n\\end{align*}\n\nIn the context of power analysis as presented in Theorem 4, it is noteworthy that most CIT methods lack theoretical results pertaining to testing power. Examples include works by Runge et al. (2018), Berrett et al. (2020), and Li et al. (2023). Even when theoretical analyses of testing power are conducted, they never give the convergence rate of testing power. For instance, the method proposed by Shi et al. (2021) achieves consistency only against a subset of alternatives in $H_1$. Wang et al. (2022) specifically consider the high-dimensional linear model and analyze the asymptotic power under local alternatives.\n\nIn our work, we achieve the consistency of our test against all alternatives stated in $H_1$ when analyzing testing power. This consistency heavily relies on the consistency of the CMI estimator. The CMI estimator's consistency hinges on whether the joint density of $(X,Y',Z)$ generated by the 1-NN sampling (Algorithm 1), denoted as $\\phi$, approximates $p(x,z)p(y|z)$ (denoted as $g$) well, in terms of TV distance. To bound the TV distance between $\\phi$ and $g$, we utilize $\\beta$ and $c\\_{d_z}$. Precisely, $d_{TV}(\\phi,g)\\leq b(n)$, where\n$$b(n)=\\frac{1}{2}\\sqrt{\\frac{\\beta}{4} \\frac{C_5 2^{1/d_z} \\Gamma (1/d_z)}{(n\\gamma_{d_z})^{1/d_z}d_z}+\\frac{\\beta \\epsilon_1 G(2c_{d_z}\\epsilon_1^2)}{4}}+\\exp \\bigg{(}-\\frac{1}{2}n\\gamma\\_{d_z}c\\_{d_z}\\epsilon_1^{d_z+2}\\bigg{)}+G(2c\\_{d_z}\\epsilon_1^2),$$ \nwhere $\\epsilon_1$ is small enough, $\\Gamma(\\cdot)$ is the gamma function, and $\\gamma_{d}$ is the volume of the unit radius $l_2$ ball in $\\mathbb{R}^d$. As indicated in Theroem 1 of Sen et al. (2017), as long as $\\beta$ and $c\\_{d_z}$ are finite (Gao et al., 2017),  $d\\_{TV}(\\phi,g)$ can be small enough when $n$ is large enough. This ensures the consistency of our test against all alternatives stated in $H_1$.\n\n[1] Wang W., Janson, L. A high-dimensional power analysis of the conditional randomization test and knockoffs, Biometrika, 2022.\n\n**Question 4**.  We have proved that $\\widehat{\\mbox{CMI}}$ is a consistent estimator of $\\mbox{CMI}$  based on the original data set $(\\pmb{X},\\pmb{Y},\\pmb{Z})$. Because the  $\\mbox{CMI}$ is greater than zero under $H_1$, $\\widehat{\\mbox{CMI}}$ is positive whp. \n\nWhen the distribution of $X|Z$ is known, we can  draw  pseudo samples $\\pmb{X}^{(b)}$. In this case,  the estimator $\\widetilde{\\mbox{CMI}}^{(b)}$ based on $(\\pmb{X}^{(b)},\\pmb{Y},\\pmb{Z})$ converges to zero, but it is not necessarily positive. As a result, the $p$-value calculated using $p:=\\\\{1+\\sum_{b=1}^{B}{1}(\\widetilde{\\mbox{CMI}}{}^{(b)}\\geq \\widehat{\\mbox{CMI}})\\\\}/(1+B)$ is highly likely to be very small, which indicates the consistency of our test against all alternatives stated in $H_1$.\n\nIn practice, the distribution of $X|Z$ is  unknown.  We use the $k$-NN local sampling strategy (Algorithm 3) to draw pseudo samples $\\widetilde{\\pmb{X}}^{(b)}$. In the proof of Theorem 4, we have demonstrated that $\\widehat{\\mbox{CMI}}^{(b)}$ calculated based on $(\\widetilde{\\pmb{X}}^{(b)},\\pmb{Y},\\pmb{Z})$ converges to zero, though it is not necessarily a positive value. Consequently, the $p$-value calculated using $p:=\\\\{1+\\sum_{b=1}^{B}{1}(\\widehat{\\mbox{CMI}}{}^{(b)}\\geq \\widehat{\\mbox{CMI}})\\\\}/(1+B)$ is very small with high probability, indicating our test can achieve high power under $H_1$.  \n\n**Limitation**: Please see the global response for details regarding  the claim “our method operates without assumptions about distribution forms or feature dependencies\". "}, {"review_id": "m4Ql3t3uHT", "reviewer": "Reviewer_xTib", "summary": "In this paper, the authors propose a novel CMI estimator with a classifier-based approach. The estimated CMI is then used for conditional independence testing. This approach achieves asymptotic type I & II control and does not require prior assumptions on the distribution, types of correlations, or knowledge of the null distribution. Empirical results show promising performance. ", "questions": "1. Why is the CMI estimated using classifiers rather than directly optimising (4) like f-GAN? Why not directly use the maximum of (4) to approximate the CMI? This way, you would not need to split the dataset into training and testing sets right? ", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper is well-written. \n2. The proposed CMI is very natural for testing CI. The algorithm is easy to implement given any classifier.  \n3.  The proposed method is backed up by theory (Theorem 3 & 4). \n", "weaknesses": "1.  As authors have pointed out, prior works for estimating MI (e.g. Suzuki et al., 2008 and [3]) using density ratio do exist, and a conditional version of such an estimator is not exactly a huge leap forward (but no small feat either)\n\n2. The control of type I & II is asymptotic, and it is unclear how useful they are when the sample size is small. The claim of Theorem 3 and 4 is unsurprising (due to the consistency of KNN and smoothness of conditional distribution). If the proof is non-trivial, please point out.  \n\nSuzuki et al., 2008 Approximating Mutual Information by Maximum Likelihood Density Ratio Estimation", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the insightful review!\n\n**Weakness 1**. The Maximum Likelihood Mutual Information (MLMI) method proposed by Suzuki et al. (2008) directly models the density ratio and avoids density estimation. MLMI estimates the density ratio using maximum likelihood and formulates it as a convex optimization problem using Legendre-Fenchel duality. We may extend this to the conditional version. Nevertheless, the MLMI-based method heavily depends on selecting appropriate basis functions that can effectively capture the information in $Z$, and it can face challenges due to the curse of dimensionality (Mukherjee et al., 2020). In contrast, our approach employs the classifier-based estimation of the likelihood ratio, offering a solution capable of addressing high-dimensional challenges. \n\n**Weakness 2**. First, in the simulation study, we present  some results when the sample size is $n=300$ as can be seen in Figures 1 and 4. Even when the sample size is small, our proposed test  controls the type I error well and achieves adequate testing power. Note that, GCIT and CMIknn consider a minimum sample size of 500, and DGCIT, NNSCIT and CCIT all consider a minimum sample size of 1000.\n\nSecond, we discuss the non-trivial points in the type-I error control of our method in Theorem 3. If the conditional distribution of $X|Z$ is known, the p-value is calculated by $$p:=\\frac{1+\\sum_{b=1}^{B}{1}(\\widetilde{\\mbox{CMI}}{}^{(b)}\\geq \\widehat{\\mbox{CMI}})} {1+B},$$ as can be seen in Equation (8) in the paper. Because of the exchangeability of $B+1$ triples $(\\pmb{X}, \\pmb{Y}, \\pmb{Z}),(\\pmb{X}^{(1)}, \\pmb{Y}, \\pmb{Z}),\\cdots,(\\pmb{X}^{(B)}, \\pmb{Y}, \\pmb{Z})$ under $H_0$, this p-value is valid and  $P(p\\leq \\alpha |H_{0})\\leq \\alpha$ holds for any given $\\alpha \\in (0,1)$.  However, in our case, $X|Z$ is unknown. $\\widetilde{\\pmb{X}}$ is drawn from $\\widehat{p}(\\cdot|\\pmb{Z})$ instead of $p(\\cdot|\\pmb{Z})$. Therefore, the existing argument based on the assumption of a known conditional distribution $X|Z$ does not apply in this context. So, in the proof of Theorem 3, we  introduce  an additional sample $\\acute{\\pmb{X}}$ drawn from $\\widehat{p}(\\cdot|\\pmb{Z})$ independently of $\\pmb{Y}$. We proceed to independently and identically draw $\\acute{\\pmb{X}}^{(1)}, \\cdots, \\acute{\\pmb{X}}^{(B)}$ using the $k$-nearest neighbor local sampling mechanism based on $(\\acute{\\pmb{X}}, \\pmb{Y}, \\pmb{Z})$. Define $\\chi\\_{\\alpha}^B:=\\\\{(\\pmb{x},\\pmb{x}^{(1)},\\ldots , \\pmb{x}^{(B)})\\Big| \\left[1+\\sum\\_{b=1}^B1(T(\\pmb{x}^{(b)},\\pmb{Y},\\pmb{Z})\\geq T(\\pmb{x},\\pmb{Y},\\pmb{Z}))\\right]\\big/(1+B)\\leq \\alpha \\\\}$. The type I error of our test  conditionally  on $\\pmb{Y}$ and $\\pmb{Z}$ can be expressed as $P((\\pmb{X},\\widetilde{\\pmb{X}}^{(1)},\\ldots, \\widetilde{\\pmb{X}}^{(B)})\\in \\chi_{\\alpha}^B|\\pmb{Y},\\pmb{Z})$. This  expression  can be decomposed into the sum of $P((\\acute{\\pmb{X}},\\acute{\\pmb{X}}^{(1)},\\ldots ,\\acute{\\pmb{X}}^{(B)})\\in \\chi_{\\alpha}^B|\\pmb{Y},\\pmb{Z})$ and $P((\\pmb{X},\\widetilde{\\pmb{X}}^{(1)},\\ldots ,\\widetilde{\\pmb{X}}^{(B)})\\in \\chi_{\\alpha}^B|\\pmb{Y},\\pmb{Z})-P((\\acute{\\pmb{X}},\\acute{\\pmb{X}}^{(1)},\\ldots ,\\acute{\\pmb{X}}^{(B)})\\in \\chi_{\\alpha}^B|\\pmb{Y},\\pmb{Z})$. On one hand, the exchangeability of $B+1$ triples $(\\acute{\\pmb{X}}, \\pmb{Y}, \\pmb{Z}),(\\acute{\\pmb{X}}^{(1)}, \\pmb{Y}, \\pmb{Z}),\\cdots,(\\acute{\\pmb{X}}^{(B)}, \\pmb{Y}, \\pmb{Z})$ under $H_0$, conditioned on $\\pmb{Y}$ and $\\pmb{Z}$, leads to the first term being smaller than $\\alpha$. On the other hand, $(\\acute{\\pmb{X}}^{(1)},\\ldots ,\\acute{\\pmb{X}}^{(B)})$, conditioned on $\\acute{\\pmb{X}},\\pmb{Y}$ and $\\pmb{Z}$, is generated from the same mechanism as $(\\widetilde{\\pmb{X}}^{(1)},\\ldots ,\\widetilde{\\pmb{X}}^{(B)})$, conditioned on $\\pmb{X},\\pmb{Y}$ and $\\pmb{Z}$. In other words, for all $\\pmb{x}\\in \\mathbb{R}^n$, $((\\widetilde{\\pmb{X}}^{(1)},\\ldots ,\\widetilde{\\pmb{X}}^{(B)})|\\pmb{X}=\\pmb{x},\\pmb{Y},\\pmb{Z})\\overset{\\mathrm{d}}{=}((\\acute{\\pmb{X}}^{(1)},\\ldots ,\\acute{\\pmb{X}}^{(B)})|\\acute{\\pmb{X}}=\\pmb{x},\\pmb{Y},\\pmb{Z})$. By Lemma 6 in the Supplementary Materials, we can bound the second term by the total variation distance between $\\widehat{p}(\\cdot|\\pmb{Z})$ and $p(\\cdot|\\pmb{Z})$.\n\nThird, we discuss the non-trivial points in the power analysis of our method in Theorem 4. We utilize the Markov inequality to establish a connection between power analysis and the consistency of the $\\mbox{CMI}$ estimator.  Mukherjee et al. (2020) considered the scenario where $p(x,z)p(y|z)$ (denoted as $g$) is known for the consistency.  However, this distribution is  unknown in practice.  We investigate the consistency  when incorporating the  $1$-NN sampling (Algorithm 1) to approximate it. In this extension, we denote the joint density of $(X,Y',Z)$ generated through $1$-NN sampling as $\\phi$. The CMI estimator $\\widehat{\\mbox{CMI}}$ based on Algorithm 2 is referred to as $\\widehat{D}\\_{KL}^{(n)}(f||\\phi)$, and the true CMI of $(X,Y,Z)$ is denoted as $\\mbox{CMI}:=I(X;Y|Z)=D\\_{KL}(f||g)$. The key aspect of this extension involves decomposing $\\widehat{D}\\_{KL}^{(n)}(f||\\phi)-D\\_{KL}(f||g)$ into two components: $\\widehat{D}\\_{KL}^{(n)}(f||\\phi)-D\\_{KL}(f||\\phi)$ and $D\\_{KL}(f||\\phi)-D\\_{KL}(f||g)$. The first term is shown to be $o\\_p(1)$ according to Mukherjee et al. (2020). By leveraging properties of point-wise minimizers of binary cross-entropy loss, Jensen's inequality, and Lagrange's mean value theorem, we bound the second term by $d\\_{TV}(\\phi,g)$, which is also $o(1)$. The combination of these two results leads to the conclusion that $\\widehat{D}\\_{KL}^{(n)}(f||\\phi)-D\\_{KL}(f||g)=o_p(1)$, which confirms the consistency of our CMI estimator.\n\n**Question 1**:\nPlease refer to the global response for details regarding the question \"Why is the CMI estimated using classifiers rather than directly optimising (4) like f-GAN? Why not directly use the maximum of (4) to approximate the CMI?\". \n"}, {"review_id": "yraNCF5XPK", "reviewer": "Reviewer_671y", "summary": "In this paper, the authors present a method for conditional independence (CI) testing based on a conditional mutual information (CMI) estimator. \nThis CMI estimator is based on a classifier over 1-NN samples. From the CMI estimates, the authors present a hypothesis testing procedure for CI based on multiple iterations over estimates of CMI on K-NN samples.\n\nThe authors present a theoretical analysis of their approach and the empirical results show that their approach has low Type 1 error and high power. The empirical section shows also that their approach remains consistent as the number of dimensions increase for the conditioning variables.", "questions": "I might have missed it but it wasn't clear to me, which classifier you used in your experiments?\n\nWhat do you mean by $V_0^{test}$ in line 7 of alg. 2? is that any instance over the union of test sets in f and g?\n\nHow long does your approach take compared to the others? After all, a large number of B repetitions could significantly add up.\n", "limitations": "I don't expect any potential negative societal impact from this work", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is well written and technically sound. The contribution builds on previous approaches but is an original method for the important task of CI. \nThe method is very promising given the theoretical and empirical results providing balanced detection of TP and TN links in real-world datasets as shown by the top F-scores. \nThe use of a single classifier for the CMI estimation is advantageous as the CI test is usually executed multiple times for graph detection.", "weaknesses": "I find the paper really compelling and useful in many other ML tasks. However, I believe the empirical evaluation could have gone a step further.\nIn particular, I believe a distillation procedure would be helpful to show which part of your approach contributes to the significant results. More precisely, it would be great to see how good your CMI estimator is. A comparison against two-classifier approaches and even a synthetic dataset where you have access to the ground truth conditional distributions would be very useful.\nThis could show how tight the bound is given your estimator. Furthermore, it would be nice to see the performance of your approach given different classifiers, e.g., DNNs or Random forest methods as mentioned in your paper.\n\nIn the presentation, I would make the figures larger. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the important review!\n\n**Weakness 1**.  To show which parts of our approach contribute to the significant results, we have done some additional experiments as per your suggestions. Specifically, the classifiers we evaluate include DNNs, XGBoost, and Random forest. The simulation dataset is generated according to 'Scenario I: The post-nonlinear model' in the paper. We maintain the distribution as the standard Gaussian, with the sample size $n=500$ and dimensions of $Z$ set to $10$ and $50$.\n\nFirst, we assess the effects of different classifiers on CMI estimation. Under the null hypothesis $H_0$, the true CMI value is zero. We calculate the mean squared errors (MSEs) of the CMI estimators using different classifiers under $H_0$. These MSE values are averaged over 100 independent trials. It is observed from Table 3 in the pdf file (attached in the global response) that the MSEs of the CMI estimators based on the three classifiers exhibit remarkable similarity. Furthermore, these three classifier-based CMI estimation methods all yield estimators that closely approximate the true value of $0$ under the null hypothesis.\n\nHowever, calculating the true value of CMI under $H_1$ is highly complex. Therefore, we do not compute the MSE in this case. Instead, we further utilize the thresholding approach to evaluate the performances. To be specific, we conduct the experiment independently 100 times under the assumption of $H_0$. While $H_0$ being true implies $I(X;Y|Z)=0$, the CMI estimator from a finite-sized dataset seldom equals zero due to statistical variation. We simply set the threshold at $0.1$ for the CMI estimator: if $\\hat{I}>0.1$, we reject $H_0$, and we then count the instances of falsely rejecting $H_0$. Similarly, under $H_1$, we threshold the CMI estimator at $0.1$ and record the number of correct reject $H_0$. Table 4 in the pdf file (attached in the global response) records the results obtained. We observed that, XGBoost-CMI  has the highest empirical type I errors $0.12 (0.19)$ and  empirical testing power $0.85 (0.83)$, and DNNs-CMI  has the lowest  empirical type I error $0.05 (0.07)$ and  empirical testing power $0.38 (0.12)$. In practice,  choosing an appropriate threshold value presents a significant challenge due to the inherent difficulty in accurately estimating the statistical variance of the CMI estimator. Therefore, relying on a threshold-based scheme for testing conditional independence is considered unwise. To enhance the effectiveness of CI testing based on the classifier-CMI estimator,  we propose in this paper applying the k-nearest-neighbor local sampling strategy to approximate the conditional distribution $X|Z$ that encodes the null hypothesis. This scheme takes into account the statistical variation inherent in the CMI estimation. It is not only simple and effective, but also has theoretical guarantee (See Theorem 2). \n\nSecond, we employ the same dataset utilized for CMI estimation to evaluate the performance of our CMI estimation method using various classifiers equipped with the k-nearest-neighbor local sampling approach. We present the type I error rate, testing power, and timing performance for a single test for each method in Table 5 in the pdf file (attached in the global response). These results are averaged over 100 independent trials. Upon examination of the results in Tables 4 and 5 in the pdf file (attached in the global response), we observe that when using the threshold value as the basis for testing conditional independence, the XGBoost-CMI method exhibits inflated type I error, while the DNNs-CMI method demonstrates low testing power. However, when incorporating the k-nearest-neighbor local sampling scheme, the performance of both methods improves. Specifically, the  empirical type I error rate of XGBoost-CMI decreases from $0.12 (0.17)$ to $0.06 (0.06)$, and its  empirical testing power increases from $0.85 (0.83)$ to $0.97 (0.96)$. Moreover, the  empirical testing power of DNNs-CMI  increases from $0.38 (0.12)$ to $0.85 (0.34)$. These results indicate that both a good classifier-based CMI estimator and the k-nearest-neighbor local sampling strategy are key factors contributing to the significant outcomes achieved by our test.\n\n**Weakness 2**. Thanks for your suggestion. We will make the figures larger in the final version.\n\n**Question 1**. The XGBoost classifier was used in all of our experiments.\n\n**Question 2**. After completing the classifier training in Algorithm 2 (line 7), we proceed to compute the prediction probability $P(l=1|w)$ for each feature $w$ in the test set $V^{test}$. The label of the sample $w$ (i.e., $l=1$ or $l=0$) is not required when computing $P(l=1|w)$. However, since $V^{test}$ contains labels for all samples, we introduce a new symbol $V_{0}^{test}$, which includes all the features in $V^{test}$ but excludes the labels.  $V_{0}^{test}$ is  essentially the union of the test sets in $f$ and $g$.\n\n**Question 3**. We report the timing performance of all methods for a single test in Figure 7 of the Supplementary Materials. Our test is found to be highly computationally efficient even when dealing with large sample sizes and high-dimensional conditioning sets. In contrast, CMIknn and CCIT for sample sizes exceeding 1000, and LPCIT for dimension of $Z$ higher than 50 are impractical due to their prohibitively long running time. For example, when $n=2000$ and $dz=80$, a single test of our method takes 42 seconds, GCIT 12 seconds, CCIT 173 seconds, KCIT 13 seconds, NNSCIT 31 seconds, CMIknn 305 seconds and LPCIT 188 seconds. In our approach, we set $B=200$. "}, {"review_id": "EV07UQ07Um", "reviewer": "Reviewer_P7Sn", "summary": "This paper argues that conditional independence (CI) testing becomes challenging because of high-dimensional conditioning variables and limited data samples. To address these issues, the authors propose a testing approach incorporating a class-based conditional mutual information (CMI) and a $k$-nearest-neighbor local sampling strategy. Theoretical analysis demonstrates its asymptotic control of type I error and consistency against all alternative hypotheses. Extensive empirical results on synthetic and real data show that the proposed method achieves computational efficiency, decent performance under different scenarios, and robustness towards heavy-tailed data.", "questions": "1. Can the authors add more discussions on the difference between the $k$-nearest-neighbor local sampling strategy and the binning strategy? What are the advantages of the $k$-nearest-neighbor local sampling strategy?\n2. Can the authors provide the optimal $k$ on real data?", "limitations": "Yes, the authors have addressed the societal impact of their work.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "This idea is intrigued by two significant challenges of CI testing, i.e., high-dimensional conditioning variables and limited data samples. \nThe authors do a good job of motivating this work, and the proposed test approach is simple but effective.\nExtensive empirical results show that the proposed method achieves a better trade-off between type I error rate and testing power.\nIn addition, the paper is very well written and easy to understand.", "weaknesses": "My main doubts/concerns regarding the paper are the following:\n- Line 104-105 claim that the $k$-nearest-neighbor local sampling strategy is an alternative to the binning strategy. However, this paper does not provide a comparison between the proposed testing approach and the methods employing a binning strategy [1,2].\n- In Figure 3, the results on synthetic data show that the performance is sensitive to the selection of $k$. Therefore, it is unreasonable to apply $k=7$ in all experiments, especially in real data.\n\nMinor:\n- Empirical analysis shows that the proposed method achieves a better trade-off between type I error rate and testing power. It is hard to say that the proposed method outperforms existing SOTA methods.\n\n[1] I. Kim, et al., \"Local Permutation Tests for Conditional Independence\", The Annals of Statistic 2022.\n\n[2] D. Margaritis, \"Distribution-free Learning of Bayesian Network Structure in Continuous Domains\", AAAI 2005.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the helpful review! \n\n**Weakness 1**. The binning strategy can only handle conditioning variables $Z$  with very few dimensions. Please see the response for Question 1. We therefore did not consider the comparison with the methods employing a binning strategy in the paper.\n\n**Weakness 2**.  Thank you for this important question. Conditional independence testing is an unsupervised problem, since the ground truth information regarding conditional independence or non-conditional independence is unavailable. Therefore, the cross-validation method does not work for determining the optimal $k$ in practical scenarios. The task of selecting the optimal $k$ appears to be quite challenging in the context of CIT. We will delve into developing methods to identify the optimal parameter $k$ in our future research. Nevertheless, utilizing $k=7$ in the experiments and the two real data analyses yields very good performances,  even though it may not represent the optimal choice for $k$. It is noted that,  in the works related to CIT, such as Runge (2018), Bellot and van der Schaar (2019), and Scetbon et al. (2022), the hyperparameters are chosen using synthetic data and then applied to other experiments and real data analyses.\n\n**Weakness (Minor)**: When assessing the effectiveness of a test, it's crucial that it showcases a higher statistical power compared to other methods, while simultaneously effectively controlling the type I error. On one hand, if a test fails to control the type I error well, it cannot be considered as reliable. On the other hand, when a test successfully controls the type I error, its power is then compared to that of other tests. If it exhibits superior power, it can be confidently regarded as more favorable and superior. Based on this perspective, it is observed from the experiments that our method outperforms the commonly used competitive SOTA methods. For this statement, we will replace “existing SOTA methods” with “the commonly used competitive SOTA methods”. \n\n**Question 1**. The binning strategy can only handle conditioning variables $Z$  with very few dimensions. Kim et al. (2022) proposed two binning strategies. The first binning strategy is to partition the conditional variable $Z$ into $M$ bins of equal size. Then, we permute the X-samples in each bin randomly and perform local permutation tests. The second binning strategy is to partition $Z$ into $M$ coarser bins and further partition each coarser bin into fine bins in which permutation occurred. Kim et al. (2022) used both strategies to process univariate $Z$ in all of the experiments. Margaritis (2005) adopted a recursive-median binning strategy and  can not handle the high-dimensional conditional variable $Z$ either. In fact, we have conducted simulation studies and confirmed that the binning strategy can only handle the condition variable $Z$ with very  few dimensions.\n\nWe give explanations for the limitations of the binning strategy in addressing the high-dimensional conditional variable $Z$.\nBinning serves as a technique to partition samples that exhibit notable similarity with respect to variable $Z$ into distinct bins. Shuffling $X$ within each bin maintains the relationship between $X$ and $Z$, as opposed to directly permuting $X$ throughout the entire dataset. Nevertheless, employing the binning method in high-dimensional scenarios poses challenges (Berrett et al., 2020). One challenge stems from the intricate task of selecting bins where the samples within each bin exhibit notable similarity concerning the variable $Z$. Another challenge arises from the generation of an excessive number of bins in high-dimensional $Z$. For instance, consider a scenario where the dimension of $Z$ is $100$. By applying the first binning strategy described in Kim et al. (2022) to partition each dimension of $Z$, we would generate an excessive number of bins for $Z$. Specifically, the number of bins is $M^{100}$. Even with a small value of $M=2$, the binary binning technique would produce a staggering $2^{100}$ bins. As highlighted by Kim et al. (2022), redundant bins can result in a reduction of test power. More specifically, a substantial number of bins within the high-dimensional space of Z consist of just one or two samples. On one hand, as per the definition of $U$-statistics in the computation of test statistics in Kim et al. (2022), only bins comprising a minimum of $4$ samples are included in the computation of $U$-statistics. Consequently, those bins containing just one or two samples do not contribute to the calculation of the test statistics. This leads to an underutilization of samples in these bins. On the other hand, significant similarity among various permutations can arise when numerous bins consist of merely one or two samples. This similarity can result in insufficient statistical power, and may also undermine effective control of the type I error rate. \n\nDifferent from the binning strategy involving permutation, the k-nearest-neighbor local sampling strategy serves as an alternative sampling approach. This approach equips us with pseudo samples $\\widetilde{\\pmb{X}}^{(b)}$ ($b=1,\\cdots,B$) that exhibit considerable dissimilarity while closely approximating the true conditional distribution in terms of the total variation distance, as evidenced by our Lemma 1. Especially in scenarios with high-dimensional $Z$, the k-nearest-neighbor local sampling strategy, compared to the permutation-based binning approach, effectively resolves concerns like maintaining the relationship between $X$ and $Z$ for pseudo sample $\\tilde{X}$, as well as addressing issues of sample underutilization and pronounced similarity among permutations. This results in adequate statistical power and good control of the type I error rate within the context of high-dimensional $Z$.\n\n**Question 2**. Please see our response to Weakness 2. "}]

Author rebuttal:
We thank the reviewers for their detailed and considerate feedback. We offer the following clarifications.  **A PDF file including 5 tables is attached.** 

1. About the claim "**our method avoids sample splitting**", we  emphasize that 
our proposed method  **allows the entire dataset to be used in computing
$\widehat{\mbox{CMI}}$**. The method in Li et al. (2023) requires splitting the dataset into two parts, and thus the dataset used for calculating the test statistics comprises only one-third of the total samples. This will reduce the statistical power of the  test, particularly when working with small datasets. In contrast, our proposed procedure  **allows the entire dataset to be used in computing $\widehat{\mbox{CMI}}$**, thereby avoiding the loss of testing power. In the computation of the CMI estimator using the entire dataset, we adopt a classifier-based approach. This approach necessitates the utilization of two datasets for a  supervised classification task. As a result, the data splitting in this context occurs naturally and does not result in any loss of testing power.  

2. Regarding the claim  "**high-dimensional conditioning variables**", we have the following response. Recently, Cai et al.(2022), Kim et al. (2022), and Ai et al. (2022) have primarily concentrated on low-dimensional cases where the dimension of $Z$ does not exceed 2. Additionally, Runge (2018) employed a sample size of 500, while restricting the dimension of $Z$ to a maximum of 8. CI testing is a really challenging task when handling high-dimensional conditioning variable sets (Scetbon et al., 2022). The concept "high dimensional" in CI testing refers to the scenario where the dimension of $Z$ is not low or is relatively high given the sample size (Scetbon et al., 2022). This concept is different from that of the high-dimensional regression, where some sparsity assumptions may be needed. 

For example, Sen et al. (2017) and Scetbon et al. (2022) both claimed they can handle the high-dimensionality of $Z$. In their simulation studies, Sen et al. (2017) set the sample size to be 1000 with the dimension of $Z$ being at most 150, and Scetbon et al. (2022) set the sample size to be 1000 with the dimension of $Z$ being at most 50. They also did not consider that the dimension of $Z$ converges to infinite with the sample size in their theoretical results. In Figure 1 of our paper, we present the favorable performance of our tests for $n=300$ and the dimension of $Z$ equal to 80. In response to Reviewer BxQV, we also consider $Z$ dimension as high as $500$ when $n=500$. In contrast to previous works, the dimension of $Z$ is relatively high.

[1] Cai et al., A Distribution Free Conditional Independence Test with Applications to Causal Discovery, JMLR 2022.

[2] Ai et al., Testing Unconditional and Conditional Independence Via Mutual Information, Journal of econometrics 2022.

3.  Regarding the claim "**our method operates without assumptions about distribution forms or feature dependencies**", we emphasize that, unlike Candes et al. (2018), our method operates successfully in practical applications without requiring the true distribution of $X|Z$ to be provided. We propose using the k-nearest-neighbor local sampling strategy to approximate  $p\_{X|Z}(x|z)$.  We establish in Section 4 that the total variation distance between the true distribution of $X|Z$ and the distribution of samples generated by the k-nearest-neighbor local sampling strategy tends to zero in probability as $n$ goes to infinity. 

However, in the theoretical part, we need Assumptions 1 and 2 to facilitate our proof. The two assumptions  have been included in Gao et al. (2016, 2017) and Sen et al. (2017), even though they may not be the weakest conditions. Moreover, we can validate Assumption 1 when ($X, Z$) follows the multivariate Gaussian distribution  (MVD) and Assumption 2 when $Z$ follows the MVD. 

4. About the question "**Why is the CMI estimated using classifiers rather than directly optimising (4) like f-GAN? Why not directly use the maximum of (4) to approximate the CMI?**", we have the following response. Belghazi et al. (2018) introduced a mutual information neural estimator (MINE) by maximizing equation (4) directly. We obtained the conditional version of this estimator to estimate CMI, referred to as C-MINE. This estimator does not require dataset splitting into training and testing sets. However, through simulations, we observed that C-MINE has a high variance and its optimization becomes unstable in high-dimensional settings. Similar conclusions have been drawn by Poole et al. (2019) and Mukherjee et al. (2020). Furthermore, C-MINE exhibits more computational complexity compared to our classifier-based CMI estimator. As outlined in Algorithm 3, to obtain a single p-value, we need to estimate $B+1$ CMI values. Therefore, any method that can offer fast and accurate estimation will yield substantial benefits.

Furthermore, we note that  Mondal et al. (2020) directly optimized (4) like f-GAN to estimate CMI. However, the training of GANs is often  challenging, with the risk of collapse if hyperparameters and regularizers are not carefully chosen (Dhariwal and Nichol, 2021). In contrast, our classifier-based CMI estimation method relies on a simple binary classifier and performs effectively even in high-dimensional scenarios. Importantly, our classifier-based approach demonstrates remarkable computational efficiency compared to the direct optimization of equation(4) by f-GAN.

[1] Belghazi et al., Mutual Information Neural Estimation, ICML 2018.

[2] Poole et al., On Variational Bounds of Mutual Information, ICML 2019.

[3] Kumar Mondal et al., C-MI-GAN: Estimation of Conditional Mutual Information Using MinMax Formulation, UAI 2020.

[4] Prafulla Dhariwal and Alex Nichol, Diffusion Models Beat GANs on Image Synthesis, NeurIPS 2021.

Meta-review: The paper presents a novel approach to Conditional Independence testing, leveraging kNN as a sampling method, which enhances computational efficiency. The work is theoretically sound, backed up well by theorems, and the algorithm is straightforward to implement. The paper is well-written and sufficiently readable.  While there are some concerns about the extent of its advance to the existing works, these do not overshadow the paper's strengths. Therefore, the recommendation is to accept the paper.
Decision: Accept (poster)