###### Abstract

This paper is concerned with the problem of conditional independence testing
for discrete data. In recent years, researchers have shed new light on this
fundamental problem, emphasizing finite-sample optimality. The non-asymptotic
viewpoint adapted in these works has led to novel conditional independence
tests that enjoy certain optimality under various regimes. Despite their
attractive theoretical properties, the considered tests are not necessarily
practical, relying on a Poissonization trick and unspecified constants in
their critical values. In this work, we attempt to bridge the gap between
theory and practice by reproving optimality without Poissonization and
calibrating tests using Monte Carlo permutations. Along the way, we also prove
that classical asymptotic $\chi^{2}$\- and $G$-tests are notably sub-optimal
in a high-dimensional regime, which justifies the demand for new tools. Our
theoretical results are complemented by experiments on both simulated and
real-world datasets. Accompanying this paper is an R package UCI that
implements the proposed tests.

Conditional Independence Testing for Discrete Distributions: Beyond
$\chi^{2}$\- and $G$-tests

Ilmun Kim† Matey Neykov‡ Sivaraman Balakrishnan‡ Larry Wasserman‡

††footnotetext: The last three authors‡ are listed randomly. †Department of
Statistics and Data Science, Yonsei University
---
‡Department of Statistics and Data Science, Carnegie Mellon University

August 31, 2024

##  1 Introduction

Conditional independence (CI) is the backbone of diverse fields in statistics,
including graphical models (De Campos and Huete,, 2000; Koller and Friedman,,
2009) and causal inference (Spohn,, 1994; Pearl,, 2014; Imbens and Rubin,,
2015). Among several benefits, this fundamental assumption allows us to
simplify the structure of a model, thereby increasing interpretability and
reducing computational costs. To justify the use of CI assumption, it is of
considerable interest to test whether two random variables $X$ and $Y$ are
independent after accounting for the effect of another random variable $Z$.
Due to its important role, the problem of CI testing has received much
attention in the past decade, resulting in numerous exciting new developments
(e.g., Candés et al.,, 2018; Shah and Peters,, 2020; Berrett et al.,, 2020;
Neykov et al.,, 2021; Petersen and Hansen,, 2021; Liu et al.,, 2022; Lundborg
et al.,, 2022). See Li and Fan, (2020) and Chatterjee, (2022) for recent
reviews. However, most of the recent work is dedicated to continuous data and
the importance of discrete CI testing is relatively overlooked.

In discrete settings, two commonly used methods are the $\chi^{2}$-test
(Pearson,, 1900) and the $G$-test (McDonald,, 2014), and their asymptotic
equivalence is well-known under regularity conditions (e.g., Chapter 14 of
Bishop et al.,, 2007). When $X$ and $Y$ are binary, the
Cochran–Mantel–Haenszel test (Agresti,, 2003) is another popular method for CI
testing. Despite their popularity, these methods are asymptotic in nature,
frequently calibrated by their limiting null distributions. Therefore their
finite-sample validity remains questionable. This miscalibration issue becomes
more serious in high-dimensional regimes where the number of categories can be
significantly larger than the sample size. Besides, the power of these methods
is not well-understood except in classical fixed-dimensional settings.

For discrete CI testing, Canonne et al., (2018) put forward two testing
algorithms and analyze their sample complexity from a non-asymptotic
perspective. Their sample complexity results are further complemented by
matching lower bounds, demonstrating optimality of their procedures in some
regimes. In spite of these technical advances, their approach poses several
practical challenges. First, their results rely on a Poissonization trick
where the sample size is treated as a Poisson random variable. This assumption
greatly simplifies the theoretical analysis, but is untenable in practice.
Another issue worth highlighting is the dependence of the test on unspecified
constants in their critical values. In many statistical applications, the type
I error is a greater concern than the type II error. It is therefore desirable
to set a critical value in such a way as to maximize the power, while
_tightly_ controlling the type I error. However, it is unclear from Canonne et
al., (2018) how to modify their tests to meet this criteria, thereby leaving
room for improvement from a practical perspective. Indeed, this issue was the
main motivation of recent work of Kim et al., 2022a ; Kim et al., 2022b  that
advocates the use of permutation methods in two-sample and (both unconditional
and conditional) independence testing problems.

With these issues in mind, our work makes the following contributions: (i) In
Theorem 1, we depoissonize the sample complexity results of Canonne et al.,
(2018) and establish the same theoretical guarantees under the standard
sampling setting. On a technical level, the challenge lies in dealing with the
complicated dependence structure of multinomial samples. We overcome this
difficulty using the negative association property of multinomial
distributions (Joag-Dev and Proschan,, 1983). (ii) We further make the
algorithms of Canonne et al., (2018) practical by leveraging the permutation
method to calibrate test statistics. This resampling approach completely
removes the issue arising from unspecified constants, and provably controls
the finite-sample type I error. In Theorem 2, we prove that Monte Carlo
permutation tests achieve the same sample complexity as the theoretical tests
of Canonne et al., (2018). (iii) The considered test statistics are linear
combinations of fourth order U-statistics, which can be daunting
computationally. We address this computational concern by presenting
alternative linear time expressions in Proposition 1. (iv) We also prove an
independent result that demonstrates sub-optimality of asymptotic $\chi^{2}$\-
and $G$-tests in their power performance. This negative result naturally
inspires efforts to develop new CI tests that perform better than the
classical ones. (v) Finally, we provide extensive simulation results that
demonstrate the practical value of the proposed methods in Section 5, and the
algorithms are available in the R package UCI.111publicly available at github
repository: https://github.com/ilmunk/UCI

Our work is related to Tsamardinos and Borboudakis, (2010) who warn about the
risk of asymptotic calibration for $\chi^{2}$\- and $G$-tests, and further
highlight benefits of the permutation procedure in type I error control. The
risk of asymptotic calibration has also been discussed in other testing
problems, such as those studied in Balakrishnan and Wasserman, (2018, 2019);
Kim et al., 2022a . In line with this research, we prove the negative result
of asymptotic $\chi^{2}$\- and $G$-tests, and demonstrate attractive
properties of the permutation method both in type I and II error control.
Another related work is Berrett and Samworth, (2021) where the authors propose
a permutation test based on a U-statistic for unconditional independence
testing. Concurring with our view, Berrett and Samworth, (2021) put an
emphasis on the permutation approach for practical calibration and demonstrate
the competitive performance of their proposal, coined USP test, over
$\chi^{2}$\- and $G$-tests. In fact, when the conditional variable is
degenerate (i.e., $Z$ takes a single value), one of our practical proposals
becomes exactly the same as that of Berrett and Samworth, (2021). In this
sense, our work can be considered as an extension of Berrett and Samworth,
(2021) to CI testing. We also refer to Agresti, (1992); Yao and Tritchler,
(1993) that discuss exact inference methods for contingency tables. It is
worth pointing out that the current paper builds on our prior work (Kim et
al., 2022b, ), which proves that the sample complexity results of Canonne et
al., (2018) continue to hold using permutation tests. However, the analysis of
Kim et al., 2022b  relies on Poissonization and also makes use of
(computationally expensive) full permutation tests. The current work deviates
from Kim et al., 2022b  by removing Poissonization and employing a more
computationally efficient permutation test via Monte Carlo sampling. We also
propose a new permutation test, called wUCI-test, that avoids sample splitting
and illustrate its competitive finite sample performance under a variety of
settings.

##### Organization.

The rest of this paper is organized as follows. In Section 2, we set the stage
by presenting some background information on sample complexity and
Poissonization. Section 3 describes the test statistics that we study, and
verifies that they can be computed in linear time. Section 4 contains our main
theoretical results including depoissonization and sub-optimality of
$\chi^{2}$\- and $G$-tests. In Section 5, we demonstrate the empirical
performance of the proposed methods based on simulated and real-world
datasets, before concluding in Section 6. All the proofs of our results are
relegated to the Appendix.

##### Notation.

For a positive integer $a$, we use the shorthand $[a]=\\{1,\ldots,a\\}$. The
conditional independence of $X$ and $Y$ given $Z$ is denoted as
$X\,\rotatebox[origin={c}]{90.0}{$\models$}\,Y\,|\,Z$. Given two discrete
distributions $p$ and $q$, we write the $L_{1}$ distance between $p$ and $q$
as $\|p-q\|_{1}$. We say that random variables $X_{1},\ldots,X_{n}$ are i.i.d.
when they are independent and identically distributed. For two positive
sequences $a_{n}$ and $b_{n}$, we write $a_{n}\asymp b_{n}$ if it holds that
$C_{1}\leq a_{n}/b_{n}\leq C_{2}$ for some positive constants $C_{1}$ and
$C_{2}$, and for all $n$. We also write $a_{n}=O(b_{n})$ or $a_{n}\lesssim
b_{n}$ to indicate that $a_{n}\leq Cb_{n}$ for some positive constant $C$
independent of $n$.

##  2 Background

Before presenting our main results, we start by building some background
knowledge on sample complexity and Poissonization.

###  2.1 Setting the stage

Consider the set of discrete distributions of $(X,Y,Z)$ on a domain
$[\ell_{1}]\times[\ell_{2}]\times[d]$, denoted by $\mathcal{P}$. Let
$\mathcal{P}_{0}$ be the subset of $\mathcal{P}$ such that
$X\,\rotatebox[origin={c}]{90.0}{$\models$}\,Y\,|\,Z$. Given $n$ i.i.d. random
vectors $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}$ drawn from
$p_{X,Y,Z}\in\mathcal{P}$, our goal is to distinguish

$\displaystyle H_{0}:p_{X,Y,Z}\in\mathcal{P}_{0}\quad\text{versus}\quad
H_{1}:p_{X,Y,Z}\in\mathcal{P}_{1}(\varepsilon)=\Big{\\{}p\in\mathcal{P}:\inf_{q\in\mathcal{P}_{0}}\|p-q\|_{1}\geq\varepsilon\Big{\\}},$
(1)

where $\varepsilon>0$ is a distance parameter (Figure 1 for a pictorial
description). When the sample size $n$ is too small, no valid test can
reliably differentiate the null from the alternative. On the other hand, when
the sample size $n$ is too large, the problem becomes trivial, resulting in
many successful tests. A natural question is then to determine the smallest
$n$ required to achieve the desired level of testing accuracy for an optimal
test. This concept is known as the optimal sample complexity formally defined
below.

Figure 1: Schematic of the hypotheses of CI testing. The space of null
distributions is $\varepsilon$ far from the space of alternative distributions
where the distance parameter $\varepsilon$ controls the difficulty of the
problem.

##### Optimal sample complexity.

We define a test $\phi$ which maps from the samples to a binary outcome as
$\phi:\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}\mapsto\\{0,1\\}$. For some fixed
$\alpha\in(0,1)$, let $\Phi_{\alpha}$ denote the set of level $\alpha$ tests
such that for each $\phi\in\Phi_{\alpha}$,

$\displaystyle\sup_{p\in\mathcal{P}_{0}}\mathbb{P}_{p}(\phi=1)\leq\alpha.$

The minimax risk is the worst-case type II error of an optimal level $\alpha$
test defined as

$\displaystyle
R_{n}(\varepsilon,\alpha)=\inf_{\phi\in\Phi_{\alpha}}\sup_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(\phi=0).$

The difficulty of the problem can be characterized as the minimum number of
samples that makes the minimax risk bounded by some fixed constant
$\beta\in(0,1-\alpha)$222Throughout this paper, we treat the target type I and
II errors $\alpha$ and $\beta$ as universal constants, e.g.,
$\alpha=\beta=0.05$.. This minimum number of samples is called the optimal
sample complexity given as

$\displaystyle
n^{\star}=\inf\Big{\\{}n:R_{n}(\varepsilon,\alpha)\leq\beta\Big{\\}}.$

We say that a level $\alpha$ test $\phi$ is rate optimal in sample complexity
if

$\displaystyle
n^{\star}\asymp\inf\bigg{\\{}n:\sup_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(\phi=0)\leq\beta\bigg{\\}}.$

In practice, an optimal test whose risk is exactly equal to
$R_{n}(\varepsilon,\alpha)$ is mostly inaccessible. For this reason, we
instead aim to find a rate optimal test.

######  Remark 1 (Choice of metrics).

As is well-known, one half of the $L_{1}$ distance between two distributions
$p$ and $q$ is the same as their total variation (TV) distance:

$\displaystyle\frac{1}{2}\|p-q\|_{1}=\mathrm{TV}(p,q)=\sup_{A}|p(A)-q(A)|,$

where the supremum is taken over all possible measurable sets. The optimal
sample complexity depend crucially on the choice of metrics in the definition
of the alternative hypothesis (1). Since our aim is simply to depoissonize the
previous results, we focus on the $L_{1}$ distance or equivalently the TV
distance throughout this paper.

###  2.2 Poissonization

Poissonization is now a standard technique to study the sample complexity in
hypothesis testing (e.g., Valiant,, 2011; Chan et al.,, 2014; Valiant and
Valiant,, 2017; Balakrishnan and Wasserman,, 2019). The idea itself is old in
statistics and probability theory, dating back at least to Kac, (1949). See
also Chapter 3.5 of Vaart and Wellner, (1996). It is used in theoretical works
on hypothesis testing as a trick to simplify several calculations in dealing
with categorical data. In particular, it is well-known that when the data are
generated by Poisson sampling (Algorithm 1), the number of samples falling
into disjoint sets are mutually independent. This independence property lies
at the heart of deriving various existing results of the sample complexity in
distribution testing. See Canonne et al., (2022) for a recent review.

Algorithm 1 Poisson Sampling

Input: For a fixed $n\in\mathbb{N}$ and a distribution $p_{X,Y,Z}$

  1. 1.

Draw $N\sim\mathrm{Poisson}(n)$.

  2. 2.

Generate i.i.d. $(X_{1},Y_{1},Z_{1}),\ldots,(X_{N},Y_{N},Z_{N})$ random
variables from $p_{X,Y,Z}$.

Return: $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{N}$

###  2.3 General recipe and related issue

If constant factors are not the main concern, there are straightforward ways
of transferring the sample complexity obtained from Poisson sampling to the
usual sampling scenario with a fixed sample size. One concrete procedure,
described in Neykov et al., (2021), is as follows. Given $n$ i.i.d. copies of
$(X,Y,Z)$,

  1. 1.

Draw $N\sim\mathrm{Poisson}\bigl{(}\frac{n}{2}\bigr{)}$.

  2. 2.

If $N>n$, then accept the null hypothesis.

  3. 3.

If $N\leq n$, perform a test based on
$(X_{1},Y_{1},Z_{1}),\ldots,(X_{N},Y_{N},Z_{N})$ and return a result.

Let $\phi$ be a generic test function using
$(X_{1},Y_{1},Z_{1}),\ldots,(X_{N},Y_{N},Z_{N})$ and then the resulting test
from the above procedure can be concretely written as
$\phi^{\ast}=\mathds{1}(N\leq n)\phi$. Suppose that the test $\phi$ has the
type I error as well as the type II error bounded by $\alpha$ and $\beta$,
respectively. Then the corresponding test $\phi^{\ast}$ has the type I and II
error bounds as

$\displaystyle\sup_{p\in\mathcal{P}_{0}}\mathbb{E}_{p}[\phi^{\ast}]\leq\alpha\quad\text{and}\quad\sup_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{E}_{p}[1-\phi^{\ast}]\leq\beta+\mathbb{P}(N>n).$

Since the Poisson distribution is tightly concentrated around its mean, the
additional term in the type II error of $\phi^{\ast}$ can be made small when
$n$ is relatively large compared to $\beta$.333More precisely, an exponential
tail bound for a Poisson random variable ensures that $\mathbb{P}(N>n)\leq
e^{-n/8}$ (e.g., Theorem A.8 of Canonne et al.,, 2022) so that if $n\geq
8\log(1/\beta)$, the type II error is bounded above by $2\beta$. Therefore one
can transfer the sample complexity obtained under Poissonization to the usual
sampling scheme up to a constant factor.

Nevertheless, due to an inefficient use of the data as well as a non-trivial
chance of accepting the null irrelevant to the data, practitioners may not
necessarily follow this general recipe. To address this concern, there has
been an effort to depoissonize the results of distribution testing under
Poisson sampling. For instance, Kim, (2020) revisits the truncated $\chi^{2}$
test for goodness-of-fit testing proposed by Balakrishnan and Wasserman,
(2019) and establishes the same optimality without Poissonization. We also
refer to Jacquet and Szpankowski, (1998) and Chapter 2.5 of Penrose, (2003)
that present useful depoissonization tools, but mostly for asymptotic results.
In this work, we analyze the tests proposed by Canonne et al., (2018) under
the usual i.i.d. sampling model with a fixed sample size and show that the
same sample complexity can be achieved without any further assumption. At the
heart of our technique is the negative association of multinomial
distributions (Joag-Dev and Proschan,, 1983), which would be potentially
useful for depoissonizing other sample complexity results.

Multinomial sampling. To fix terminology, we refer to the usual sampling with
fixed sample size $n$ as _multinomial sampling_ in what follows.

##  3 Test statistics

To describe our main results, we first need to recall the test statistics
introduced by Canonne et al., (2018). As noted by Neykov et al., (2021), their
test statistics can be viewed as linear combinations of U-statistics
constructed using the observations of $(X,Y)$ in the same category of $Z$. We
describe two kinds of U-statistics considered in Canonne et al., (2018).

###  3.1 Unweighted U-statistic

Suppose that we observe $\sigma\geq 4$ i.i.d. observations of $(X,Y)$
supported on $[\ell_{1}]\times[\ell_{2}]$. We let $p_{X,Y}$ denote the joint
discrete distribution of $(X,Y)$ and let $p_{X}$ and $p_{Y}$ denote the
marginal distribution of $X$ and $Y$, respectively. The first U-statistic is
an unbiased estimator of the squared $L_{2}$ distance between $p_{XY}$ and
$p_{X}p_{Y}$. In more detail, borrowing the notation from Neykov et al.,
(2021), let

$\displaystyle\phi_{ij}(qr)=\mathds{1}(X_{i}=q)\mathds{1}(Y_{i}=r)-\mathds{1}(X_{i}=q)\mathds{1}(Y_{j}=r).$

For four distinct observations indexed by $i,j,k,l\in[\sigma]$, the kernel of
the unweighted U-statistic is defined as

$\displaystyle
h_{ijkl}=\frac{1}{4!}\sum_{(\pi_{1},\pi_{2},\pi_{3},\pi_{4})\in\Pi}\sum_{q\in[\ell_{1}],r\in[\ell_{2}]}\phi_{\pi_{1}\pi_{2}}(qr)\phi_{\pi_{3}\pi_{4}}(qr),$

where $\Pi$ is the set of all possible permutations of $(i,j,k,l)$. By the
linearity of expectations, we see that $h_{ijkl}$ is an unbiased estimator of
the squared $L_{2}$ distance between $p_{XY}$ and $p_{X}p_{Y}$. Given this
kernel and the dataset $D=\\{(X_{1},Y_{1}),\ldots,(X_{\sigma},Y_{\sigma})\\}$,
the unweighted U-statistic is computed as

$\displaystyle
U(D)=\binom{\sigma}{4}^{-1}\sum_{i<j<k<l:(i,j,k,l)\in[\sigma]}h_{ijkl}.$ (2)

It is worth pointing out that $U(D)$ is equivalent to the U-statistic for
unconditional independence testing considered in Berrett et al., (2021);
Berrett and Samworth, (2021); Kim et al., 2022a . Denote the datasets
$D_{m}=\\{(X_{i},Y_{i}):Z_{i}=m,i\in[n]\\}$ and write the sample size within
$D_{m}$ by $\sigma_{m}$. The final test statistic for CI testing based on
$U(D)$ is then constructed as

$\displaystyle T=\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}U(D_{m}).$

###  3.2 Weighted U-statistic

The previous unweighted U-statistic may suffer from high variance especially
when the $L_{2}$ norms of $p_{XY}$ and $p_{X}p_{Y}$ are large, resulting in
sub-optimal performance (see the simulation result in Figure 2 under Scenario
1). To address this issue, Canonne et al., (2018) employ the flattening idea
proposed by Diakonikolas and Kane, (2016). This approach involves partitioning
heavier bins into multiple smaller pieces in a way to reduce the $L_{2}$ norms
of the transformed distributions, leading to a smaller variance of the test
statistic. As noted in Balakrishnan and Wasserman, (2018) and further
elaborated by Neykov et al., (2021), the flattening procedure is equivalent to
using a carefully designed weighted kernel for a U-statistic.

To describe the weighted U-statistic, assume that the sample size
$\sigma=4+4t$ for some $t\in\mathbb{N}$, and let $t_{1}=\min(t,\ell_{1})$ and
$t_{2}=\min(t,\ell_{2})$. The construction of the weighted U-statistic
involves sample splitting where the dataset
$D=\\{(X_{1},Y_{1}),\ldots,(X_{\sigma},Y_{\sigma})\\}$ is split into
$D_{X}=\\{X_{i}:i\in[t_{1}]\\}$, $D_{Y}=\\{Y_{i}:t_{1}+1\leq i\leq
t_{1}+t_{2}\\}$ and $D_{X,Y}=\\{(X_{i},Y_{i}):2t+1\leq i\leq\sigma\\}$ of
sizes $t_{1}$, $t_{2}$ and $4+2t$, respectively. Define $a_{q}$ (and
$a_{r}^{\prime}$) as the number of observations equal to $q$ (and $r$) in
$D_{X}$ (and $D_{Y}$). For four distinct observations indexed by $i,j,k,l$ in
$D_{X,Y}$, consider a weighted kernel given as

$\displaystyle
h_{ijkl}^{\bm{a}}=\frac{1}{4!}\sum_{(\pi_{1},\pi_{2},\pi_{3},\pi_{4})\in\Pi}\sum_{q\in[\ell_{1}],r\in[\ell_{2}]}\frac{\phi_{\pi_{1}\pi_{2}}(qr)\phi_{\pi_{3}\pi_{4}}(qr)}{(1+a_{q})(1+a_{r}^{\prime})}.$

Due to the independence among the split datasets, the conditional expectation
of $h_{ijkl}^{\bm{a}}$ taken over $D_{X,Y}$ is the square of the $L_{2}$
distance between $p_{XY}$ and $p_{X}p_{Y}$ weighted by
$(1+a_{q})(1+a_{r}^{\prime})$. Given the kernel, the weighted U-statistic is
computed as

$\displaystyle
U_{W}^{\bm{a}}(D)=\binom{2t+4}{4}^{-1}\sum_{i<j<k<l:(i,j,k,l)\in
D_{X,Y}}h_{ijkl}^{\bm{a}},$ (3)

where $(i,j,k,l)\in D_{X,Y}$ indicates taking four observations indexed by
$(i,j,k,l)$ from the dataset $D_{X,Y}$. As before, denote the datasets
$D_{m}=\\{(X_{i},Y_{i}):Z_{i}=m,i\in[n]\\}$ and write the sample size within
$D_{m}$ by $\sigma_{m}$. By further writing
$\omega_{m}=\sqrt{\min(\sigma_{m},\ell_{1})\min(\sigma_{m},\ell_{2})}$, the
final test statistic for CI testing based on $U_{W}^{\bm{a}}(D)$ is then
constructed as

$\displaystyle T_{W}=\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{a}}(D_{m}).$

##### Practical considerations.

As shown in Canonne et al., (2018), the test based on $T_{W}$ achieves the
optimal sample complexity in broader regimes than that based on $T$. Despite
its attractive theoretical properties, the resulting test may experience a
loss of practical power due to an inefficient use of the data arising from
sample splitting. To mitigate this issue, we introduce another weighted
U-statistic without sample splitting. First, it is worth trying to understand
the weight $(1+a_{q})(1+a_{r}^{\prime})$ at a population level. In particular,
the conditional expectation of the weight given the random sample size $t$ is
$\bigl{(}1+\min\\{t,\ell_{1}\\}p_{X}(q)\bigr{)}\bigl{(}1+\min\\{t,\ell_{1}\\}p_{Y}(r)\bigr{)}$.
This motivates us to consider a weight kernel

$\displaystyle
h_{ijkl}^{\bm{b}}=\frac{1}{4!}\sum_{(\pi_{1},\pi_{2},\pi_{3},\pi_{4})\in\Pi}\sum_{q\in[\ell_{1}],r\in[\ell_{2}]}\frac{\phi_{\pi_{1}\pi_{2}}(qr)\phi_{\pi_{3}\pi_{4}}(qr)}{(1+b_{q})(1+b_{r}^{\prime})},$

where
$b_{q}=\min\\{\sigma,\ell_{1}\\}\sigma^{-1}\sum_{i=1}^{\sigma}\mathds{1}(X_{i}=q)$
and
$b_{r}^{\prime}=\min\\{\sigma,\ell_{2}\\}\sigma^{-1}\sum_{i=1}^{\sigma}\mathds{1}(Y_{i}=r)$.
The resulting weighted U-statistic is then computed as

$\displaystyle
U_{W}^{\bm{b}}=\binom{\sigma}{4}^{-1}\sum_{i<j<k<l:(i,j,k,l)\in[\sigma]}h_{ijkl}^{\bm{b}}.$
(4)

Similarly as before, the final CI test statistic based on $U_{W}^{\bm{b}}(D)$
is given as

$\displaystyle T_{W}^{\dagger}=\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{b}}(D_{m}).$

We coin the permutation test based on $T_{W}^{\dagger}$ as wUCI-test, and
formally introduce the procedure in Algorithm 4 of Appendix A.

######  Remark 2 (Connection with the truncated $\chi^{2}$-test).

As pointed out by several authors (Haberman,, 1988; Balakrishnan and
Wasserman,, 2019; Kim,, 2020), the classical $\chi^{2}$-test for goodness-of-
fit testing can easily break down for sparse multinomial data. To address this
problem, Balakrishnan and Wasserman, (2019) introduce a modification of the
$\chi^{2}$-test by using a truncated weight function and prove its minimax
optimality. Interestingly, the weight $(1+b_{q})(1+b_{r}^{\prime})$ that we
consider is closely connected to the truncated weight of Balakrishnan and
Wasserman, (2019) and may be regarded as an empirical counterpart for
independence testing. To explain, let us simply focus on the first term in the
product weight and notice that

$\displaystyle\ell_{1}^{-1}(1+b_{q})=\frac{1}{\ell_{1}}+\min\biggl{(}\frac{\sigma}{\ell_{1}},1\biggr{)}\widehat{p}_{X}(q)\asymp\max\biggl{\\{}\frac{1}{\ell_{1}},\min\biggl{(}\frac{\sigma}{\ell_{1}},1\biggr{)}\widehat{p}_{X}(q)\biggr{\\}},$

where $\widehat{p}_{X}(q)=\sigma^{-1}\sum_{i=1}^{\sigma}\mathds{1}(X_{i}=q)$.
When $\sigma\ell_{1}^{-1}$ is large, the right-hand side of the above equation
approximates $\max\\{\ell_{1}^{-1},p_{X}(q)\\}$, which is exactly the same as
the truncated weight in Balakrishnan and Wasserman, (2019) for goodness-of-fit
testing.

###  3.3 Linear time expression

The original forms of the aforementioned U-statistics take
$O(\sigma^{4}\ell_{1}\ell_{2})$ time to compute, which can be prohibitive for
large $\sigma,\ell_{1},\ell_{2}$. Luckily, this computational complexity can
be reduced to $O(\sigma)$ by exploiting a contingency table representation. A
computationally convenient form of the unweighted U-statistic $U(D)$ is
already given by Canonne et al., (2018); Berrett et al., (2021); Berrett and
Samworth, (2021). We also note that an alternative form of $U_{W}^{\bm{a}}(D)$
is provided in Canonne et al., (2018), but a naive calculation of their
expression takes at least $O(\sigma\ell_{1}^{2}\ell_{2}^{2})$ time. Here we
present a general expression for the U-statistics and explain that it can be
run in linear time _independent of_ $\ell_{1}$ and $\ell_{2}$.

To this end, let us set some notation. Let
$\bm{\eta}=\\{\eta_{1},\ldots,\eta_{\ell_{1}}\\}$ and
$\bm{\upsilon}=\\{\upsilon_{1},\ldots,\upsilon_{\ell_{2}}\\}$ be some weight
vectors with non-zero components. Given the dataset
$D=\\{(X_{1},Y_{1}),\ldots,(X_{\sigma},Y_{\sigma})\\}$, consider four distinct
observations indexed by $i,j,k,l\in[\sigma]$ and define the weighted kernel
associated with $\bm{\eta}$ and $\bm{\upsilon}$ as

$\displaystyle
h_{ijkl}^{\bm{\eta},\bm{\upsilon}}=\frac{1}{4!}\sum_{(\pi_{1},\pi_{2},\pi_{3},\pi_{4})\in\Pi}\sum_{q\in[\ell_{1}],r\in[\ell_{2}]}\frac{\phi_{\pi_{1}\pi_{2}}(qr)\phi_{\pi_{3}\pi_{4}}(qr)}{\eta_{q}\upsilon_{r}}.$

It is clear that the above kernel is equivalent to $h_{ijkl}$ when $\bm{\eta}$
and $\bm{\upsilon}$ are vectors with each entry equal to one. Similarly, when
$\bm{\eta}$ and $\bm{\upsilon}$ are defined with $1+a_{q}$ and
$1+a_{r}^{\prime}$, respectively, then the above kernel corresponds to
$h_{ijkl}^{\bm{a}}$. The U-statistic based on
$h_{ijkl}^{\bm{\eta},\bm{\upsilon}}$ is denoted by
$U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$, which is similarly computed as in (2).
For $q\in[\ell_{1}]$ and $r\in[\ell_{2}]$, define
$o_{qr}=\sum_{i=1}^{n}\mathds{1}(X_{i}=q)\mathds{1}(Y_{i}=r)$,
$o_{q+}=\sum_{r=1}^{\ell_{2}}o_{qr}$ and $o_{+r}=\sum_{q=1}^{\ell_{1}}o_{qr}$.
With this notation in place, we give an alternative expression of
$U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ as follows.

######  Proposition 1 (Alternative expression).

The U-statistic $U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ can be written as

$\displaystyle U^{\bm{\eta},\bm{\upsilon}}_{W}(D)\leavevmode\nobreak\
=\leavevmode\nobreak\
\frac{1}{\sigma(\sigma-3)}\Biggl{[}A_{1}+\frac{1}{(\sigma-1)(\sigma-2)}A_{2}-\frac{2}{\sigma-2}A_{3}\Biggr{]}$

where

$\displaystyle
A_{1}=\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\Biggl{(}\frac{o_{qr}^{2}-o_{qr}}{\eta_{q}\upsilon_{r}}\Biggr{)},\quad
A_{2}=\sum_{q=1}^{\ell_{1}}\Biggl{(}\frac{o_{q+}^{2}-o_{q+}}{\eta_{q}}\Biggr{)}\cdot\sum_{r=1}^{\ell_{2}}\Biggl{(}\frac{o_{+r}^{2}-o_{+r}}{\upsilon_{r}}\Biggr{)},$
$\displaystyle
A_{3}=\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}(o_{q+}o_{+r}-o_{q+}-o_{+r}+1)}{\eta_{q}\upsilon_{r}}.$

We now discuss the average time complexity of
$U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ by assuming that the weight vectors
$\bm{\eta},\bm{\upsilon}$ are given in advance. First of all, we note that the
$\ell_{1}\times\ell_{2}$ contingency table of $D$ is sparse in a sense that it
has at most $\sigma$ non-zero entries. Importantly, the zero entries do not
affect the calculation of $A_{1},A_{2},A_{3}$. Hence we only focus on the non-
zero entries of the contingency table, which can be computed in linear time,
for instance, by using hash tables (e.g., Cormen et al.,, 2022). Similarly,
the non-zero row sums and the non-zero column sums of the contingency table
can be computed in linear time independent of $\ell_{1}$ and $\ell_{2}$. Given
the non-zero entries of
$\\{o_{qr},o_{q+},o_{+r}:q\in[\ell_{1}],r\in[\ell_{2}]\\}$, the computational
complexity of the terms $A_{1},A_{2},A_{3}$ is linear as their non-zero
summands are at most $O(\sigma)$.

For the unweighted U-statistic $U(D)$, there is no additional cost for
computing $\bm{\eta},\bm{\upsilon}$ as they are vectors with each entry equal
to one. For the weighted U-statistics $U_{W}^{\bm{a}}(D)$ and
$U_{W}^{\bm{b}}(D)$, the weight vectors $\bm{\eta},\bm{\upsilon}$ are
functions of $o_{q+}$ and $o_{+r}$ (computed on a separate dataset for
$U_{W}^{\bm{a}}(D)$) and they only require an additional $O(\sigma)$ time to
compute. Thus the overall time complexity is still linear.

##  4 Theoretical results

Having introduced the test statistics, we are now ready to provide the main
theoretical results of this paper. In Section 4.1, we establish the same
sample complexity of the tests using $T$ and $T_{W}$ under multinomial
sampling. In Section 4.2, we provide and analyze more practical tests based on
permutation procedures.

###  4.1 Sample complexity without Poissonization

In this subsection, we revisit two main results of Canonne et al., (2018),
namely Theorem 1.1 and Theorem 1.3 concerning with the sample complexity of a
test using $T$ and $T_{W}$, respectively.

##### Sample complexity of a test based on $T$:

Suppose that the test statistic $T$ is constructed using $N$ i.i.d. samples
from $p_{X,Y,Z}$ where $N\sim\mathrm{Poisson}(n)$. We reject the null of CI
when $T>\zeta\sqrt{\min(n,d)}$ for a sufficiently large constant $\zeta>0$.
Then for $\ell_{1}=\ell_{2}=2$, Theorem 1.1 of Canonne et al., (2018) proves
that the resulting test has the sample complexity

$\displaystyle O\Biggl{(}\max\Bigg{\\{}\frac{d^{1/2}}{\varepsilon^{2}},\
\min\bigg{\\{}\frac{d^{7/8}}{\varepsilon},\
\frac{d^{6/7}}{\varepsilon^{8/7}}\bigg{\\}}\Bigg{\\}}\Biggr{)}.$ (5)

They also prove that this sample complexity is rate optimal by presenting a
matching lower bound.

##### Sample complexity of a test based on $T_{W}$:

The test based on $T$ is not necessarily optimal in the high-dimensional
regime where $\ell_{1}$ and $\ell_{2}$ can vary with other parameters. As
shown in Theorem 1.3 of Canonne et al., (2018), a more general result of the
sample complexity can be derived by using $T_{W}$. In particular, given $N$
i.i.d. samples from $p_{X,Y,Z}$ where $N\sim\mathrm{Poisson}(n)$, we reject
the null of CI when $T_{W}>\zeta^{\prime}\sqrt{\min(n,d)}$ for a sufficiently
large constant $\zeta^{\prime}>0$. The sample complexity of the resulting test
satisfies

$\displaystyle
O\Biggl{(}\max\Biggl{\\{}\min\biggl{\\{}\frac{d^{7/8}\ell_{1}^{1/4}\ell_{2}^{1/4}}{\varepsilon},\
\frac{d^{6/7}\ell_{1}^{2/7}\ell_{2}^{2/7}}{\varepsilon^{8/7}}\biggr{\\}},\
\frac{d^{3/4}\ell_{1}^{1/2}\ell_{2}^{1/2}}{\varepsilon},\
\frac{d^{2/3}\ell_{1}^{2/3}\ell_{2}^{1/3}}{\varepsilon^{4/3}},\
\frac{d^{1/2}\ell_{1}^{1/2}\ell_{2}^{1/2}}{\varepsilon^{2}}\Biggr{\\}}\Biggr{)}.$
(6)

Moreover, this upper bound is shown to be optimal, up to constant factors, in
a number of regimes. See Canonne et al., (2018) for a discussion. We note in
passing that there was an error in the original proof of Theorem 1.3 in
Canonne et al., (2018), which has been recently fixed by Kim, (2022).

We now depoissonize the previous results and establish the same sample
complexity under multinomial sampling.

######  Theorem 1 (Multinomial sampling).

Suppose that we observe
$D_{n}=\\{(X_{1},Y_{1},Z_{1}),\ldots,(X_{n},Y_{n},Z_{n})\\}$ i.i.d. samples
from $p_{X,Y,Z}$ with nonrandom sample size $n$. Then the following two
statements hold:

  1. 1.

Compute $T$ based on $D_{n}$ and reject the null if $T>\zeta\sqrt{\min(n,d)}$
for a sufficiently large constant $\zeta>0$. Then the resulting test has the
sample complexity as in (5) when $\ell_{1},\ell_{2}=O(1)$.

  2. 2.

Compute $T_{W}$ based on $D_{n}$ and reject the null if
$T_{W}>\zeta^{\prime}\sqrt{\min(n,d)}$ for a sufficiently large constant
$\zeta^{\prime}>0$. Then the resulting test has the sample complexity as in
(6).

A few remarks are in order.

######  Remark 3.

  * •

The above theorem essentially says that the tests based on $T$ and $T_{W}$
have the same performance in sample complexity under Poisson sampling and
multinomial sampling. This result may not come as a surprise given that a
Poisson random variable is strongly concentrated around its mean. See
empirical evidence in Kim et al., 2022b . However, the proof under multinomial
sampling turns out to be non-trivial, requiring a careful analysis.

  * •

One of the main technical hurdles in the proof is to overcome a lack of
independence between random variables in different bins when the sample size
is no longer Poisson. The independence property is useful in analyzing the
variance of the sum of U-statistics as it leads to zero covariance terms. We
address the lack of independence by employing the negative association (NA)
property of multinomial random vectors (Joag-Dev and Proschan,, 1983). This NA
property ensures that the covariance terms are non-positive, which turns out
to be enough to ensure the same theoretical guarantees hold under multinomial
sampling.

  * •

Even though we remove Poissonization, the resulting tests are not necessarily
practical. In particular, their critical values depend on unspecified
constants $\zeta$ and $\zeta^{\prime}$. The choice of these constants,
resulting in tight control of the type I error, is challenging in practice. We
take a further step to address this issue via the permutation method in
Section 4.2, and demonstrate their empirical performance in Section 5.

###  4.2 Calibration via permutations

As mentioned earlier, the tests used in Theorem 1 depend on unspecified
constants, which raises the issue of practicality. This section attempts to
address this problem by presenting more practical tests calibrated by the
permutation method, and examine their sample complexity under multinomial
sampling. In particular, we prove that their sample complexity remains the
same as the corresponding (theoretical) tests in Theorem 1. As briefly
mentioned earlier, a similar result was established in Theorem 5 of Kim et
al., 2022b  but under Poisson sampling. In contrast, we do not assume that the
sample size follows a Poisson distribution and therefore reduce the gap
between theory and practice. We start by describing the testing procedures
that we analyze.

##### Permutation test using $T$.

This first test compares the test statistic $T$ with its permutation
correspondences, and rejects the null if the resulting permutation $p$-value
is less than or equal to significance level $\alpha$. To further explain, let
$\Pi_{\sigma_{m}}$ denote the set of all permutations of $[\sigma_{m}]$ for
each $m\in[d]$. Given $\pi$ drawn from $\Pi_{\sigma_{m}}$, we define
$D_{m}^{\pi}$ by rearranging $Y$ values in $D_{m}$ according to $\pi$. More
specifically, suppose that we have
$D_{m}=\\{(X_{1},Y_{1}),\ldots,(X_{\sigma_{m}},Y_{\sigma_{m}})\\}$. Then the
corresponding permuted dataset becomes
$D_{m}^{\pi}=\\{(X_{1},Y_{\pi_{1}}),\ldots,(X_{\sigma_{m}},Y_{\pi_{\sigma_{m}}})\\}$.
Equipped with this notation, we implement Algorithm 2 and make a decision
based on the output.

Algorithm 2 UCI: U-statistic permutation CI test

Input: Sample $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}$, the number of
permutations $B$, significance level $\alpha$ For $j\in[B]$ do For $m\in[d]$
do

  * •

Generate $\pi\sim\mathrm{Uniform}(\Pi_{\sigma_{m}})$ independent of everything
else.

  * •

Compute $U(D_{m}^{\pi})$ as in (2) based on the permuted dataset
$D_{m}^{\pi}$.

End Set $T_{j}\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}U(D_{m}^{\pi})$. End

  * •

Set $T\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq 4)\sigma_{m}U(D_{m})$
computed without permutations.

  * •

Compute the permutation $p$-value

$\displaystyle
p_{\mathrm{perm}}=\frac{1}{B+1}\Biggl{[}\sum_{j=1}^{B}\mathds{1}(T_{j}\geq
T)+1\Biggr{]}.$

Output: Reject $H_{0}$ if $p_{\mathrm{perm}}\leq\alpha$; otherwise, accept
$H_{0}$.

##### Permutation test using $T_{W}$.

The second test that we analyze calculates its $p$-value by comparing $T_{W}$
with its permutation correspondences. The overall procedure is similar to the
previous one except that it utilizes the half-permutation procedure for a
technical reason (see Remark 6 Kim et al., 2022b, ). To explain the procedure,
we decompose $D_{m}$ into $D_{X,m}$, $D_{Y,m}$ and $D_{X,Y,m}$ of size
$t_{1,m}$, $t_{2,m}$ and $2t_{m}+4$, respectively, as in Section 3.2. Unlike
Algorithm 2, we only permute $Y$ values within $D_{X,Y,m}$ for each $m\in[d]$
and then evaluate the significance of $T_{W}$. A more detailed procedure is
described in Algorithm 3.

Algorithm 3 wUCI$\\_$split: weighted U-statistic permutation CI test using
sample splitting

Input: Sample $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}$, the number of
permutations $B$, significance level $\alpha$ For $j\in[B]$ do For $m\in[d]$
do

  * •

Generate $\pi\sim\mathrm{Uniform}(\Pi_{2t_{m}+4})$ independent of everything
else.

  * •

Define $D_{m}^{\pi}=D_{X,m}\cup D_{Y,m}\cup D_{X,Y,m}^{\pi}$.

  * •

Compute $U_{W}^{\bm{a}}(D_{m}^{\pi})$ as in (3) based on the permuted dataset
$D_{m}^{\pi}$.

End Set $T_{j,W}\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{a}}(D_{m}^{\pi})$. End

  * •

Set $T_{W}\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{a}}(D_{m})$ computed without permutations.

  * •

Compute the permutation $p$-value

$\displaystyle
p_{\mathrm{perm}}=\frac{1}{B+1}\Biggl{[}\sum_{j=1}^{B}\mathds{1}(T_{j,W}\geq
T_{W})+1\Biggr{]}.$

Output: Reject $H_{0}$ if $p_{\mathrm{perm}}\leq\alpha$; otherwise, accept
$H_{0}$.

Having outlined the permutation procedures, we now discuss their sample
complexity. First of all, it is noteworthy that both permutation tests are
exact level $\alpha$ in any finite sample scenarios. This simply follows by
the fact that the original test statistic and their permutation
correspondences are exchangeable under the null. Using this observation, it
can be seen that that the resulting $p$-values in Algorithm 2 and 3 are super-
uniform (e.g., Lemma 1 of Romano and Wolf,, 2005). The next theorem turns to
the type II error and establishes their sample complexity.

######  Theorem 2 (Permutation tests).

Suppose that we observe
$D_{n}=\\{(X_{1},Y_{1},Z_{1}),\ldots,(X_{n},Y_{n},Z_{n})\\}$ i.i.d. samples
from $p_{X,Y,Z}$ with nonrandom sample size $n$. We also assume that the
number of random permutations $B$ satisfies
$B\geq\max\\{4(1-\alpha)\alpha^{-1},\,8\alpha^{-2}\log(4\beta^{-1})\\}$ where
$\alpha$ and $\beta$ are pre-specified type I and II errors, respectively.
Then the following two statements hold:

  1. 1.

The test from Algorithm 2 has the sample complexity as in (5) when
$\ell_{1},\ell_{2}=O(1)$.

  2. 2.

The test from Algorithm 3 has the sample complexity as in (6).

We highlight that the above theorem studies the random permutation tests where
$B$ is not required to increase with the sample size $n$. This is in contrast
to full permutation tests considered in Kim et al., 2022b  that enumerate all
possible permutations. We also note that the constant factors in the condition
on $B$ is not tight and can be improved with more effort.

We next turn our attention to $\chi^{2}$\- and $G$-tests and discuss their
sub-optimality.

###  4.3 Sub-optimality of $\chi^{2}$\- and $G$-tests

Practitioners frequently use $\chi^{2}$\- and $G$-tests for conditional
independence, which have nice asymptotic properties in fixed dimensional
settings. In this subsection, we move beyond this fixed dimensional case and
prove that these classical tests are markedly sub-optimal in terms of sample
complexity. To define the $\chi^{2}$-test and $G$-test formally, let us write
$o_{qrs}=\sum_{i=1}^{n}\mathds{1}(X_{i}=q)\mathds{1}(Y_{i}=r)\mathds{1}(Z_{i}=s)$
and $e_{qrs}=o_{q+s}o_{+rs}/o_{++s}$ where
$o_{q+s}=\sum_{r\in[\ell_{2}]}o_{qrs}$, $o_{+rs}=\sum_{q\in[\ell_{1}]}o_{qrs}$
and $o_{++s}=\sum_{q\in[\ell_{1}],r\in[\ell_{2}]}o_{qrs}$, respectively, for
$q\in[\ell_{1}],r\in[\ell_{2}],s\in[d]$. Given this notation, the
$\chi^{2}$-test and $G$-test are based on the following test statistics

$\displaystyle\chi^{2}=\sum_{q\in[\ell_{1}],r\in[\ell_{2}],s\in[d]}\frac{(o_{qrs}-e_{qrs})^{2}}{e_{qrs}}\quad\text{and}\quad
G=2\sum_{q\in[\ell_{1}],r\in[\ell_{2}],s\in[d]}o_{qrs}\log\biggl{(}\frac{o_{qrs}}{e_{qrs}}\biggr{)}.$
(7)

These test statistics converge to a $\chi^{2}$ distribution with
$(\ell_{1}-1)\times(\ell_{2}-1)\times d$ degrees of freedom under the null of
conditional independence and under some regularity conditions (Tsamardinos and
Borboudakis,, 2010). Based on this asymptotic result, $\chi^{2}$-test and
$G$-test reject the null when $\chi^{2}$ and $G$ are larger than the
$1-\alpha$ quantile of the $\chi^{2}$ distribution with
$(\ell_{1}-1)\times(\ell_{2}-1)\times d$ degrees of freedom. We first
emphasize that these classical tests do not control the type I error uniformly
over the null distributions and their validity guarantee requires that the
sample size go to infinity. This is even true for the simplest case where
$d=1$, which corresponds to the unconditional independence problem (see
Berrett and Samworth,, 2021, for details). Moreover, their asymptotic power
can be exactly equal to zero in some regimes where the sample size is much
larger than the bound in (5) as shown below.

######  Proposition 2 (Sub-optimality of $\chi^{2}$\- and $G$-tests).

Assume that $\ell_{1}=\ell_{2}=2$, $\varepsilon=0.25$ and $\alpha\in(0,1)$ is
some fixed constant. Further assume that $d=n\times r_{n}$ where $r_{n}$ is an
arbitrary positive sequence that increases to infinity as
$n\rightarrow\infty$. In this scenario, the worse case power of $\chi^{2}$\-
and $G$-tests approach zero as

$\displaystyle\lim_{n\rightarrow\infty}\inf_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(\chi^{2}>q_{1-\alpha,d})=0\quad\text{and}\quad\lim_{n\rightarrow\infty}\inf_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(G>q_{1-\alpha,d})=0,$

where $q_{1-\alpha,d}$ denotes the $1-\alpha$ quantile of the $\chi^{2}$
distribution with $d$ degrees of freedom.

We provide some remarks on this result.

######  Remark 4.

  * •

As shown in Theorem 2, the proposed tests can achieve significant power
(indeed rate optimal when $\ell_{1}=\ell_{2}=2$) under the same scenario and
$r_{n}\lesssim n^{1/6}$. On the other hand, $\chi^{2}$\- and $G$-tests have
asymptotically zero power for any $r_{n}$ that increases to infinity, which
highlights sub-optimality of these classical tests. Moreover, the UCI and wUCI
tests are valid over the entire class of null distributions unlike asymptotic
$\chi^{2}$\- and $G$-tests.

  * •

At a high-level, the reason behind this negative result is that the critical
values of $\chi^{2}$\- and $G$-tests are not adaptive to the underlying
distribution. More specifically, we can think of a setting where most of
conditional bins are empty with high probability, i.e., the intrinsic
dimension of $Z$ is much smaller than $d$. In this case, it is more natural to
use a critical value that reflects the intrinsic dimension rather than the
ambient dimension. However, the $\chi^{2}$\- and $G$-tests do not take this
intuition into account, and their test statistics can be much smaller than
$q_{1-\alpha,d}$ under the alternative. This leads to asymptotically zero
power as we formally prove in Appendix C.4.

  * •

This issue can be alleviated by using the permutation approach where empty
bins are ignored in calibration (Tsamardinos and Borboudakis,, 2010).
Nevertheless, it is unknown whether the permutation-based $\chi^{2}$\- and
$G$-tests are optimal in terms of sample complexity. We leave this important
question for future work.

##  5 Numerical analysis

In this section, we provide numerical results that compare the proposed tests
(UCI-test in Algorithm 2 and wUCI-test in Algorithm 4) with $\chi^{2}$-test
and $G$-test under various scenarios. For a fair comparison, we calibrate both
$\chi^{2}$-test and $G$-test using the permutation method (as in Algorithm 2)
and reject the null when their permutation $p$-values are less than or equal
to the significance level $\alpha$. Throughout our simulations, we set
$\alpha=0.05$ and the number of permutations $B=199$. All the power values
reported in this section are estimated by Monte Carlo simulation with 10,000
repetitions.

###  5.1 Simulated data examples

We start by comparing the power of the considered tests based on synthetic
datasets. We only focus on the power results given that all of the tests are
calibrated by the permutation procedure, resulting in valid type I error
control in any finite sample sizes. There are eight different scenarios that
we consider under the alternative where the domain sizes of $X,Y,Z$ are set to
$\ell_{1}=\ell_{2}=20$ and $d=10$, respectively. The considered scenarios are
described as follows.

  * •

Scenario 1. Set $p_{Z}$ to be uniform over $[d]$. For each $z\in[d]$, (i)
first let $p_{X,Y|Z}(x,y\,|\,z)\propto x^{-2}y^{-2}$, (ii) then replace
$p_{X,Y|Z}(\ell_{1},\ell_{2}\,|\,z)$ with $0.015$, and (iii) finally normalize
$p_{X,Y|Z}$ to have its sum to be one. This setting results in a strong signal
in $\chi^{2}$ divergence but relatively weaker signal in the $L_{2}$ distance
over bins.

  * •

Scenario 2. Set $p_{Z}$ to be uniform over $[d]$. For each $z\in[d]$, (i) let
$p_{X,Y|Z}(x,y\,|\,z)\propto x^{-2}y^{-2}$, (ii) set
$\delta=\min\\{p_{X,Y|Z}(x,y\,|\,z):x\in[2],y\in[2]\\}$, and (iii) perturb
$p_{X,Y|Z}$ by replacing $p_{X,Y|Z}(x,y\,|\,z)$ with
$p_{X,Y|Z}(x,y\,|\,z)+(-1)^{x+y}\delta$ for $x\in[2]$ and $y\in[2]$. The
resulting probability vector has a small signal in $\chi^{2}$ divergence, but
relatively stronger signal in the $L_{2}$ distance over bins.

  * •

Scenario 3. Set $p_{Z}$ to be uniform over $[d]$. For each $z\in[d]$, (i) set
$p_{X,Y|Z}(x,y\,|\,z)=0$ for all $x\in[\ell_{1}]$ and $y\in[\ell_{2}]$, (ii)
set $p_{X,Y|Z}(1,1\,|\,z)=(1-q)^{2}$,
$p_{X,Y|Z}(1,y\,|\,z)=(1-q)q(\ell_{1}-1)^{-1}$ for
$y\in[\ell_{2}]\\!\setminus\\!\\{1\\}$,
$p_{X,Y|Z}(x,1\,|\,z)=(1-q)q(\ell_{1}-1)^{-1}$ for
$x\in[\ell_{1}]\\!\setminus\\!\\{1\\}$,
$p_{X,Y|Z}(x,y\,|\,z)=q^{2}(\ell_{1}-1)^{-1}$ for
$x=y\in[\ell_{1}]\\!\setminus\\!\\{1\\}$ where $q=0.2$. This simulation
setting is borrowed from Zhang and Zhang, (2022).

  * •

Scenario 4. Set $p_{Z}$ to be uniform over $[d]$. For each $z\in[d]$, let
$p_{X,Y|Z}(x,y\,|\,z)=\\{1+(-1)^{x+y}\\}\ell_{1}^{-1}\ell_{2}^{-1}$ be a
perturbed uniform distribution. This is the setting where $\chi^{2}$-test,
$G$-test and UCI-test perform similarly for unconditional independence
testing. See Figure 5 of Berrett and Samworth, (2021).

  * •

Scenario 5. Set $p_{Z}$ to be uniform over $[d]$. For $z=1$, let
$p_{X,Y|Z}(x,y\,|\,z)=0.25$ for $x\in[2],y\in[2]$ and zero otherwise. In other
words, there is no signal in the first category of $Z$. On the other hand, for
$z\in[d]\\!\setminus\\!\\{1\\}$, set
$p_{X,Y|Z}(x,y\,|\,z)=\\{1+(-1)^{x+y}\\}\ell_{1}^{-1}\ell_{2}^{-1}$ as in
Scenario 4.

  * •

Scenario 6. Set $p_{Z}$ to be uniform over $[d]$. For $z=1$,
$p_{X,Y|Z}(1,1\,|\,z)=p_{X,Y|Z}(2,2\,|\,z)=0.4$,
$p_{X,Y|Z}(1,2\,|\,z)=p_{X,Y|Z}(2,1\,|\,z)=0.1$ and zero otherwise. On the
other hand, for $z\in[d]\\!\setminus\\!\\{1\\}$, set
$p_{X,Y|Z}(x,y\,|\,z)=\ell_{1}^{-1}\ell_{2}^{-1}$, i.e.,
$X\,\rotatebox[origin={c}]{90.0}{$\models$}\,Y$ for
$z\in[d]\\!\setminus\\!\\{1\\}$, resulting in a sparse alternative.

  * •

Scenario 7. Set $p_{Z}(z)\propto z^{-1}$. For $z\in[d]$, let
$p_{X,Y|Z}(x,y\,|\,z)=\\{1+(-1)^{x+y}z^{-1}\\}\ell_{1}^{-1}\ell_{2}^{-1}$. By
construction, the signal becomes weaker as $z$ increases and the sample size
$\sigma_{z}$ tends to be smaller as $z$ increases.

  * •

Scenario 8. $p_{Z}(z)\propto z^{-1}$. For $z\in[d]$, let
$p_{X,Y|Z}(x,y\,|\,z)=\\{1+(-1)^{x+y}(d-z+1)^{-1}\\}\ell_{1}^{-1}\ell_{2}^{-1}$.
Note that the signal becomes stronger as $z$ increases, and the sample size
$\sigma_{z}$ tends to be smaller as $z$ increases. To put it in another way,
this is the reverse setting of Scenario 7.

|
---|---
|
Figure 2: Power comparisons of the considered tests in Scenario 1–Scenario 4 described in Section 5.1. |
---|---
|
Figure 3: Power comparisons of the considered tests in Scenario 5–Scenario 8
described in Section 5.1.

The results are presented in Figure 2 and Figure 3. It is clear from the
results that no test outperforms the others over all scenarios. In particular,
$\chi^{2}$-test has the highest power when the underlying distribution has a
strong signal in $\chi^{2}$ divergence such as Scenario 1, and similarly, UCI-
test performs well when there is a strong signal in the $L_{2}$ distance such
as Scenario 2. It is also interesting to observe that wUCI-test shows an
impressive performance compared to the others in Scenario 3, and all of the
tests perform similarly in Scenario 4 (except $\chi^{2}$-test with large
sample sizes) where the corresponding null distribution of $(X,Y)\,|\,Z$ is
uniform over bins. In Scenario 5 and Scenario 6, we are essentially in a
situation where $\ell_{1}=\ell_{2}=2$ for $z=1$ and $\ell_{1}=\ell_{2}=20$ for
$z\in[d]\\!\setminus\\!\\{1\\}$. In these scenarios, the first bin of $Z$
plays a less important role than the other bins in determining the value of
$\chi^{2}$ and $G$ statistics (these statistics have higher variance when
$\ell_{1}$ and $\ell_{2}$ are large). In contrast, the test statistic for the
UCI-test is mostly dominated by the data from the first bin of $Z$ in the same
scenarios. This explains the outstanding performance of $\chi^{2}$-test and
$G$-test in Scenario 5 where signals are spread out over bins except the first
one. On the other hand, $\chi^{2}$-test and $G$-test have low power in
Scenario 6 where only the first bin of $Z$ has a signal. Under the same
scenarios, UCI-test behaves in the opposite way, attaining high power when the
first bin is significant. Another interesting observation is that wUCI-test
performs as powerful as $G$-test in Scenario 5 whereas it outperforms both
$\chi^{2}$\- and $G$-tests in Scenario 6 when the sample size is large. This
may be explained by the choice of weights in its statistic that roughly
interpolate $\chi^{2}$ weights and uniform weights as explained in Remark 2.
We also note that the test statistic for UCI-test is a linear combination of
U-statistics weighted by sample sizes over bins. This explains the relatively
lower power of UCI-test than $\chi^{2}$\- and $G$-tests in Scenario 7 where
the bins with smaller sample sizes tend to have stronger signals. In contrast,
we observe the opposite behavior in Scenario 8. On the other hand, wUCI
performs the second best in both Scenario 7 and Scenario 8.

To summarize, we observe that different tests perform better than the others
under different scenarios. The proposed tests often dominate the classical
ones when there are strong signals, especially in the $L_{2}$ distance, over
bins with large sample sizes. On the other hand, it is possible to design
situations such as Scenario 1 where the proposed tests attain lower power than
the classical ones. Nevertheless, the classical tests, especially
$\chi^{2}$-test, can fail badly in terms of the worse-case performance. By
contrast, our proposals demonstrate robust performance across different
scenarios, indicating that they can work as practical tools that complement
classical $\chi^{2}$-test and $G$-test under various scenarios.

###  5.2 Real-world data examples

We next provide numerical illustrations based on real-world datasets.

Admission dataset. The first dataset that we look at is the Berkeley
admissions dataset, which is a well-known example of Simpson’s paradox (Bickel
et al.,, 1975). As summarized in Table 1, the dataset consists of 4,526
applications with 3 variables $(X,Y,Z)$ where $X$ and $Y$ are binary
variables, representing the gender (male or female) and the admission status
(admitted or rejected), respectively. The conditional variable $Z$ takes the
department name among $\\{$A, B, C, D, E, F$\\}$. When the dataset is
aggregated over the departments, it appears that male applicants are more
likely to be admitted than woman applicants. However, as reported by Bickel et
al., (1975), there seems to exist a bias in favor of women when looking at the
individual departments, indicating the existence of conditional dependence. We
assess this claim of conditional dependence by implementing the considered
permutation tests. For this dataset, the corresponding $p$-values are computed
as 0.005 for both $\chi^{2}$\- and $G$-tests, 0.04 for UCI-test and 0.03 for
wUCI-test, respectively. All the $p$-values are significant at level
$\alpha=0.05$, revealing evidence of conditional dependence.

Table 1: Admissions data at University of California, Berkeley from the six largest departments in 1973 | Men | Women
---|---|---
Major | Applicants | Admitted | Applicants | Admitted
A | 825 | 62 | 108 | 82
B | 560 | 63 | 25 | 68
C | 325 | 37 | 593 | 34
D | 417 | 33 | 375 | 35
E | 191 | 28 | 393 | 24
F | 373 | 6 | 341 | 7

Diamonds dataset. Next we consider the diamonds dataset available in R package
ggplot2. The dataset contains the information of 53,940 diamonds including
their price, clarity, color, quality of the cut, etc. In our analysis, the
price variable is partitioned into 100 intervals of equal size. We set the
corresponding categorized price variable as $X$ and set the clarity variable
as $Y$. The clarity variable has 8 categories (I1, SI2, SI1, VS2, VS1, VVS2,
VVS1, IF) and it measures the purity of a diamond. The conditional variable
$Z$ is chosen to be either the cut variable or the color variable in our
analysis. Both variables are discrete with 5 (Fair, Good, Very Good, Premium,
Ideal) and 7 (D, E, F, G, H, I, J) categories, respectively. In the
experiments, we treat the entire dataset as the population (thereby the ground
truth is known to us)444We verified numerically that $X$ and $Y$ are
conditionally dependent on $Z$ by comparing the distribution of $p_{X,Y,Z}$
and $p_{X|Z}p_{Y|Z}p_{Z}$. from which we randomly draw $n$ observations
without replacement. Based on this subsample of size $n$, we compute the
permutation $p$-values of the considered tests and we repeat this process
$10,000$ times to estimate their power. The results can be found in Figure 4
where we collect the power of
$\\{\texttt{UCI},\texttt{wUCI},\chi^{2},G\\}$-tests by changing the sample
size $n$. The left panel of Figure 4 provides the power results when the
conditional variable is set to be the color variable. As can be seen,
$\chi^{2}$-test has the significantly lower power than the others. Among the
other three tests, wUCI-test has the highest power followed by UCI-test while
the difference is minor. We can see a similar pattern from the right panel of
Figure 4 where the conditional variable is set to be the cut variable. These
results highlight the practical value of the proposed tests in analyzing real-
world datasets where classical tests potentially suffer from low power.

|
---|---
Figure 4: Power comparisons of the considered tests based on the diamonds
dataset. Both panels analyze independence between the (categorized) price and
clarity variables conditional on the color variable and the cut variable,
respectively. All of the tests have increasing power as the sample size
increases. Markedly, $\chi^{2}$-test has significantly lower power than the
other tests, whereas wUCI-test seems to perform the best for this dataset.

##  6 Discussion

In this paper, we have revisited recent developments of CI testing for
discrete data. Despite attractive theoretical properties, these recent tests
have limited practical value, relying on Poissonization and unspecified
constants in their critical values. In this work, we have made an attempt to
bridge the gap between theory and practice by removing Poissonization and
utilizing the Monte Carlo permutation method to calibrate test statistics. We
have also complemented our theoretical results with a thorough numerical
analysis and demonstrated certain benefits of the proposed tests over
classical $\chi^{2}$\- and $G$-tests. Finally, we have developed R package UCI
that implements the proposed methods.

Our work leaves several important avenues for future research. One prominent
direction is to depoissonize other sample complexity results in the literature
using the tools developed in this paper. For instance, one can reproduce the
results of Neykov et al., (2021); Kim et al., 2022b  for continuous CI testing
without Poissonization. Another direction which may be fruitful to pursue is
to devise a CI test that incorporates prior information about potential
alternative distributions. For example, suppose that we are in an alternative
setting where only a handful of conditional categories are significant. In
this case, it is possible to obtain a substantial power gain by using sparse
weights in the proposed statistics, and one could analyze the resulting tests.
Additionally, it would be interesting to see whether wUCI-test or other tests
without sample splitting can achieve the sample complexity (6) in general. We
leave these interesting questions for future work.

##### Acknowledgements.

This work was partially supported by funding from the NSF grants DMS-2113684
and DMS-2310632, as well as an Amazon AI and a Google Research Scholar Award
to SB. MN acknowledges support from the NSF grant DMS-2113684.

## References

  * Agresti, (1992) Agresti, A. (1992).  A Survey of Exact Inference for Contingency Tables.  Statistical Science, 7(1):131–153.
  * Agresti, (2003) Agresti, A. (2003).  Categorical Data Analysis.  John Wiley & Sons.
  * Balakrishnan and Wasserman, (2018) Balakrishnan, S. and Wasserman, L. (2018).  Hypothesis testing for high-dimensional multinomials: A selective review.  The Annals of Applied Statistics, 12(2):727–749.
  * Balakrishnan and Wasserman, (2019) Balakrishnan, S. and Wasserman, L. (2019).  Hypothesis testing for densities and high-dimensional multinomials: Sharp local minimax rates.  The Annals of Statistics, 47(4):1893–1927.
  * Berrett et al., (2021) Berrett, T. B., Kontoyiannis, I., and Samworth, R. J. (2021).  Optimal rates for independence testing via U-statistic permutation tests.  The Annals of Statistics, 49(5):2457–2490.
  * Berrett and Samworth, (2021) Berrett, T. B. and Samworth, R. J. (2021).  USP: an independence test that improves on Pearson’s chi-squared and the G-test.  Proceedings of the Royal Society A, 477(2256):20210549.
  * Berrett et al., (2020) Berrett, T. B., Wang, Y., Barber, R. F., and Samworth, R. J. (2020).  The conditional permutation test for independence while controlling for confounders.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):175–197.
  * Bickel et al., (1975) Bickel, P. J., Hammel, E. A., and O’Connell, J. W. (1975).  Sex Bias in Graduate Admissions: Data from Berkeley: Measuring bias is harder than is usually assumed, and the evidence is sometimes contrary to expectation.  Science, 187(4175):398–404.
  * Bishop et al., (2007) Bishop, Y. M., Fienberg, S. E., and Holland, P. W. (2007).  Discrete Multivariate Analysis: Theory and Practice.  Springer Science & Business Media.
  * Boucheron et al., (2013) Boucheron, S., Lugosi, G., and Massart, P. (2013).  Concentration inequalities: A nonasymptotic theory of independence.  Oxford university press.
  * Candés et al., (2018) Candés, E., Fan, Y., Janson, L., and Lv, J. (2018).  Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3):551–577.
  * Canonne et al., (2018) Canonne, C. L., Diakonikolas, I., Kane, D. M., and Stewart, A. (2018).  Testing Conditional Independence of Discrete Distributions.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2018, page 735–748, New York, NY, USA. Association for Computing Machinery.
  * Canonne et al., (2022) Canonne, C. L. et al. (2022).  Topics and Techniques in Distribution Testing: A Biased but Representative Sample.  Foundations and Trends® in Communications and Information Theory, 19(6):1032–1198.
  * Chan et al., (2014) Chan, S.-O., Diakonikolas, I., Valiant, P., and Valiant, G. (2014).  Optimal algorithms for testing closeness of discrete distributions.  In Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms, pages 1193–1203. SIAM.
  * Chatterjee, (2022) Chatterjee, S. (2022).  A survey of some recent developments in measures of association.  arXiv preprint arXiv:2211.04702.
  * Cormen et al., (2022) Cormen, T. H., Leiserson, C. E., Rivest, R. L., and Stein, C. (2022).  Introduction to Algorithms.  MIT press.
  * De Campos and Huete, (2000) De Campos, L. M. and Huete, J. F. (2000).  A new approach for learning belief networks using independence criteria.  International Journal of Approximate Reasoning, 24(1):11–37.
  * Diakonikolas and Kane, (2016) Diakonikolas, I. and Kane, D. M. (2016).  A New Approach for Testing Properties of Discrete Distributions.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 685–694. IEEE.
  * Gretton et al., (2005) Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005).  Measuring Statistical Dependence with Hilbert-Schmidt norms.  In International Conference on Algorithmic Learning Theory, pages 63–77. Springer.
  * Haberman, (1988) Haberman, S. J. (1988).  A Warning on the Use of Chi-squared Statistics with Frequency Tables with Small Expected Cell Counts.  Journal of the American Statistical Association, 83(402):555–560.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal Inference in Statistics, Social, and Biomedical Sciences.  Cambridge University Press.
  * Jacquet and Szpankowski, (1998) Jacquet, P. and Szpankowski, W. (1998).  Analytical Depoissonization and its Applications.  Theoretical Computer Science, 201(1-2):1–62.
  * Joag-Dev and Proschan, (1983) Joag-Dev, K. and Proschan, F. (1983).  Negative association of random variables with applications.  The Annals of Statistics, 11(1):286–295.
  * Kac, (1949) Kac, M. (1949).  On deviations between theoretical and empirical distributions.  Proceedings of the National Academy of Sciences, 35(5):252–257.
  * Kim, (2020) Kim, I. (2020).  Multinomial goodness-of-fit based on U-statistics: High-dimensional asymptotic and minimax optimality.  Journal of Statistical Planning and Inference, 205:74–91.
  * Kim, (2022) Kim, I. (2022).  Comments on “Testing Conditional Independence of Discrete Distributions”.  arXiv preprint arXiv:2207.02819.
  * (27) Kim, I., Balakrishnan, S., and Wasserman, L. (2022a).  Minimax optimality of permutation tests.  The Annals of Statistics, 50(1):225–251.
  * (28) Kim, I., Neykov, M., Balakrishnan, S., and Wasserman, L. (2022b).  Local permutation tests for conditional independence.  The Annals of Statistics, 50(6):3388–3414.
  * Koller and Friedman, (2009) Koller, D. and Friedman, N. (2009).  Probabilistic graphical models: principles and techniques.  MIT press.
  * Li and Fan, (2020) Li, C. and Fan, X. (2020).  On nonparametric conditional independence tests for continuous variables.  Wiley Interdisciplinary Reviews: Computational Statistics, 12(3):e1489.
  * Liu et al., (2022) Liu, M., Katsevich, E., Janson, L., and Ramdas, A. (2022).  Fast and powerful conditional randomization testing via distillation.  Biometrika, 109(2):277–293.
  * Lundborg et al., (2022) Lundborg, A. R., Shah, R. D., and Peters, J. (2022).  Conditional independence testing in Hilbert spaces with applications to functional data analysis.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(5):1821–1850.
  * McDonald, (2014) McDonald, J. H. (2014).  G-test of goodness-of-fit.  Handbook of Biological Statistics, 487:53–58.
  * Neykov et al., (2021) Neykov, M., Balakrishnan, S., and Wasserman, L. (2021).  Minimax optimal conditional independence testing.  The Annals of Statistics, 49(4):2151–2177.
  * Pearl, (2014) Pearl, J. (2014).  Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference.  Elsevier.
  * Pearson, (1900) Pearson, K. (1900).  On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 50(302):157–175.
  * Penrose, (2003) Penrose, M. (2003).  Random Geometric Graphs, volume 5.  OUP Oxford.
  * Petersen and Hansen, (2021) Petersen, L. and Hansen, N. R. (2021).  Testing Conditional Independence via Quantile Regression Based Partial Copulas.  Journal of Machine Learning Research, 22(70):1–47.
  * Romano and Wolf, (2005) Romano, J. P. and Wolf, M. (2005).  Exact and approximate stepdown methods for multiple hypothesis testing.  Journal of the American Statistical Association, 100(469):94–108.
  * Shah and Peters, (2020) Shah, R. D. and Peters, J. (2020).  The hardness of conditional independence testing and the generalised covariance measure.  The Annals of Statistics, 48(3):1514–1538.
  * Song et al., (2012) Song, L., Smola, A., Gretton, A., Bedo, J., and Borgwardt, K. (2012).  Feature Selection via Dependence Maximization.  Journal of Machine Learning Research, 13(5):1393–1434.
  * Spohn, (1994) Spohn, W. (1994).  On the Properties of Conditional Independence.  In Patrick suppes: Scientific philosopher, pages 173–196. Springer.
  * Tsamardinos and Borboudakis, (2010) Tsamardinos, I. and Borboudakis, G. (2010).  Permutation Testing Improves Bayesian Network Learning.  In Joint European conference on machine learning and knowledge discovery in databases, pages 322–337. Springer.
  * Vaart and Wellner, (1996) Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence.  In Weak convergence and empirical processes, pages 16–28. Springer.
  * Valiant and Valiant, (2017) Valiant, G. and Valiant, P. (2017).  An Automatic Inequality Prover and Instance Optimal Identity Testing.  SIAM Journal on Computing, 46(1):429–455.
  * Valiant, (2011) Valiant, P. (2011).  Testing Symmetric Properties of Distributions.  SIAM Journal on Computing, 40(6):1927–1968.
  * Yao and Tritchler, (1993) Yao, Q. and Tritchler, D. (1993).  An exact analysis of conditional independence in several $2\times 2$ contingency tables.  Biometrics, pages 233–236.
  * Zhang and Zhang, (2022) Zhang, J. and Zhang, Z. (2022).  A Normal Test for Independence via Generalized Mutual Information.  arXiv preprint arXiv:2207.09541.

##  Appendix A Algorithm

The below is the algorithm for wUCI-test used in our simulation study.

Algorithm 4 wUCI: weighted U-statistic permutation CI test

Input: Sample $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}$, the number of
permutations $B$, significance level $\alpha$ For $j\in[B]$ do For $m\in[d]$
do

  * •

Generate $\pi\sim\mathrm{Uniform}(\Pi_{\sigma_{m}})$ independent of everything
else.

  * •

Compute $U_{W}^{\bm{b}}(D_{m}^{\pi})$ as in (4) based on the permuted dataset
$D_{m}^{\pi}$.

End Set $T_{W,j}^{\dagger}\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{b}}(D_{m}^{\pi})$. End

  * •

Set $T_{W}^{\dagger}\leftarrow\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\omega_{m}U_{W}^{\bm{b}}(D_{m})$ computed without permutations.

  * •

Compute the permutation $p$-value

$\displaystyle
p_{\mathrm{perm}}=\frac{1}{B+1}\Biggl{[}\sum_{j=1}^{B}\mathds{1}(T_{W,j}^{\dagger}\geq
T_{W}^{\dagger})+1\Biggr{]}.$

Output: Reject $H_{0}$ if $p_{\mathrm{perm}}\leq\alpha$; otherwise, accept
$H_{0}$.

##  Appendix B Lemmas

In this section, we prove several lemmas in Canonne et al., (2018) under
multinomial sampling. These results are main building blocks that lead to
Theorem 1 and Theorem 2 in the main text. We start by showing that the result
in Lemma 3.1 of Canonne et al., (2018) holds without Poissonization.

######  Lemma 1.

Suppose that $X\sim\mathrm{Binomial}(n,p)$ where $n\geq 4$. Then there exists
an absolute constant $\gamma>0.85$ such that

$\displaystyle\mathbb{E}[X\mathds{1}(X\geq
4)]\geq\gamma\min\\{np,(np)^{4}\\}.$

The second lemma that we prove corresponds to Claim 2.1 of Canonne et al.,
(2018).

######  Lemma 2.

Suppose that $X\sim\mathrm{Binomial}(n,p)$ where $n\geq 4$. Then there exists
a constant $\gamma>0$ such that

$\displaystyle\mathrm{Var}[X\mathds{1}(X\geq
4)]\leq\gamma\mathbb{E}[X\mathds{1}(X\geq 4)].$

The following lemma proves Claim 2.2 of Canonne et al., (2018) for a binomial
random variable with parameters $n,p$. It is worth pointing out that the the
original statement of Claim 2.2 in Canonne et al., (2018) contains an error,
which has been corrected by Kim, (2022).

######  Lemma 3.

Suppose that $X\sim\mathrm{Binomial}(n,p)$ where $n\geq 9$ and let $a,b\geq
2$. Then there exists a constant $C>0$ such that

$\displaystyle\mathrm{Var}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
C\min\\{np,\sqrt{ab}\\}\mathbb{E}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}.$

The following lemma proves Claim 2.3 of Canonne et al., (2018) for a binomial
random variable with parameters $n,p$.

######  Lemma 4.

Suppose that $X\sim\mathrm{Binomial}(n,p)$ where $n\geq 9$. Then there exists
a constant $C>0$ such that for any $a,b\geq 2$ and $\lambda=np$,

$\displaystyle\mathbb{E}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}\geq
C\min\bigl{\\{}\lambda\sqrt{\min(\lambda,a)\min(\lambda,b)},\lambda^{4}\bigr{\\}}.$

We say that random variables $X_{1},\ldots,X_{k}$ are negatively associated if

$\displaystyle\mathrm{Cov}\\{f_{1}(X_{i},i\in A_{1}),f_{2}(X_{j},j\in
A_{2})\\}\leq 0,$

for all increasing functions $f_{1},f_{2}$ and for every pair of disjoint
subsets $A_{1},A_{2}$ of $\\{1,2,\ldots,k\\}$. It is well-known that a
multinomial sample is negatively associated, which we record below.

######  Lemma 5 (Joag-Dev and Proschan, (1983)).

Let $(\sigma_{1},\ldots,\sigma_{d})$ have a multinomial distribution with
parameters $n$ and $(p_{1},\ldots,p_{d})$. Then for every pair of disjoint
subsets $A_{1},A_{2}$ of $\\{1,\ldots,d\\}$, it holds that

$\displaystyle\mathrm{Cov}\\{f_{1}(\sigma_{i},i\in
A_{1}),f_{2}(\sigma_{j},j\in A_{2})\\}\leq 0,$

whenever $f_{1}$ and $f_{2}$ are increasing.

The next lemma builds on the negative association of a multinomial random
vector, which is proved in Appendix C.9.

######  Lemma 6.

Suppose that $(\sigma_{1},\ldots,\sigma_{d})$ has a multinomial distribution
with parameters $n$ and $(p_{1},\ldots,p_{d})$. Then for any $1\leq i\neq
j\leq d$

$\displaystyle\mathrm{Cov}\bigl{\\{}\sigma_{i}\mathds{1}(\sigma_{i}\geq
4),\sigma_{j}\mathds{1}(\sigma_{j}\geq 4)\bigr{\\}}\leq 0.$

Moreover, let
$\omega_{i}=\sqrt{\min\\{\sigma_{i},\ell_{1}\\}\min\\{\sigma_{i},\ell_{2}\\}}$
for $\ell_{1},\ell_{2}\geq 0$. Then for any $1\leq i\neq j\leq d$

$\displaystyle\mathrm{Cov}\bigl{\\{}\sigma_{i}\mathds{1}(\sigma_{i}\geq
4)\omega_{i},\sigma_{j}\mathds{1}(\sigma_{j}\geq 4)\omega_{j}\bigr{\\}}\leq
0.$

As shown in Lemma 6 of Kim et al., 2022b , rejecting the null when the
permutation $p$-value is less than or equal to $\alpha$ is equivalent to
rejecting the null when the corresponding test statistic is greater than the
$1-\alpha$ quantile of the permutation distribution. We state this result
below for completeness, and its proof can be found in Appendix C.10.

######  Lemma 7 (Quantile).

Let $q_{1-\alpha}$ be the $1-\alpha$ quantile of the empirical distribution of
$V,V_{1},\ldots,V_{B}$. Then it holds that

$\displaystyle\mathds{1}\Biggl{(}\frac{1}{B+1}\Biggl{[}\sum_{i=1}^{B}\mathds{1}(V_{i}\geq
V)+1\Biggr{]}>\alpha\Biggr{)}=\mathds{1}(V\leq q_{1-\alpha}).$

Finally, we state the Efron–Stein inequality for completeness (e.g., Theorem
3.1 of Boucheron et al.,, 2013).

######  Lemma 8 (Efron–Stein inequality).

Let $X_{1},\ldots,X_{n}$ be independent random variables and let $Z=f(\bm{X})$
be a square-integrable function of $\bm{X}=(X_{1},\ldots,X_{n})$. Moreover,
$X_{1}^{\prime},\ldots,X_{n}^{\prime}$ are independent copies of
$X_{1},\ldots,X_{n}$ and define

$\displaystyle
Z^{\prime}_{i}=f(X_{1},\ldots,X_{i-1},X_{i}^{\prime},X_{i+1},\ldots,X_{n}).$

Then

$\displaystyle\mathrm{Var}[Z]\leq\frac{1}{2}\sum_{i=1}^{n}\mathbb{E}\bigl{[}(Z-Z_{i}^{\prime})^{2}\bigr{]}.$

##  Appendix C Proofs

This section collects the proof of the main theorems and lemmas.

###  C.1 Proof of Theorem 1

Given the lemmas in Appendix B together with the results in Canonne et al.,
(2018), the claims in Theorem 1 are almost immediate. We present more details
below. Since the permutation test controls the type I error in any finite
sample sizes, we only focus on the alternative case where the underlying
distribution is $\varepsilon$-far from the null in the $L_{1}$ distance, i.e.,
$\inf_{q\in\mathcal{P}_{0}}\|p-q\|_{1}\geq\varepsilon$, for both $T$ and
$T_{W}$.

####  C.1.1 Claim for $T$

Starting with the first claim, we analyze the expectation and the variance of
$T$ under multinomial sampling. Conditional on $\sigma_{m}$, $U(D_{m})$ is an
unbiased estimator of
$\delta_{m}^{2}:=\|p_{X,Y\,|\,Z=m}-p_{X|Z=m}p_{Y|Z=m}\|_{2}^{2}$. Moreover,
since $\sigma_{m}\sim\mathrm{Binomial}\bigl{(}n,p_{Z}(m)\bigr{)}$, Lemma 1
proves that

$\displaystyle\mathbb{E}[T]\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathbb{E}\bigl{[}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}U(D_{m})\bigr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathbb{E}[\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}]\cdot\delta_{m}^{2}$ $\displaystyle\gtrsim\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\min\bigl{\\{}np_{Z}(m),\bigl{(}np_{Z}(m)\bigr{)}^{4}\bigr{\\}}\cdot\delta_{m}^{2}.$

Set the sample size $n=\beta^{\prime}\max\bigl{(}\sqrt{d}/\varepsilon^{\prime
2},\min(d^{7/8}/\varepsilon^{\prime},d^{6/7}\varepsilon^{\prime 8/7})\bigr{)}$
where $\varepsilon^{\prime}:=\varepsilon/\sqrt{\ell_{1}\ell_{2}}$ and
$\beta^{\prime}\geq 1$ is a sufficiently large constant. Using the inequality,
we follow the same lines of the proof of Proposition 3.1 in Canonne et al.,
(2018) and show that

$\displaystyle\mathbb{E}[T]\leavevmode\nobreak\ \gtrsim\leavevmode\nobreak\
\beta^{\prime}\sqrt{\min(n,d)}.$

Turning to the variance of $T$ and letting
$\bm{\sigma}:=(\sigma_{1},\ldots,\sigma_{d})$, the law of total variance
yields

$\displaystyle\mathrm{Var}[T]\leavevmode\nobreak\ =\leavevmode\nobreak\
\mathbb{E}\bigl{[}\mathrm{Var}[T\,|\,\bm{\sigma}]\bigr{]}+\mathrm{Var}\bigl{[}\mathbb{E}[T\,|\,\bm{\sigma}]\bigr{]}.$

The analysis of the first term is essentially the same as Proposition 3.2 of
Canonne et al., (2018). First note that conditional on $\bm{\sigma}$,
$D_{1},\ldots,D_{d}$ are mutually independent (given
$\sigma_{1},\ldots,\sigma_{d}$, we independently generate
$D_{1},\ldots,D_{d}$), and therefore the following identity holds

$\displaystyle\mathrm{Var}[T\,|\,\bm{\sigma}]\leavevmode\nobreak\
=\leavevmode\nobreak\ \sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}^{2}\mathrm{Var}[U(D_{m})\,|\,\sigma_{m}].$

We then use the variance bound for $U(D_{m})$ in equation (4) of Canonne et
al., (2018) and see

$\displaystyle\mathrm{Var}[T\,|\,\bm{\sigma}]\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}^{2}\biggl{(}\frac{\|p_{X,Y\,|\,Z=m}-p_{X|Z=m}p_{Y|Z=m}\|_{2}^{2}}{\sigma_{m}}+\frac{1}{\sigma_{m}^{2}}\biggr{)}$
$\displaystyle\lesssim\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)\sigma_{m}\|p_{X,Y\,|\,Z=m}-p_{X|Z=m}p_{Y|Z=m}\|_{2}^{2}+\sum_{m\in[d]}\mathds{1}(\sigma_{m}\geq
4)$ $\displaystyle\lesssim\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}[T\,|\,\bm{\sigma}]+\min(n,d).$

This reveals that
$\mathbb{E}\bigl{[}\mathrm{Var}[T\,|\,\bm{\sigma}]\bigr{]}\lesssim\mathbb{E}[T]+\min(n,d)$.
For the second term, the analysis is again similar to that in Proposition 3.2
of Canonne et al., (2018). The only difference is that
$\sigma_{1},\ldots,\sigma_{d}$ are no longer independent but negatively
associated. In Canonne et al., (2018), the independence between
$\sigma_{1},\ldots,\sigma_{d}$ is used to show that

$\displaystyle\mathrm{Var}\bigl{[}\mathbb{E}[T\,|\,\bm{\sigma}]\bigr{]}=\mathrm{Var}\Biggl{[}\sum_{m\in[d]}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\delta_{m}^{2}\Biggr{]}=\sum_{m\in[d]}\delta_{m}^{4}\mathrm{Var}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}.$

In our case, we make use of Lemma 6 (negative association) and prove that

$\displaystyle\mathrm{Var}\Biggl{[}\sum_{m\in[d]}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\delta_{m}^{2}\Biggr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\delta_{m}^{4}\mathrm{Var}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}+\sum_{\begin{subarray}{c}m,m^{\prime}\in[d]\\\ m\neq
m^{\prime}\end{subarray}}\mathrm{Cov}\bigl{\\{}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4),\sigma_{m^{\prime}}\mathds{1}(\sigma_{m^{\prime}}\geq 4)\bigr{\\}}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\delta_{m}^{4}\mathrm{Var}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}.$

This does not change the overall conclusion as our main interest is in
bounding the variance (up to constants). Additionally, Lemma 2 yields that
$\mathrm{Var}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}\lesssim\mathbb{E}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}$ and thus

$\displaystyle\mathrm{Var}\bigl{[}\mathbb{E}[T\,|\,\bm{\sigma}]\bigr{]}\lesssim\sum_{m\in[d]}\delta_{m}^{4}\mathbb{E}\bigl{[}\sigma_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}\lesssim\mathbb{E}[T].$

Combining the previous bounds, we have
$\mathrm{Var}[T]\lesssim\min(n,d)+\mathbb{E}[T]$. Given the bounds on the
expectation and the variance, the claim for $T$ follows by Chebyshev’s
inequality and adjusting the value of $\beta^{\prime}$ and $\zeta$ as in
Section 3.1.3 of Canonne et al., (2018).

####  C.1.2 Claim for $T_{W}$

Next we prove the sample complexity result for $T_{W}$. The type II error
result based on $T_{W}$ under Poisson sampling is given in Lemma 5.6 of
Canonne et al., (2018). This lemma builds upon several preliminary results
proved under Poisson sampling. Once we depoissonize these preliminary results,
the rest of the proof remains the same. Instead of reproducing the entire
proof, we only go over places where Poissonization trick is critical and
reprove the corresponding results under multinomial sampling.

  * •

The first place where Poissonization plays a role is Lemma 5.2 of Canonne et
al., (2018). In the process of investigating the expectation of $T_{W}$,
Canonne et al., (2018) introduce a quantity $D=\sum_{m\in[d]}D_{m}$ where

$\displaystyle
D_{m}=\sigma_{m}\omega_{m}\frac{\varepsilon_{m}^{2}}{\ell_{1}\ell_{2}}\mathds{1}(\sigma_{m}\geq
4)\quad\text{and}\quad\varepsilon_{m}=\mathrm{TV}\bigl{(}p_{X,Y\,|\,Z=m},\,p_{X|Z=m}p_{Y|Z=m}\bigr{)}.$

Then their lemma 5.2 claims that

$\displaystyle\mathbb{E}[D]\leavevmode\nobreak\ \gtrsim\leavevmode\nobreak\
\sum_{m\in[d]}\frac{\varepsilon_{m}^{2}}{\ell_{1}\ell_{2}}\min(\alpha_{m}\beta_{m},\alpha_{m}^{4}),$

where $\alpha_{m}=np_{Z}(m)$ and
$\beta_{m}=\sqrt{\min(\alpha_{m},\ell_{1})\min(\alpha_{m},\ell_{2})}$. We note
that the same result holds under multinomial sampling since

$\displaystyle\mathbb{E}[D]\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\frac{\varepsilon_{m}^{2}}{\ell_{1}\ell_{2}}\mathbb{E}\bigl{[}\sigma_{m}\omega_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}\leavevmode\nobreak\ \gtrsim\leavevmode\nobreak\
\sum_{m\in[d]}\frac{\varepsilon_{m}^{2}}{\ell_{1}\ell_{2}}\min(\alpha_{m}\beta_{m},\alpha_{m}^{4}),$

where the inequality follows by Lemma 4.

  * •

The second place where Poissonization is important is Lemma 5.3 of Canonne et
al., (2018), which claims that

$\displaystyle\mathrm{Var}[D]\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
\mathbb{E}[D].$

This result holds under multinomial sampling as well. To see this, let
$\varepsilon_{m}^{\prime}:=\varepsilon_{m}/\sqrt{\ell_{1}\ell_{2}}$ and
observe

$\displaystyle\mathrm{Var}[D]\leavevmode\nobreak\ =\leavevmode\nobreak\
\mathrm{Var}\biggl{[}\sum_{m\in[d]}\sigma_{m}\omega_{m}\varepsilon_{m}^{\prime
2}\mathds{1}(\sigma_{m}\geq 4)\biggr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathrm{Var}\bigl{[}\sigma_{m}\omega_{m}\varepsilon_{m}^{\prime
2}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}+\sum_{\begin{subarray}{c}m,m^{\prime}\in[d]\\\ m\neq
m^{\prime}\end{subarray}}\mathrm{Cov}\bigl{\\{}\sigma_{m}\omega_{m}\mathds{1}(\sigma_{m}\geq
4),\sigma_{m^{\prime}}\omega_{m^{\prime}}\mathds{1}(\sigma_{m^{\prime}}\geq
4)\bigr{\\}}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\mathrm{Var}\bigl{[}\sigma_{m}\omega_{m}\varepsilon_{m}^{\prime
2}\mathds{1}(\sigma_{m}\geq 4)\bigr{]},$

where the inequality uses the negative association result established in Lemma
6. We then use Lemma 3 and see that

$\displaystyle\sum_{m\in[d]}\mathrm{Var}\bigl{[}\sigma_{m}\omega_{m}\varepsilon_{m}^{\prime
2}\mathds{1}(\sigma_{m}\geq 4)\bigr{]}\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\varepsilon_{m}^{\prime
4}\sqrt{\ell_{1}\ell_{2}}\mathbb{E}\bigl{[}\sigma_{m}\omega_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}$ $\displaystyle\lesssim\leavevmode\nobreak\ $
$\displaystyle\sum_{m\in[d]}\varepsilon_{m}^{\prime
2}\mathbb{E}\bigl{[}\sigma_{m}\omega_{m}\mathds{1}(\sigma_{m}\geq
4)\bigr{]}=\mathbb{E}[D].$

Therefore, Lemma 5.3 of Canonne et al., (2018) holds under multinomial
sampling as well.

  * •

Lemma 5.4 of Canonne et al., (2018) provides a lower bound for the
(conditional) expected value of the test statistic and an upper bound for the
(conditional) variance of the test statistic, which simultaneously hold with
high probability. In this result, the authors use a concentration property of
$N\sim\mathrm{Poisson}(n)$ to its mean in order to connect $\min(N,d)$ with
$\min(n,d)$. This step can be skipped under multinomial sampling as $n$ is
fixed in this case.

Other than Lemma 5.2, Lemma 5.3 and Lemma 5.4 of Canonne et al., (2018), the
other steps of the proof of Lemma 5.6 in Canonne et al., (2018) are
essentially the same under multinomial sampling. Hence we omit the details and
finish the proof.

###  C.2 Proof of Theorem 2

Our proof of Theorem 2 leverages the technical tools developed in Kim et al.,
2022a  and Kim et al., 2022b  for controlling the randomness of permuted
statistics.

####  C.2.1 Claim for $T$

We begin with the first part that proves the sample complexity of UCI-test in
Algorithm 2 based on the test statistic $T$. By Lemma 7, the type II error of
UCI-test is equal to

$\displaystyle\mathbb{P}(p_{\mathrm{perm}}>\alpha)=\mathbb{P}(T\leq
q_{1-\alpha}),$

where $q_{1-\alpha}$ is the $1-\alpha$ quantile of the empirical distribution
of $T,T_{1},\ldots,T_{B}$. Our main strategy is to show that $q_{1-\alpha}$ is
upper bounded by the critical value used in Theorem 1. More specifically,
suppose that there exists a sufficiently large $\zeta>0$ (not necessarily the
same $\zeta$ in Theorem 1) such that

$\displaystyle\mathbb{P}\bigl{(}q_{1-\alpha}\geq\zeta\sqrt{\min(n,d)}\bigr{)}\leq\frac{\beta}{2}.$
(8)

Then by denoting the event of $q_{1-\alpha}<\zeta\sqrt{\min(n,d)}$ by
$\mathcal{A}$, the type II error of UCI-test is bounded by

$\displaystyle\mathbb{P}(p_{\mathrm{perm}}>\alpha)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ $ $\displaystyle\mathbb{P}(T\leq
q_{1-\alpha},\,\mathcal{A})+\mathbb{P}(\mathcal{A}^{c})$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\mathbb{P}\bigl{(}T\leq\zeta\sqrt{\min(n,d)}\bigr{)}+\frac{\beta}{2}$
$\displaystyle\leq\leavevmode\nobreak\ $ $\displaystyle\beta,$

where the last inequality uses the result of Theorem 1. Having this
inequality, we only need to prove inequality (8). To this end, we consider the
permutation distribution of $T$ based on all possible local permutations. To
explain it further, recall that $\Pi_{\sigma_{m}}$ is the set of all
permutations of $[\sigma_{m}]$ for each $m\in[d]$, and $D_{m}^{\pi}$ is the
set of the data with $Z=m$ whose $Y$ values are permuted based on
$\pi\in\Pi_{\sigma_{m}}$. Let $\pi_{(m)}$ be a random permutation drawn
uniformly from $\Pi_{\sigma_{m}}$ for each $m\in[d]$. Further assume that
$\pi_{(1)},\ldots,\pi_{(m)}$ are mutually independent and write
$\bm{\pi}=\\{\pi_{(1)},\ldots,\pi_{(m)}\\}$. Let $T^{\bm{\pi}}$ be the test
statistic computed as $T$ but based on
$\\{D_{1}^{\pi_{(1)}},\ldots,D_{d}^{\pi_{(d)}}\\}$. Then the permutation
distribution of $T$ based on all possible local permutations is given as

$\displaystyle
F_{\bm{\pi},T}(x):=\mathbb{E}_{\pi_{(1)},\ldots,\pi_{(d)}}\bigl{[}\mathds{1}\bigl{(}T^{\bm{\pi}}\leq
x\bigr{)}\,\big{|}\,\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}\bigr{]}.$

In other words, $F_{\bm{\pi},T}$ is the conditional distribution of $T_{j}$ in
Algorithm 2 given $\\{(X_{i},Y_{i},Z_{i})\\}_{i=1}^{n}$. We denote the
$1-\alpha$ quantile of $F_{\bm{\pi},T}$ by $q_{1-\alpha}^{\dagger}$. Unlike
$q_{1-\alpha}$, the quantile of $F_{\bm{\pi},T}$ does not involve randomness
from Monte Carlo simulation and thus is relatively easier to handle
theoretically. In fact, when $B$ is sufficiently large, $q_{1-\alpha}$ and
$q_{1-\alpha}^{\dagger}$ are close to each other by the law of large numbers.
For our purpose, it is enough to show that $q_{1-\alpha}^{\dagger}$ is larger
than $q_{1-\alpha}$ with high probability. To prove this statement, note that
the Dvoretzky–Kiefer–Wolfowitz inequality guarantees that

$\displaystyle\sup_{x\in\mathbb{R}}\bigg{|}F_{\bm{\pi},T}(x)-\frac{1}{B}\sum_{i=1}^{B}\mathds{1}\big{(}T_{i}\leq
x\big{)}\bigg{|}\leq\sqrt{\frac{1}{2B}\log\left(\frac{4}{\beta}\right)}$

holds with probability at least $1-\beta/4$. Under this good event, it holds
that

$\displaystyle q_{1-\alpha}\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\inf\bigg{\\{}x\in\mathbb{R}:1-\alpha\leq\frac{1}{B+1}\sum_{i=0}^{B}\mathds{1}\big{(}T_{i}\leq
x\big{)}\bigg{\\}}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\inf\bigg{\\{}x\in\mathbb{R}:1-\alpha\leq\frac{1}{B+1}\sum_{i=1}^{B}\mathds{1}\big{(}T_{i}\leq
x\big{)}\bigg{\\}}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\inf\bigg{\\{}x\in\mathbb{R}:(1-\alpha)\frac{B+1}{B}\leq\frac{1}{B}\sum_{i=1}^{B}\mathds{1}\big{(}T_{i}\leq
x\big{)}\bigg{\\}}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\inf\bigg{\\{}x\in\mathbb{R}:(1-\alpha)\frac{B+1}{B}+\sqrt{\frac{1}{2B}\log\left(\frac{4}{\beta}\right)}\leq
F_{\bm{\pi},T}(x)\bigg{\\}}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle q_{1-\alpha/2}^{\dagger},$

where we utilize the assumption on
$B\geq\max\\{4(1-\alpha)\alpha^{-1},\,8\alpha^{-2}\log(4/\beta)\\}$ to ensure
that

$\displaystyle(1-\alpha)\frac{B+1}{B}+\sqrt{\frac{1}{2B}\log\left(\frac{4}{\beta}\right)}\leq
1-\frac{\alpha}{2}.$

Here constant $2$ in $\alpha/2$ is chosen for convenience and can be an
arbitrary positive number greater than one. This result combined with the
union bound yields that

$\displaystyle\mathbb{P}\bigl{(}q_{1-\alpha}\geq\zeta\sqrt{\min(n,d)}\bigr{)}\leavevmode\nobreak\
\leq\leavevmode\nobreak\
\mathbb{P}\bigl{(}q_{1-\alpha/2}^{\dagger}\geq\zeta\sqrt{\min(n,d)}\bigr{)}+\frac{\beta}{4}.$

Hence condition (8) is ensured if

$\displaystyle\mathbb{P}\bigl{(}q_{1-\alpha/2}^{\dagger}\geq\zeta\sqrt{\min(n,d)}\bigr{)}\leq\frac{\beta}{4}.$
(9)

Indeed, the above inequality is essentially proved in the proof of Theorem 5
in Kim et al., 2022b . Strictly speaking, Kim et al., 2022b  prove this under
Poissonization, but this is not critical in the proof of (9) and the same
lines of the proof go through under multinomial sampling. This completes the
proof of the first part of Theorem 2.

####  C.2.2 Claim for $T_{W}$

The proof of the second part is similar to the first part. Again, we consider
the permutation distribution of $T_{W}$ where the permutation only applies to
$D_{X,Y,1},\ldots,D_{X,Y,d}$. More specifically, the permutation distribution
function of $T_{W}$ at $x\in\mathbb{R}$ is given as the conditional
expectation of $\mathds{1}(T_{j,W}\leq x)$ given everything except
permutations. We denote the $1-\alpha$ quantile of the resulting permutation
distribution by $q_{1-\alpha,W}^{\dagger}$ and the $1-\alpha$ quantile of
$T_{W},T_{1,W},\ldots,T_{B,W}$ by $q_{1-\alpha,W}$. By following the same
logic in the first part of the proof, when
$B\geq\max\\{4(1-\alpha)\alpha^{-1},\,8\alpha^{-2}\log(4/\beta)\\}$, it can be
ensured that $q_{1-\alpha,W}\leq q_{1-\alpha,W}^{\dagger}$ with probability at
least $1-\beta$. Therefore, the type II error of the permutation test in
Algorithm 3 is bounded by

$\displaystyle\mathbb{P}(T_{W}\leq q_{1-\alpha,W})\leq\mathbb{P}(T_{W}\leq
q_{1-\alpha,W}^{\dagger})+\frac{\beta}{4}.$

Under Poissonization, Theorem 5 of Kim et al., 2022b  shows that

$\displaystyle\mathbb{P}(T_{W}\leq
q_{1-\alpha,W}^{\dagger})\leq\frac{3}{4}\beta.$

In fact, the same lines of the proof of Theorem 5 in Kim et al., 2022b
combined with the results in Appendix C.1.2 prove that the same conclusion
holds under multinomial sampling. Therefore we omit the details and finish the
proof.

###  C.3 Proof of Proposition 1

We start by making an important observation that
$U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ is equivalent to the U-statistic of
Hilbert–Schmidt Independence Criterion (Gretton et al.,, 2005) based on
kernels

$\displaystyle
k(x,x^{\prime})=\sum_{q=1}^{\ell_{1}}\frac{\mathds{1}(x=q)\mathds{1}(x^{\prime}=q)}{\eta_{q}}\quad\text{and}\quad
l(y,y^{\prime})=\sum_{r=1}^{\ell_{2}}\frac{\mathds{1}(y=r)\mathds{1}(y^{\prime}=r)}{\upsilon_{r}},$

for $x,x^{\prime}\in[\ell_{1}]$ and $y,y^{\prime}\in[\ell_{2}]$. Leveraging
this observation along with the computational trick introduced in Song et al.,
(2012), $U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ can be computed via convenient
matrix operations. Specifically, we let $\bm{\mathrm{K}},\bm{\mathrm{L}}$ be
$\sigma\times\sigma$ kernel matrices whose entries are
$\bm{\mathrm{K}}_{i,j}=k(X_{i},X_{j})$ and
$\bm{\mathrm{L}}_{i,j}=l(Y_{i},Y_{j})$, respectively. Define
$\widetilde{\bm{\mathrm{K}}}$ by letting
$\widetilde{\bm{\mathrm{K}}}_{i,j}=\bm{\mathrm{K}}_{i,j}$ if $i\neq j$ and
$\widetilde{\bm{\mathrm{K}}}_{i,j}=0$ if $i=j$. That is,
$\widetilde{\bm{\mathrm{K}}}$ is the same as $\bm{\mathrm{K}}$ except having
zero diagonal components. We similarly define $\widetilde{\bm{\mathrm{L}}}$ by
relating it to $\bm{\mathrm{L}}$. Denote the vector of size $\sigma$ with each
entry equaling one by $\bm{\mathrm{1}}$. Then equation (5) of Song et al.,
(2012) gives a computationally efficient form of
$U^{\bm{\eta},\bm{\upsilon}}_{W}(D)$ as

$\displaystyle U^{\bm{\eta},\bm{\upsilon}}_{W}(D)\leavevmode\nobreak\
=\leavevmode\nobreak\
\frac{1}{\sigma(\sigma-3)}\Biggl{[}\text{tr}(\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}})+\frac{\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\bm{\mathrm{1}}\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}}{(\sigma-1)(\sigma-2)}-\frac{2}{\sigma-2}\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}\Biggr{]}.$

It is therefore enough to show that
$A_{1}=\text{tr}(\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}})$,
$A_{2}=\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\bm{\mathrm{1}}\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}$
and
$A_{3}=\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}$.
We prove these in order.

  * •

For the claim
$A_{1}=\text{tr}(\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}})$, we
observe that

$\displaystyle\text{tr}(\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}})\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\bm{\mathrm{K}}_{i,j}\bm{\mathrm{L}}_{i,j}\mathds{1}(i\neq
j)$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\Biggl{[}\sum_{q=1}^{\ell_{1}}\frac{\mathds{1}(X_{i}=q)\mathds{1}(X_{j}=q)}{\eta_{q}}\Biggr{]}\Biggl{[}\sum_{r=1}^{\ell_{2}}\frac{\mathds{1}(Y_{i}=r)\mathds{1}(Y_{j}=r)}{\upsilon_{r}}\Biggr{]}\mathds{1}(i\neq
j)$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\mathds{1}(X_{i}=q)\mathds{1}(X_{j}=q)\mathds{1}(Y_{i}=r)\mathds{1}(Y_{j}=r)}{\eta_{q}\upsilon_{r}}\mathds{1}(i\neq
j)$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\Biggl{[}\frac{o_{qr}^{2}-o_{qr}}{\eta_{q}\upsilon_{r}}\Biggr{]}=A_{1}.$

  * •

To prove the second claim
$A_{2}=\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\bm{\mathrm{1}}\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}$,
notice that

$\displaystyle\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\bm{\mathrm{1}}\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{q=1}^{\ell_{1}}\frac{\mathds{1}(X_{i}=q)\mathds{1}(X_{j}=q)}{\eta_{q}}\mathds{1}(i\neq
j)$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\Biggl{[}\frac{o_{q+}^{2}-o_{q+}}{\eta_{q}}\Biggr{]}.$

Similarly, we have

$\displaystyle\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\sum_{r=1}^{\ell_{2}}\Biggl{[}\frac{o_{+r}^{2}-o_{+r}}{\upsilon_{r}}\Biggr{]}.$

Combining these two proves the identity
$A_{2}=\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\bm{\mathrm{1}}\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}$.

  * •

Finally, the claim
$A_{3}=\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}$
can be proved by following a series of identities:

$\displaystyle\bm{\mathrm{1}}^{\top}\widetilde{\bm{\mathrm{K}}}\widetilde{\bm{\mathrm{L}}}\bm{\mathrm{1}}\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\sum_{i=1}^{n}\Biggl{[}\sum_{j=1}^{n}\bm{\mathrm{K}}_{i,j}\mathds{1}(i\neq
j)\Biggr{]}\Biggl{[}\sum_{j^{\prime}=1}^{n}\bm{\mathrm{L}}_{i,j^{\prime}}\mathds{1}(i\neq
j^{\prime})\Biggr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{i=1}^{n}\Biggl{[}\sum_{j=1}^{n}\sum_{q=1}^{\ell_{1}}\frac{\mathds{1}(X_{i}=q)\mathds{1}(X_{j}=q)}{\eta_{q}}\mathds{1}(i\neq
j)\Biggr{]}\Biggl{[}\sum_{j^{\prime}=1}^{n}\sum_{r=1}^{\ell_{2}}\frac{\mathds{1}(Y_{i}=r)\mathds{1}(Y_{j^{\prime}}=r)}{\upsilon_{r}}\mathds{1}(i\neq
j^{\prime})\Biggr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\sum_{i=1}^{n}\Biggl{[}\sum_{j=1}^{n}\frac{\mathds{1}(X_{i}=q)\mathds{1}(X_{j}=q)}{\eta_{q}}\mathds{1}(i\neq
j)\Biggr{]}\Biggl{[}\sum_{j^{\prime}=1}^{n}\frac{\mathds{1}(Y_{i}=r)\mathds{1}(Y_{j^{\prime}}=r)}{\upsilon_{r}}\mathds{1}(i\neq
j^{\prime})\Biggr{]}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}o_{q+}o_{+r}}{\eta_{q}\upsilon_{r}}-\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}o_{q+}}{\eta_{q}\upsilon_{r}}-\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}o_{+r}}{\eta_{q}\upsilon_{r}}+\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}}{\eta_{q}\upsilon_{r}}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{q=1}^{\ell_{1}}\sum_{r=1}^{\ell_{2}}\frac{o_{qr}(o_{q+}o_{+r}-o_{q+}-o_{+r}+1)}{\eta_{q}\upsilon_{r}}$
$\displaystyle=\leavevmode\nobreak\ $ $\displaystyle A_{3}.$

This completes the proof of Proposition 1.

###  C.4 Proof of Proposition 2

Our strategy is to find a sequence of distributions under the given conditions
for which the asymptotic power of $\chi^{2}$\- and $G$-tests are zero. This
result implies the desired claim since we are concerned with the worse case
power. Let us start by forming a sequence of discrete distributions on
$[2]\times[2]\times[d]$.

Worst case scenario. Consider the setting where
$p_{X,Y|Z}(1,1\,|\,z)=p_{X,Y|Z}(2,2\,|\,z)=1/2$ and
$p_{X,Y|Z}(1,2\,|\,z)=p_{X,Y|Z}(2,1\,|\,z)=0$ for all $z\in[d]$. For each $n$,
set $p_{Z}(1)=(1-1/n)^{1/n}$ so that the probability of having all of
$Z_{1},\ldots,Z_{n}$ equal to one is $1-1/n$. We then set
$p_{Z}(2)=\ldots=p_{Z}(d)=\bigl{(}1-p_{Z}(1)\bigr{)}/(d-1)$. Notice that by
construction, $p_{Z}$ takes non-zero values over $[d]$ and thus the dimension
of $Z$ is $d$. We then let $p=p_{X,Y|Z}p_{Z}$ and claim that

$\displaystyle\inf_{q\in\mathcal{P}_{0}}\|p-q\|_{1}\geq\frac{1}{4}.$ (10)

Note that for any $q=q_{X|Z}q_{Y|Z}q_{Z}\in\mathcal{P}_{0}$,

$\displaystyle\|p-q\|_{1}\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\|p_{X,Y|Z}p_{Z}-q_{X|Z}q_{Y|Z}q_{Z}\|_{1}$
$\displaystyle\overset{\mathrm{(i)}}{\geq}\leavevmode\nobreak\ $
$\displaystyle\|p_{X,Y|Z}p_{Z}-q_{X|Z}q_{Y|Z}p_{Z}\|_{1}-\|q_{X|Z}q_{Y|Z}p_{Z}-q_{X|Z}q_{Y|Z}q_{Z}\|_{1}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\sum_{z\in[d]}p_{Z}(z)\cdot\|p_{X,Y|Z=z}-q_{X|Z=z}q_{Y|Z=z}\|_{1}-\|p_{Z}-q_{Z}\|_{1}$
$\displaystyle\overset{\mathrm{(ii)}}{\geq}\leavevmode\nobreak\ $
$\displaystyle\sum_{z\in[d]}p_{Z}(z)\cdot\|p_{X,Y|Z=z}-q_{X|Z=z}q_{Y|Z=z}\|_{1}-\|p-q\|_{1},$

where step (i) uses the triangle inequality and step (ii) uses the following
inequality (which can be proved using the triangle inequality):

$\displaystyle\|p_{Z}-q_{Z}\|_{1}=\sum_{z\in[d]}|p_{Z}(z)-q_{Z}(z)|=\sum_{z\in[d]}\Bigg{|}\sum_{x\in[\ell_{1}],y\in[\ell_{2}]}p(x,y,z)-q(x,y,z)\Biggr{|}\leq\|p-q\|_{1}.$

Therefore we can conclude that

$\displaystyle\inf_{q\in\mathcal{P}_{0}}\|p-q\|_{1}\geq\frac{1}{2}\inf_{q\in\mathcal{P}_{0}}\Biggl{\\{}\sum_{z\in[d]}p_{Z}(z)\cdot\|p_{X,Y|Z=z}-q_{X|Z=z}q_{Y|Z=z}\|_{1}\Biggr{\\}}.$

In order to further lower bound the above quantity, we first prove a
preliminary result. Let $\mathcal{Q}_{0}$ be the collection of joint
distributions of $(X,Y)$ on $[2]\times[2]$ such that for any
$q_{X,Y}\in\mathcal{Q}_{0}$, $q_{X,Y}=q_{X}q_{Y}$, i.e.,
$X\,\rotatebox[origin={c}]{90.0}{$\models$}\,Y$. Then for a joint distribution
$p_{X,Y}$ such that $p_{X,Y}(1,1)=p_{X,Y}(2,2)=1/2$ and
$p_{X,Y}(1,2)=p_{X,Y}(2,1)=0$, we have

$\displaystyle\inf_{q\in\mathcal{Q}_{0}}\|p_{X,Y}-q\|_{1}\geq\inf_{q\in\mathcal{Q}_{0}}\|p_{X,Y}-q\|_{2}=\frac{1}{2},$

where the last equality can be shown by solving the following optimization
problem:

$\displaystyle\inf_{a,b\in[0,1]}\sqrt{\Bigl{(}\frac{1}{2}-ab\Bigr{)}^{2}+a^{2}(1-b)^{2}+(1-a)^{2}b^{2}+\Bigl{(}\frac{1}{2}-(1-a)(1-b)\Bigr{)}^{2}}=\frac{1}{2},$

where the infimum is achieved at $a=b=1/2$. Consequently, given that
$(X,Y)\,\rotatebox[origin={c}]{90.0}{$\models$}\,Z$ in our construction,

$\displaystyle\inf_{q\in\mathcal{P}_{0}}\|p-q\|_{1}\leavevmode\nobreak\
\geq\leavevmode\nobreak\ $
$\displaystyle\frac{1}{2}\inf_{q\in\mathcal{P}_{0}}\Biggl{\\{}\sum_{z\in[d]}p_{Z}(z)\cdot\|p_{X,Y|Z=z}-q_{X|Z=z}q_{Y|Z=z}\|_{1}\Biggr{\\}}$
$\displaystyle\geq\leavevmode\nobreak\ $ $\displaystyle\frac{1}{4},$

which proves claim (10).

Main proof. Now let $\mathcal{A}_{n}$ be an event that all of
$Z_{1},\ldots,Z_{n}$ are equal to one, and $\mathcal{A}^{c}$ be its
complement. Starting with $\chi^{2}$-test, its worst case power is bounded by

$\displaystyle\inf_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(\chi^{2}>q_{1-\alpha,d})\leavevmode\nobreak\
\leq\leavevmode\nobreak$
$\displaystyle\mathbb{P}_{p^{\ast}}(\chi^{2}>q_{1-\alpha,d})$ (11)
$\displaystyle\leq\leavevmode\nobreak$
$\displaystyle\mathbb{P}_{p^{\ast}}(\chi^{2}>q_{1-\alpha,d}\,|\,\mathcal{A}_{n})\mathbb{P}_{p^{\ast}}(\mathcal{A}_{n})+\mathbb{P}_{p^{\ast}}(\mathcal{A}_{n}^{c})$
$\displaystyle\leq\leavevmode\nobreak$
$\displaystyle\mathbb{P}_{p^{\ast}}(\chi^{2}>q_{1-\alpha,d}\,|\,\mathcal{A}_{n})+\frac{1}{n}.$

Conditional on $\mathcal{A}_{n}$, $\chi^{2}/n$ has the same distribution as

$\displaystyle
W^{\star}_{n}:=\frac{(W-W^{2}/n)^{2}}{W^{2}}+\frac{2\bigl{(}W^{2}(n-W)^{2}/n^{2}\bigr{)}}{W(n-W)}+\frac{\bigl{(}n-W-(n-W)^{2}/n\bigr{)}^{2}}{(n-W)^{2}},$

where $W\sim\mathrm{Binomial}(n,1/2)$. Using this notation, we have the
identity

$\displaystyle\mathbb{P}_{p^{\ast}}(\chi^{2}>q_{1-\alpha,d}\,|\,\mathcal{A}_{n})\leavevmode\nobreak\
=\leavevmode\nobreak\ \mathbb{P}\biggl{(}\frac{n}{d}\cdot
W^{\star}_{n}-1>\frac{1}{\sqrt{d}}\cdot\frac{q_{1-\alpha,d}-d}{\sqrt{d}}\biggr{)}.$

Under our choice of $d=n\times r_{n}$ where $r_{n}\rightarrow\infty$, we apply
the law of large number and continuous mapping theorem and see that

$\displaystyle\frac{n}{d}\cdot W^{\star}_{n}-1\overset{p}{\longrightarrow}-1,$

where $X_{n}\overset{p}{\longrightarrow}X$ means convergence of $X_{n}$ to $X$
in probability. In addition, by the central limit theorem along the fact that
convergence in distribution implies convergence of the quantile function at
all continuity points, we have

$\displaystyle\frac{q_{1-\alpha}-d}{\sqrt{2d}}\rightarrow
z_{1-\alpha}\quad\text{as $d\rightarrow\infty$,}$

where $z_{1-\alpha}$ is the $1-\alpha$ quantile of $N(0,1)$. Thus

$\displaystyle\frac{1}{\sqrt{d}}\cdot\frac{q_{1-\alpha,d}-d}{\sqrt{d}}\rightarrow
0\quad\text{as $n\rightarrow\infty$.}$

Therefore by Slutsky’s theorem, we have

$\displaystyle\mathbb{P}_{p^{\ast}}(\chi^{2}>q_{1-\alpha,d}\,|\,\mathcal{A}_{n})=\mathbb{P}\biggl{(}\frac{n}{d}\cdot
W^{\star}_{n}-1>\frac{1}{\sqrt{d}}\cdot\frac{q_{1-\alpha,d}-d}{\sqrt{d}}\biggr{)}\rightarrow
0.$

This in conjunction with the upper bound in (11) implies that

$\displaystyle\lim_{n\rightarrow\infty}\inf_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(\chi^{2}>q_{1-\alpha,d})=0.$

The proof of $G$-test is essentially the same. The only difference is that
conditional on $\mathcal{A}_{n}$, $G/n$ has the same distribution as

$\displaystyle
V_{n}^{\star}:=2\biggl{[}\frac{W}{n}\log\biggl{(}\frac{n}{W}\biggr{)}+\biggl{(}1-\frac{W}{n}\biggr{)}\\!\log\biggl{(}\frac{1}{1-\frac{W}{n}}\biggr{)}\biggr{]},$

from which it holds that

$\displaystyle\mathbb{P}_{p^{\ast}}(G>q_{1-\alpha,d}\,|\,\mathcal{A}_{n})\leavevmode\nobreak\
=\leavevmode\nobreak\ \mathbb{P}\biggl{(}\frac{n}{d}\cdot
V^{\star}_{n}-1>\frac{1}{\sqrt{d}}\cdot\frac{q_{1-\alpha,d}-d}{\sqrt{d}}\biggr{)}.$

Again, it can be seen that $\frac{n}{d}\cdot
V^{\star}_{n}-1\overset{p}{\longrightarrow}-1$ and
$\frac{1}{\sqrt{d}}\cdot\frac{q_{1-\alpha,d}-d}{\sqrt{d}}\rightarrow 0$.
Following the same steps as in the case of $\chi^{2}$-test, we conclude

$\displaystyle\lim_{n\rightarrow\infty}\inf_{p\in\mathcal{P}_{1}(\varepsilon)}\mathbb{P}_{p}(G>q_{1-\alpha,d})=0.$

This completes the proof of Proposition 2.

###  C.5 Proof of Lemma 1

We start by computing the expectation of $X\mathds{1}(X\geq 4)$ explicitly as

$\displaystyle\mathbb{E}[X\mathds{1}(X\geq 4)]\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}[X]-\mathbb{E}[X\mathds{1}(X\leq 3)]$
$\displaystyle=\leavevmode\nobreak\ $ $\displaystyle np-
np(1-p)^{n-1}-n(n-1)p^{2}(1-p)^{n-2}-\frac{n(n-1)(n-2)}{2}p^{3}(1-p)^{n-3}.$

Given this formula, we consider two cases: (i) $np>1$ and (ii) $np\leq 1$,
separately. Throughout the proof, we assume that $p>0$ since the claim is
trivial when $p=0$.

Case (i) $np>1$. First assume that $np>1$, equivalently $p>1/n$. At a high-
level, when $np$ is large, there is only a small possibility that
$\mathds{1}(X<4)=1$. In this case, $X\mathds{1}(X\geq 4)$ behaves like $X$ and
we expect that $\mathbb{E}[X\mathds{1}(X\geq 4)]\approx\mathbb{E}[X]=np$. We
now make this statement more precise. To this end, let us define

$\displaystyle r_{n}(p)\leavevmode\nobreak\ :=\leavevmode\nobreak\ $
$\displaystyle\frac{\mathbb{E}[X\mathds{1}(X\geq 4)]}{np}$
$\displaystyle=\leavevmode\nobreak\ $ $\displaystyle
1-(1-p)^{n-1}-(n-1)p(1-p)^{n-2}-\frac{(n-1)(n-2)}{2}p^{2}(1-p)^{n-3}.$

We will show that $r_{n}(p)$ is an increasing function of $p$ for $p>1/n$ and
thus the minimum is achieved at $r_{n}(1/n)$. To show this, we calculate its
derivative with respect to $p$, which can be seen to be non-negative for all
$n\geq 4$ as

$\displaystyle
r_{n}^{\prime}(p)=\frac{(n-1)(n-2)(n-3)(1-p)^{n}p^{2}}{2(1-p)^{4}}\geq
0\quad\text{for all $n\geq 4$.}$

Therefore, the minimum of $r_{n}(p)$ is achieved at $p=1/n$ as a function of
$p$, which yields

$\displaystyle r_{n}(1/n)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle
1-\left(1-\frac{1}{n}\right)^{n-1}-\frac{n-1}{n}\left(1-\frac{1}{n}\right)^{n-2}-\frac{(n-1)(n-2)}{2n^{2}}\left(1-\frac{1}{n}\right)^{n-3}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\biggl{(}1-\frac{1}{n}\biggr{)}^{n}\biggl{\\{}\frac{n(5n-6)}{2(n-1)^{2}}\biggr{\\}}$
$\displaystyle\overset{(\ast)}{\geq}\leavevmode\nobreak\ $ $\displaystyle
e^{-1}\biggl{(}\frac{n-1}{n}\biggr{)}\biggl{\\{}\frac{n(5n-6)}{2(n-1)^{2}}\biggr{\\}}=e^{-1}\biggl{\\{}\frac{5}{2}-\frac{1}{2(n-1)}\biggr{\\}},$

where step ($\ast$) follows by the inequality

$\displaystyle\left(1+\frac{x}{n}\right)^{n}\geq
e^{x}\left(1-\frac{x^{2}}{n}\right)\quad\text{for all $n\geq 1$ and $|x|\leq
n$}.$

From this, we conclude that $r_{n}(1/n)\geq e^{-1}\times\frac{7}{3}\approx
0.858$ for all $n\geq 4$ and thus

$\displaystyle\mathbb{E}[X\mathds{1}(X\geq 4)]\geq\gamma np\quad\text{if
$np>1$.}$

Case (ii) $np\leq 1$. Next assume $np\leq 1$, i.e., $p\leq 1/n$. Then
$\min\\{np,(np)^{4}\\}=(np)^{4}$ and the ratio function $r_{n}(p)$ is computed
as

$\displaystyle r_{n}(p)\leavevmode\nobreak\ :=\leavevmode\nobreak\ $
$\displaystyle\frac{\mathbb{E}[X\mathds{1}(X\geq 4)]}{n^{4}p^{4}}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\frac{1}{(np)^{3}}-\frac{(1-p)^{n-1}}{(np)^{3}}-\frac{(n-1)}{n^{3}p^{2}}(1-p)^{n-2}-\frac{(n-1)(n-2)}{2n^{3}p}(1-p)^{n-3}.$

Given this formula, we aim to prove that $r_{n}(p)$ is a decreasing function
of $p$ for $p\leq 1/n$ and hence the minimum is achieved at
$r_{n}(1/n)\geq\gamma$ as shown earlier. To prove this, by assuming $p\neq 0$,
it can be seen that the derivative of $r_{n}(p)$ as a function of $p$ is given
by

$\displaystyle r_{n}^{\prime}(p)\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{(1-p)^{n}\\{(n-4)p[(n-3)p((n-2)p+3)+6]+6\\}-6(1-p)^{4}}{2n^{3}(1-p)^{4}p^{4}}.$

Our goal is to show that $r_{n}^{\prime}(p)\leq 0$ for any $p\in(0,1/n]$,
which is equivalent to the conditions

$\displaystyle r_{n}^{\prime}(p)\leq 0$
$\displaystyle\overset{\mathrm{iff}}{\Longleftrightarrow}\leavevmode\nobreak\
$ $\displaystyle
6(1-p)^{n}-6(1-p)^{4}+(1-p)^{n}(n-4)p\\{(n^{2}-5n+6)p^{2}+3(n-3)p+6\\}\leq 0$
$\displaystyle\overset{\mathrm{iff}}{\Longleftrightarrow}\leavevmode\nobreak\
$ $\displaystyle 6-6(1-p)^{4-n}+(n-4)\\{(n-2)(n-3)p^{3}+3(n-3)p^{2}+6p\\}\leq
0.$ (12)

To this end, we make a key (and nontrivial) observation that the function
$h_{n}(p):=6(1-p)^{4-n}$ has a Taylor expansion around $0$ as

$\displaystyle h_{n}(p)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle 6+6(n-4)p+3(n-4)(n-3)p^{2}+(n-4)(n-3)(n-2)p^{3}$
$\displaystyle+\underbrace{\frac{(n-4)(n-3)(n-2)(n-1)}{(1-\xi)^{n}}p^{4}}_{\mathrm{remainder}},$

for some $\xi\in(0,1/n)$. Importantly, the remainder term is strictly
positive, which concludes that the inequality (12) holds as well as
$r_{n}^{\prime}(p)\leq 0$ for all $p\in(0,1/n]$. This means that $r_{n}(p)$ is
a decreasing function and $r_{n}(1/n)\leq r_{n}(p)$ for $p\in(0,1/n]$. From
the previous calculation in case (i), we know that $r_{n}(1/n)\geq
e^{-1}\times\frac{7}{3}$ for all $n\geq 4$ (when $p=1/n$, we have
$np=(np)^{4}$) and thereby prove the claim.

###  C.6 Proof of Lemma 2

We divide the proof into two parts depending on whether (i) $np>1$ or (ii)
$np\leq 1$. Throughout, we assume that $p>0$ since the claim is trivial when
$p=0$.

Case (i) $np>1$. The statement is relatively easier to prove when $np>1$.
Again, the intuition is that when $np$ is large, $X\mathds{1}(X\geq 4)$
behaves like $X$ and thus $\mathrm{Var}[X\mathds{1}(X\geq
4)]\approx\mathrm{Var}[X]=np(1-p)$ and $\mathbb{E}[X\mathds{1}(X\geq
4)]\approx\mathbb{E}[X]=np$. Thus the expectation is larger than the variance.
To make this statement more precise, note that the variance can be bounded by

$\displaystyle\mathrm{Var}[X\mathds{1}(X\geq 4)]\leavevmode\nobreak\
=\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}[X^{2}]-\mathbb{E}[X^{2}\mathds{1}(X\leq
3)]-\\{\mathbb{E}[X]-\mathbb{E}[X\mathds{1}(X\geq 3)]\\}^{2}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\mathbb{E}[X^{2}]-\\{\mathbb{E}[X]\\}^{2}+2\mathbb{E}[X]\mathbb{E}[X\mathds{1}(X\geq
3)]$ $\displaystyle\leq\leavevmode\nobreak\ $ $\displaystyle
np(1-p)+6np\mathbb{P}(X\geq 3)$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle 7np.$

From Lemma 1, we know that $\mathbb{E}[X\mathds{1}(X\geq 4)]\geq Cnp$ for some
$C>0$ when $np>1$. Therefore, by taking $\gamma$ larger than, for example,
$7/0.85$, the desired result follows.

Case (ii) $np\leq 1$.  Now we assume that $np\leq 1$, i.e., $p\leq 1/n$. In
this case, we cannot ignore the effect of the indicator function
$\mathds{1}(X\geq 4)$ as in the case of $np>1$ and thus it requires a more
delicate analysis. To simplify the problem a bit, we first upper bound the
ratio between the variance and the expectation as

$\displaystyle\frac{\mathrm{Var}[X\mathds{1}(X\geq
4)]}{\mathbb{E}[X\mathds{1}(X\geq
4)]}\lesssim\frac{\mathbb{E}[X^{2}\mathds{1}(X\geq 4)]}{(np)^{4}},$

where we use the fact that $\mathrm{Var}[X\mathds{1}(X\geq
4)]\leq\mathbb{E}[X^{2}\mathds{1}(X\geq 4)]$ and $\mathbb{E}[X\mathds{1}(X\geq
4)]\gtrsim(np)^{4}$ when $np\leq 1$ from Lemma 1. Our goal is to show that

$\displaystyle r_{n}(p):=\frac{\mathbb{E}[X^{2}\mathds{1}(X\geq
4)]}{(np)^{4}}$

is upper bounded by some positive constant for all $n\geq 4$ and $p\leq 1/n$.
In particular, we will show that $r_{n}(p)$ is a decreasing function in $p$
for all $n\geq 4$ and

$\displaystyle\lim_{p\rightarrow
0}r_{n}(p)=\frac{2(n^{3}-6n^{2}+11n-6)}{3n^{3}}<0.67\quad\text{for all $n\geq
4$.}$

To verify these claims, note that the second moment of $X\mathds{1}(X\geq 4)$:

$\displaystyle\mathbb{E}[X^{2}\mathds{1}(X\geq 4)]\leavevmode\nobreak\
=\leavevmode\nobreak\ \mathbb{E}[X^{2}]-\mathbb{E}[X^{2}\mathds{1}(X\leq 3)]$
$\displaystyle=\leavevmode\nobreak\ $ $\displaystyle
np(1-p)+(np)^{2}-np(1-p)^{n-1}-2n(n-1)p^{2}(1-p)^{n-2}-\frac{3}{2}n(n-1)(n-2)p^{3}(1-p)^{n-3},$

and therefore $r_{n}(p)$ can be written as

$\displaystyle r_{n}(p)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\frac{1-p}{(np)^{3}}+\frac{1}{(np)^{2}}-\frac{(1-p)^{n-1}}{(np)^{3}}-\frac{2(n-1)(1-p)^{n-2}}{n^{3}p^{2}}-\frac{3(n-1)(n-2)(1-p)^{n-3}}{2n^{3}p}.$

To show that $r_{n}(p)$ is a decreasing function in $p$, consider its
derivative with respect to $p$:

$\displaystyle r_{n}^{\prime}(p)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle\dfrac{3\left(n-3\right)\left(n-2\right)\left(n-1\right)\left(1-p\right)^{n-4}}{2n^{3}p}+\dfrac{7\left(n-2\right)\left(n-1\right)\left(1-p\right)^{n-3}}{2n^{3}p^{2}}+\dfrac{5\left(n-1\right)\left(1-p\right)^{n-2}}{n^{3}p^{3}}$
$\displaystyle+\dfrac{3\left(1-p\right)^{n-1}}{n^{3}p^{4}}-\dfrac{2}{n^{2}p^{3}}-\frac{1}{n^{3}p^{3}}-\dfrac{3(1-p)}{n^{3}p^{4}}.$

We would like to show that $r_{n}^{\prime}(p)$ is less than or equal to $0$
for all $0<p\leq 1/n$ and $n\geq 4$, which is equivalent to verifying the
following slightly simplified conditions:

$\displaystyle r_{n}^{\prime}(p)\leq 0$
$\displaystyle\overset{\mathrm{iff}}{\Longleftrightarrow}\quad$
$\displaystyle\frac{3}{2}\left(n-3\right)\left(n-2\right)\left(n-1\right)np^{3}\left(1-p\right)^{n-4}+\frac{7}{2}\left(n-2\right)\left(n-1\right)np^{2}\left(1-p\right)^{n-3}$
$\displaystyle+5\left(n-1\right)np\left(1-p\right)^{n-2}+3n\left(1-p\right)^{n-1}-2n^{2}p-np-3n(1-p)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0$
$\displaystyle\overset{\mathrm{iff}}{\Longleftrightarrow}\quad$ $\displaystyle
3\left(n-3\right)\left(n-2\right)\left(n-1\right)np^{3}\left(1-p\right)^{n-4}+7\left(n-2\right)\left(n-1\right)np^{2}\left(1-p\right)^{n-3}$
$\displaystyle+10\left(n-1\right)np\left(1-p\right)^{n-2}+6n\left(1-p\right)^{n-1}-4n^{2}p-2np-6n(1-p)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0$
$\displaystyle\overset{\mathrm{iff}}{\Longleftrightarrow}\quad$ $\displaystyle
3\left(n-3\right)\left(n-2\right)\left(n-1\right)np^{4}\left(1-p\right)^{n-4}+7\left(n-2\right)\left(n-1\right)np^{3}\left(1-p\right)^{n-3}$
$\displaystyle+10\left(n-1\right)np^{2}\left(1-p\right)^{n-2}+6np\left(1-p\right)^{n-1}-4n^{2}p^{2}-2np^{2}-6np(1-p)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0.$

Focusing on the last equivalent condition, we first define

$\displaystyle f_{n}(p)\leavevmode\nobreak\ :=\leavevmode\nobreak\ $
$\displaystyle
3\left(n-3\right)\left(n-2\right)\left(n-1\right)np^{4}\left(1-p\right)^{n-4}+7\left(n-2\right)\left(n-1\right)np^{3}\left(1-p\right)^{n-3}$
$\displaystyle+10\left(n-1\right)np^{2}\left(1-p\right)^{n-2}+6np\left(1-p\right)^{n-1}.$

Then the last condition simply becomes

$\displaystyle f_{n}(p)\leq 6np+4n(n-1)p^{2}.$ (13)

Somewhat surprisingly, it can be seen that $6np+4n(n-1)p^{2}$ is the second-
order Taylor expansion of $f_{n}(p)$ near $p=0$. In particular, we have

$\displaystyle f_{n}(p)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle
f_{n}(0)+f_{n}^{\prime}(0)p+\frac{f_{n}^{\prime\prime}(0)}{2}p^{2}+R_{n}$
$\displaystyle=\leavevmode\nobreak\ $ $\displaystyle 6np+4n(n-1)p^{2}+R_{n},$

where

$\displaystyle R_{n}=\frac{f_{n}^{(3)}(\zeta)}{6}p^{3}\quad\text{for some
$\zeta\in(0,1/n)$.}$

Therefore claim (13) follows once we show that $R_{n}\leq 0$ for all $n\geq 4$
and $p\leq 1/n$. Note that the third derivative of $f_{n}(p)$ is computed as

$\displaystyle f_{n}^{(3)}(p)$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\dfrac{n(n-1)(n-2)(n-3)(n-4)(1-p)^{n}p^{2}\\{n(3n-4)p^{2}+(35-29n)p+55\\}}{(p-1)^{7}}$

and since the denominator $(p-1)^{7}$ is negative, it holds that
$f_{n}^{(3)}(p)\leq 0$ once

$\displaystyle k_{n}(p):=n(3n-4)p^{2}+(35-29n)p+55$

is positive for $p\in(0,1/n]$. Note that $k_{n}(p)$ is a quadratic function of
$p$ and its global minimum is achieved at
$p^{\ast}_{n}=\frac{29n-35}{2n(3n-4)}$. Moreover,

$\displaystyle
p^{\ast}_{n}-\frac{1}{n}=\frac{23n-27}{2n(3n-4)}>0\quad\text{for $n\geq 4$,}$

which means that $k_{n}(p)$ is decreasing for all $p\leq\frac{1}{n}$. On the
other hand, $k_{n}(1/n)$ is computed as

$\displaystyle k_{n}(1/n)=29+\frac{23}{n}>0\quad\text{for all $n\geq 4$.}$

Therefore $k_{n}(p)$ is positive on $p\in(0,1/n]$ for all $n\geq 4$, which in
turn implies $R_{n}\leq 0$. In conclusion, the claim in (13) holds for all
$n\geq 4$ and therefore we complete the proof of Lemma 2.

###  C.7 Proof of Lemma 3

In contrast to a Poisson random variable that is characterized by a single
parameter $\lambda$, a binomial random variable depends on two parameters $n$
and $p$. As in the case of Lemma 2, this makes the proof significantly more
challenging than the proof for a Poisson random variable. Indeed, we need to
take a totally different approach than Canonne et al., (2018) and Kim, (2022)
to prove the given statement. The key ingredient of our proof is the
Efron–Stein inequality recalled in Lemma 8.

Let $X_{1},\ldots,X_{n}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Binomial}(1,p)$.
Then by letting $X=\sum_{i=1}^{n}X_{i}$, we can represent the quantity of
interest as a function of i.i.d. random variables as

$\displaystyle
f(X_{1},\ldots,X_{n})=X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq 4).$

For simplicity, we denote $f(X_{1},\ldots,X_{n})$ by $f$. Then the Efron–Stein
inequality (Lemma 8) yields that

$\displaystyle\mathrm{Var}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}\leq\frac{1}{2}\sum_{i=1}^{n}\mathbb{E}\big{[}\big{(}f-f^{i}\big{)}^{2}\big{]},$

where we define
$f^{i}=f(X_{1},\ldots,X_{i-1},X_{i}^{\prime},X_{i+1},\ldots,X_{n})$ and
$X_{i}^{\prime}\sim\mathrm{Binomial}(1,p)$ independent of
$X_{1},\ldots,X_{n}$. By symmetry, it holds that

$\displaystyle\mathrm{Var}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}\leq\frac{n}{2}\mathbb{E}\big{[}\big{(}f-f^{1}\big{)}^{2}\big{]},$

and we let $Y=\sum_{i=2}^{n}X_{i}$ for simplicity so that $X=X_{1}+Y$. By the
law of total expectation,

$\displaystyle
n\mathbb{E}\big{[}\big{(}f-f^{1}\big{)}^{2}\big{]}\leavevmode\nobreak\
=\leavevmode\nobreak\ $ $\displaystyle
n\mathbb{E}_{Y}\big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\big{]}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\underbrace{n\mathbb{E}_{Y}\big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\mathds{1}(Y<4)\big{]}}_{(\mathrm{I})}+\underbrace{n\mathbb{E}_{Y}\big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\mathds{1}(Y\geq
4)\big{]}}_{(\mathrm{II})}.$

We analyze two terms (I) and (II) separately.

Analysis of (I). We first claim that

$\displaystyle\mathbb{E}_{Y}\big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\mathds{1}(Y<4)\big{]}=\mathbb{E}_{Y}\big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\mathds{1}(Y=3)\big{]}.$

To see this, we unpack the expression of $f$ and $f^{1}$ as

$\displaystyle(f-f^{1})^{2}\mathds{1}(Y<4)$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\Big{\\{}(X_{1}+Y)\sqrt{\min\\{(X_{1}+Y),a\\}\min\\{(X_{1}+Y),b\\}}\mathds{1}(X_{1}+Y\geq
4)$
$\displaystyle-(X_{1}^{\prime}+Y)\sqrt{\min\\{(X_{1}^{\prime}+Y),a\\}\min\\{(X_{1}^{\prime}+Y),b\\}}\mathds{1}(X_{1}^{\prime}+Y\geq
4)\Big{\\}}^{2}\mathds{1}(Y<4).$

When $Y\in\\{0,1,2\\}$, the indicator functions $\mathds{1}(X_{1}+Y\geq 4)$
and $\mathds{1}(X_{1}+Y\geq 4)$ are equal to zero since $X_{1}$ and
$X_{1}^{\prime}$ are either zero or one. Therefore, we only need to focus on
the case of $Y=3$. Conditional on $Y=3$, we make use of the fact that $X_{1}$
and $X_{1}^{\prime}$ are independent and show that the expectation of
$(f-f^{1})^{2}\mathds{1}(Y=3)$ is

$\displaystyle
2p(1-p)(Y+1)\sqrt{\min\\{(Y+1),a\\}\min\\{(Y+1),b\\}}\mathds{1}(Y=3).$

Now by taking the expectation over $Y\sim\mathrm{Binomial}(n-1,p)$,

$\displaystyle\mathbb{E}_{Y}\Big{[}\mathbb{E}_{X_{1},X_{1}^{\prime}}\big{[}(f-f^{1})^{2}\,|\,Y\big{]}\mathds{1}(Y<4)\Big{]}\leavevmode\nobreak\
=\leavevmode\nobreak\ $ $\displaystyle
8p(1-p)\sqrt{\min\\{4,a\\}\min\\{4,b\\}}\binom{n-1}{3}p^{3}(1-p)^{n-4}.$

Therefore by using the condition that $n\geq 4$ and the inequality $(1+x)\leq
e^{x}$ for $x\in\mathbb{R}$,

$\displaystyle n(\mathrm{I})\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
(np)^{4}(1-p)^{n-3}\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
\lambda^{4}\exp(-\lambda/4).$

Moreover, from Lemma 4, we know that

$\displaystyle\mathbb{E}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}\gtrsim\min\bigl{\\{}\lambda\sqrt{\min(\lambda,a)\min(\lambda,b)},\lambda^{4}\bigr{\\}}.$

Hence, once we show

$\displaystyle\lambda^{4}\exp(-\lambda/4)\lesssim\min\bigl{\\{}\lambda\sqrt{\min(\lambda,a)\min(\lambda,b)},\lambda^{4}\bigr{\\}},$
(14)

then it implies that

$\displaystyle
n(\mathrm{I})\lesssim\mathbb{E}\bigl{[}X\sqrt{\min\\{X,a\\}\min\\{X,b\\}}\mathds{1}(X\geq
4)\bigr{]}.$ (15)

For $\lambda<1$, the dominating term in the upper bound on inequality (14) is
$\lambda^{4}$ and thus inequality (14) follows easily. For $\lambda\geq 1$, by
assuming $a\leq b$ without loss of generality, it is enough to show

$\displaystyle\lambda^{4}\exp(-\lambda/4)\lesssim\lambda^{2},$
$\displaystyle\lambda^{4}\exp(-\lambda/4)\lesssim\lambda\sqrt{\lambda}\sqrt{a},\quad\text{for
$\lambda\leq a$},$
$\displaystyle\lambda^{4}\exp(-\lambda/4)\lesssim\lambda\sqrt{ab},\quad\text{for
$\lambda\lesssim\min\\{a,b\\}$.}$

Starting with the first inequality, the Taylor expansion of
$e^{x}=\sum_{n=0}^{\infty}\frac{x^{n}}{n!}$ shows that
$\exp(\lambda/4)\gtrsim\lambda^{2}$ for all $\lambda\geq 1$. For the second
inequality, again, the Taylor expansion of $e^{x}$ shows that
$\exp(\lambda/4)\gtrsim\lambda^{2}\gtrsim\lambda^{5/2}/\sqrt{a}$ for
$\lambda\leq a$. The last inequality follows similarly. Hence inequality (15)
follows as desired.

Analysis of (II). Switching to the second term (II), note that when $Y\geq 4$,
it is always the case that $X_{1}+Y\geq 4$ and $X_{1}^{\prime}+Y\geq 4$.
Therefore we can write

$\displaystyle(f-f^{1})^{2}\mathds{1}(Y\geq 4)$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\bigl{\\{}(X_{1}+Y)\sqrt{\min\\{(X_{1}+Y),a\\}\min\\{(X_{1}+Y),b\\}}$
$\displaystyle-(X_{1}^{\prime}+Y)\sqrt{\min\\{(X_{1}^{\prime}+Y),a\\}\min\\{(X_{1}^{\prime}+Y),b\\}}\bigr{\\}}^{2}\mathds{1}(Y\geq
4).$

Based on the condition that
$X_{1},X_{1}^{\prime}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Binomial}(1,p)$,
the expectation of $(f-f^{1})^{2}\mathds{1}(Y\geq 4)$ conditional on $Y$ is
computed as

$\displaystyle\mathbb{E}_{X_{1},X_{1}^{\prime}}[(f-f^{1})^{2}\mathds{1}(Y\geq
4)\,|\,Y]\leavevmode\nobreak\ =\leavevmode\nobreak\ $ $\displaystyle
2p(1-p)\bigl{\\{}(Y+1)\sqrt{\min\\{(Y+1),a\\}\min\\{(Y+1),b\\}}$
$\displaystyle-Y\sqrt{\min\\{Y,a\\}\min\\{Y,b\\}}\bigr{\\}}^{2}\mathds{1}(Y\geq
4)$ $\displaystyle:=\leavevmode\nobreak\ $ $\displaystyle
2p(1-p)h(Y,a,b)\mathds{1}(Y\geq 4).$

Therefore by taking the expectation over $Y$, the second term is

$\displaystyle n(\mathrm{II})\leavevmode\nobreak\ =\leavevmode\nobreak\
2np(1-p)\mathbb{E}_{Y}[h(Y,a,b)\mathds{1}(Y\geq 4)].$

For simplicity, we denote the probability mass function of
$\mathrm{Binomial}(n-1,p)$ by $q_{n-1,p}(x)$ for $x=0,1,\ldots,n-1$. Then

$\displaystyle np\mathbb{E}_{Y}[h(Y,a,b)\mathds{1}(Y\geq
4)]\leavevmode\nobreak\ =\leavevmode\nobreak\
np\sum_{k=4}^{n-1}h(k,a,b)q_{n-1,p}(k).$

Recall that

$\displaystyle
h(y,a,b)=\bigl{[}(y+1)\sqrt{\min\\{y+1,a\\}\min\\{y+1,b\\}}-y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}\bigr{]}^{2}.$

and, by assuming $a\leq b$ without loss of generality and $y\geq 4$, bound
this function under different cases carefully. In the following analysis, we
use the fact that for $x,y\geq 4$ and some $\xi\geq 4$,

$\displaystyle(y+x)^{3/2}=y^{3/2}+\frac{3\sqrt{y}}{2}x+\frac{3}{8\sqrt{\xi+y}}x^{2}.$

Therefore, it holds that $[(y+1)^{3/2}-y^{3/2}]^{2}\lesssim y$ for all $y\geq
4$. In fact, it can be numerically verified that

$\displaystyle[(y+1)^{3/2}-y^{3/2}]^{2}\leq\underbrace{\frac{(5^{3/2}-4^{3/2})^{2}}{4}}_{:=c_{0}}y\quad\text{for
all $y\geq 4$.}$

Based on this observation, we upper bound $h(y,a,b)$ as follows:

  * •

If $a\leq b\leq y$, then

$\displaystyle h(y,a,b)=ab\leavevmode\nobreak\ \leq\leavevmode\nobreak\ $
$\displaystyle\min\\{y,a\\}\min\\{y,b\\}$
$\displaystyle\leq\leavevmode\nobreak\ $ $\displaystyle
y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $a\leq y<y+1=b$, then

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{a(y+1)}-y\sqrt{ay}]^{2}=a[(y+1)^{3/2}-y^{3/2}]^{2}\leq
c_{0}ay$ $\displaystyle\leq\leavevmode\nobreak\ $ $\displaystyle
c_{0}\min\\{y,a\\}\min\\{y,b\\}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle c_{0}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $a\leq y<y+1<b$, then

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{a(y+1)}-y\sqrt{ay}]^{2}=a[(y+1)^{3/2}-y\sqrt{y}]^{2}\leq
c_{0}ay$ $\displaystyle\leq\leavevmode\nobreak\ $ $\displaystyle
c_{0}\min\\{y,a\\}\min\\{y,b\\}$ $\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle c_{0}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $y=a<y+1=b$,

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{ab}-y\sqrt{ay}]^{2}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{y(y+1)}-y^{2}]^{2}$
$\displaystyle<\leavevmode\nobreak\ $
$\displaystyle[(y+1)^{2}-y^{2}]^{2}=(2y+1)^{2}\leq\frac{81}{16}y^{2}\quad\text{since
$y\geq 4$}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}\min\\{y,a\\}\min\\{y,b\\}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $y=a<y+1<b$,

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{a(y+1)}-y\sqrt{ay}]^{2}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle[(y+1)\sqrt{y(y+1)}-y^{2}]^{2}$
$\displaystyle<\leavevmode\nobreak\ $
$\displaystyle[(y+1)^{2}-y^{2}]^{2}=(2y+1)^{2}\leq\frac{81}{16}y^{2}\quad\text{since
$y\geq 4$.}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}\min\\{y,a\\}\min\\{y,b\\}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $y<a=y+1=b$,

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)^{2}-y^{2}]^{2}=(2y+1)^{2}<\frac{81}{16}y^{2}\quad\text{since
$y\geq 4$.}$ $\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}\min\\{y,a\\}\min\\{y,b\\}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

  * •

If $y<y+1<a\leq b$,

$\displaystyle h(y,a,b)\leavevmode\nobreak\ =\leavevmode\nobreak\ $
$\displaystyle[(y+1)^{2}-y^{2}]^{2}=(2y+1)^{2}<\frac{81}{16}y^{2}\leq
9\min\\{b^{2},y^{2}\\}\quad\text{since $y\geq 4$.}$
$\displaystyle=\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}\min\\{y,a\\}\min\\{y,b\\}$
$\displaystyle\leq\leavevmode\nobreak\ $
$\displaystyle\frac{81}{16}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

In summary, for $y\geq 4$,

$\displaystyle h(y,a,b)\lesssim\min\\{y,a\\}\min\\{y,b\\}\lesssim
y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}.$

Using this upper bound, we have

$\displaystyle np\mathbb{E}_{Y}[h(Y,a,b)\mathds{1}(Y\geq
4)]\leavevmode\nobreak\ =\leavevmode\nobreak\ $ $\displaystyle
np\sum_{k=4}^{n-1}h(k,a,b)q_{n-1,p}(k)$
$\displaystyle\overset{\mathrm{(i)}}{\lesssim}\leavevmode\nobreak\ $
$\displaystyle np\sum_{k=4}^{n-1}\min\\{y,a\\}\min\\{y,b\\}q_{n-1,p}(k)$
$\displaystyle\overset{\mathrm{(ii)}}{\lesssim}\leavevmode\nobreak\ $
$\displaystyle
np\sum_{k=4}^{n-1}y\sqrt{\min\\{y,a\\}\min\\{y,b\\}}q_{n-1,p}(k).$

The last inequality (ii) yields that

$\displaystyle n\mathrm{(II)}\leavevmode\nobreak\ \lesssim\leavevmode\nobreak$
$\displaystyle np\mathbb{E}[Y\sqrt{\min\\{Y,a\\}\min\\{Y,b\\}}\mathds{1}(Y\geq
4)]$ (16) $\displaystyle\lesssim\leavevmode\nobreak$ $\displaystyle