Limit theorems of Azadkia-Chatterjee's conditional graph correlation

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

97,697 characters

Limit theorems of Azadkia-Chatterjee's conditional graph correlation



\setlength{\abovedisplayskip}{5pt}
\setlength{\belowdisplayskip}{5pt}
\setlength{\abovedisplayshortskip}{5pt}
\setlength{\belowdisplayshortskip}{5pt}
\hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}



\title{\LARGE Limit theorems of Azadkia-Chatterjee's conditional graph correlation}

\author{Muhong Gao\thanks{School of Statistics, University of International Business and Economics, Beijing, China; e-mail: {\tt [email removed]}}, ~~Fang Han\thanks{Department of Statistics, University of Washington, Seattle, WA 98195, USA; e-mail: {\tt [email removed]}}, ~ and ~Qizhai Li\thanks{Academy of Mathematics and System Science, Chinese Academy of Sciences, Beijing, China; e-mail: {\tt [email removed]}}}

\date{\today}

\maketitle

\vspace{-1em}

\begin{abstract}
Inferring the strength of conditional dependence and testing conditional independence are fundamental problems in statistics. A recent breakthrough by Azadkia and Chatterjee introduced, for the first time, a conditional dependence measure that equals $0$ if and only if the variables under study are conditionally independent, and equals $1$ if and only if they are conditionally perfectly dependent. They further proposed a computationally efficient and strongly consistent estimator, $T_n$, based on an ingenious use of ranks and nearest neighbors. Despite these attractive features, the asymptotic theory of $T_n$ has remained largely undeveloped. This paper closes that gap. We prove that, under general dependence, $T_n$ is asymptotically normal and its limiting variance admits a closed form. We also construct consistent variance estimators that are computationally efficient and implementable in $O(n\log n)$ time. Taken together with existing bias-correction methods, these results provide a complete inferential theory for $T_n$.
\end{abstract}

{\bf Keywords}: Measure of conditional dependence, test of conditional independence, dependence measure, rank-based statistic, graph-based statistic




\section{Introduction} \label{sec:intro}
Consider the random triplet $(Y,{\boldsymbol{X}},{\boldsymbol{Z}})$, where $Y\in\mathbb{R}$ is a random scalar, and ${\boldsymbol{X}}\in\mathbb{R}^p$ and ${\boldsymbol{Z}}\in\mathbb{R}^q$ are random vectors of dimensions $p$ and $q$, respectively. Our goal is to infer the strength of conditional dependence between $Y$ and ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$, as well as to test the null hypothesis
\begin{eqnarray}
H_0:\ \text{$Y$ is conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$,}
\label{eq:null}
\end{eqnarray}
on the basis of $n$ independent copies $(Y_i,{\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)$'s of $(Y,{\boldsymbol{X}},{\boldsymbol{Z}})$.

While the problem of testing \eqref{eq:null} has been studied extensively in the literature, the problem of quantifying the strength of conditional dependence is arguably equally important. In this direction, the recent work of \cite{azadkia2019simple} constitutes a major breakthrough. Building on ideas from \cite{MR3024030} and \cite{chatterjee2020new} for measuring unconditional dependence, they introduced the following population quantity, where ${\mathrm P}_Y$ denotes the law of $Y$:


\begin{eqnarray}
T = T(Y,{\boldsymbol{X}} \mid {\boldsymbol{Z}})
=
\frac{\int {\mathrm E}\Big[\mathrm{Var}\Big\{{\mathrm P}(Y \ge y \mid {\boldsymbol{X}},{\boldsymbol{Z}})\,\big|\,{\boldsymbol{Z}}\Big\}\Big] \, {\, \mathrm{d}} {\mathrm P}_Y(y)}
{\int {\mathrm E}\Big[\mathrm{Var}\Big\{\mathbf{1}(Y \ge y)\mid {\boldsymbol{Z}}\Big\}\Big] \, {\, \mathrm{d}} {\mathrm P}_Y(y)}.
\label{eq:T}
\end{eqnarray}
Azadkia and Chatterjee proved that $T=0$ if and only if $Y$ is conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$, whereas $T=1$ if and only if $Y$ is almost surely a measurable function of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$. To the best of our knowledge, this is the \emph{first} measure of conditional dependence that captures the full range of dependence strength in this manner.

What makes the contribution of \cite{azadkia2019simple} even more striking is the accompanying statistical estimator. Let $R_i$ denote the rank of $Y_i$ among $\{Y_j\}_{j=1}^n$, let $N(i)$ index the nearest neighbor (NN) of ${\boldsymbol{Z}}_i$, and let $M(i)$ index the nearest neighbor of $({\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)$, with both nearest neighbors defined under the Euclidean metric. Azadkia and Chatterjee introduced the following rank/graph-based statistic, which we refer to as the ``Azadkia--Chatterjee conditional graph correlation'',
\begin{eqnarray}
T_n = T_n(Y,{\boldsymbol{X}} \mid {\boldsymbol{Z}})
=
\frac{\sum_{i=1}^n \bigl(\min\{R_i,R_{M(i)}\}-\min\{R_i,R_{N(i)}\}\bigr)}
{\sum_{i=1}^n \bigl(R_i-\min\{R_i,R_{N(i)}\}\bigr)},
\label{eq:T_n}
\end{eqnarray}
as a strongly consistent for $T$. Moreover, $T_n$ possesses several notable advantages:
\begin{enumerate}[label=(\roman*)]
\item it is fully nonparametric and tuning-parameter-free;
\item it completely avoids the need to estimate conditional densities, conditional characteristic functions, or mutual information;
\item it is computable in $O(n\log n)$ time.
\end{enumerate}
These features, together with the conceptual appeal of $T$, make $T_n$ an especially attractive tool for quantifying conditional dependence.

At the same time, $T_n$ has clear and important limitations. As noted in \cite{azadkia2019simple}, it generally converges to $T$ at a subparametric rate, and no limit theory has been available for $T_n$. Consequently, one cannot directly quantify the uncertainty arising from random sampling. Resolving this difficulty is by no means routine, and for years after \cite{azadkia2019simple}, the inferential theory of $T_n$ remained open.

The goal of this paper is to resolve this issue in a definitive manner. Our main contributions are threefold:
\begin{enumerate}[label=(\roman*)]
\item under conditions on the joint distribution $F = F_{{\boldsymbol{X}},Y,{\boldsymbol{Z}}}$ of $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$, we prove that $T_n$ or its bias-corrected version is asymptotically normal (Theorem~\ref{thm:CLT-main} and Corollary~\ref{cor: CLT-Tn});
\item in addition, the limiting variance of $\sqrt{n}\,T_n$ exists and admits a closed form (the same Theorem~\ref{thm:CLT-main} and Corollary~\ref{cor: CLT-Tn});
\item moreover, a consistent variance estimator exists that is computationally efficient and can be implemented in $O(n\log n)$ time (Theorem \ref{thm: est_var-main} and Proposition \ref{prop:nlogn_new}).
\end{enumerate}
Combined with the bias-correction method developed in \cite{azadkia2026biascorrection}, these results provide a complete inferential theory for $T_n$.


\subsection{Related literature}

The study of conditional dependence is intimately related to that of unconditional dependence. Indeed, when ${\boldsymbol{Z}}$ is degenerate, the Azadkia--Chatterjee conditional graph correlation $T_n(Y,{\boldsymbol{X}} \mid {\boldsymbol{Z}})$ reduces to an unconditional graph correlation between $Y$ and ${\boldsymbol{X}}$, which we denote by $\xi_n=\xi_n(Y, {\boldsymbol{X}})$.
It is therefore natural to discuss the related literature on both conditional and unconditional dependence together.

We first discuss the seminal work of \cite{chatterjee2020new} on quantifying unconditional dependence between $Y$ and $X$ (when $p=1$), together with the subsequent work of \cite{azadkia2019simple} on conditional and unconditional dependence. They have by now generated a large and steadily expanding literature, and below we try to give a brief and highly selective review.


\begin{enumerate}[label=(\roman*)]
\item The general asymptotic normality of $\xi_n$ under arbitrary dependence between ${\boldsymbol{X}}$ and $Y$ was established in \cite{Lin_Han_2025_CLT} via direct moment calculations. Subsequently, \cite{kroll2024asymptotic} and \cite{chhaibi2026martingaleapproachfluctuationsrank} developed complementary theory based on empirical-process and martingale methods, respectively, for the analysis of Chatterjee's rank correlation \citep{chatterjee2020new}.

\item From the perspective of statistical inference, \cite{Lin_Han_2024_boostrap} showed that the classical bootstrap generally fails for $\xi_n$, whereas \cite{Dette_Kroll_2025} and \cite{olivares2025powerful} established the consistency of alternative resampling procedures, namely the $m$-out-of-$n$ bootstrap and the multiplier bootstrap, under different regimes. The closed-form expression for the limiting variance of $\xi_n$ under unconditional independence was derived in \cite{Shi_Drton_Han_2024_Bernoulli} and \cite{han2024azadkia}. Large random matrix theory for matrices built from Chatterjee's rank correlation was developed in \cite{dong2025spectral}.

\item As for statistical efficiency, the (Azadkia--)Chatterjee approach has generally been found to be underpowered for testing marginal independence in regular statistical models \citep{cao2020correlations,shi2020power,Shi_Drton_Han_2024_Bernoulli}, even though it is rate-optimal for estimating the corresponding population quantity \citep{auddy2021exact,Lin_Han_2025_CLT}. See also \cite{azadkia2026kernel} for a re-examination of the kernel-based estimator of \cite{MR3024030} and a discussion of its statistical efficiency.

\item On the methodological side, the NN graph-based framework introduced in \cite{azadkia2019simple} has inspired a variety of follow-up works. These include, among many others, \cite{deb2020kernel}, \cite{huang2022kernel}, \cite{chatterjee2024kernel}, and \cite{roudaki2026kernel} on combining kernels with graph-based methods; \cite{gamboa2022global}, which extends the idea to sensitivity analysis; \cite{lin2021boosting}, which advocates incorporating multiple NNs into estimation; \cite{hormann2026azadkia}, which extends the framework to functional data; \cite{tran2024rank}, which proposes rank-based metrics for NN graph construction; and \cite{ansari2022direct} and \cite{huang2025multivariate}, which extend the setting to multivariate ${\boldsymbol{Y}}$.

\item An equally active line of research concerns the population quantity $T$ itself, as well as related alternatives in the setting of unconditional dependence. Representative contributions include \cite{strothmann2022rearranged}, \cite{bucher2024lack}, \cite{ansari2025exact}, \cite{chierichetti2025metricity}, \cite{ansari2025ordering}, \cite{fuchs2025exact}, and \cite{ansari2026quantifying}, among many others.
\end{enumerate}

Notably, the existing literature has so far focused predominantly on the setting of unconditional dependence, in which the conditioning variable ${\boldsymbol{Z}}$ is absent. To the best of our knowledge, the main exceptions are \cite{Shi_Drton_Han_2024_Bernoulli}, which studied the use of $T_n$ for testing \eqref{eq:null} within the conditional randomization test framework; \cite{huang2022kernel}, which proposed a class of conditional dependence measures by combining graph-based and kernel-based ideas; and \cite{azadkia2025new}, which introduced refined versions of $T$ and $T_n$. Even so, inferential results remain unavailable beyond the simple setting in which $Y$ is further assumed to be independent of $({\boldsymbol{X}},{\boldsymbol{Z}})$.

Concerning the task \eqref{eq:null}, the present paper is also inevitably connected to the vast literature on testing conditional independence, and the (bias-corrected) statistic $T_n$ does yield a consistent test of \eqref{eq:null}. Our work therefore also complements the broad class of nonparametric, consistent conditional independence tests developed in \cite{MR2413488}, \cite{MR3449068}, \cite{10.5555/3020548.3020641}, \cite{cai2022distribution}, and \cite{zhang2026doubly}, among many others. At the same time, it is worth noting that these methods are not designed to consistently capture conditional perfect dependence, and their implementations are typically quadratic in $n$ or more expensive.


\subsection{Technical ingredients}

The present work builds on several earlier contributions, especially \cite{Lin_Han_2025_CLT} and \cite{azadkia2026biascorrection}, which established the asymptotic normality of the unconditional version of $T_n$ and resolved the corresponding bias-correction issue, respectively. It is therefore worth clarifying more explicitly what is technically new in the current paper.

Our first main technical contribution is a central limit theorem (CLT) for the Azadkia--Chatterjee conditional correlation coefficient $T_n$ in general settings. The main difficulty here is to handle the interaction between the following two terms in the nominator of \eqref{eq:T_n}:
\[
\sum_{i=1}^n \min\{R_i,R_{M(i)}\}
\qquad\text{and}\qquad
\sum_{i=1}^n \min\{R_i,R_{N(i)}\},
\]
when $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$ is allowed to be arbitrarily dependent. This issue is in our opinion technically far more challenging than in the unconditional setting considered in \cite{Lin_Han_2025_CLT}, and took substantially additional work of ours. In fact, resolving it necessitates sharpening several results from \cite{Lin_Han_2025_CLT} and \cite{Shi_Drton_Han_2024_Bernoulli}; these improvements are highlighted in Section~\ref{sec:Lin_Han} below.

Our second main contribution is the identification of a closed-form expression for the limiting variance of $T_n$. More precisely, we show that this variance can be represented explicitly as a functional of the joint distribution of $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$; see \eqref{eq:sigma2_tTn}, \eqref{eq:han-Tn-var}, and \eqref{eq:han-Tn-var2} below. Moreover, under $H_0$ in \eqref{eq:null}, a further simplification exists; see \eqref{eq:sigma2_H0} below. These explicit characterizations in turn enable us to construct a consistent and computationally efficient variance estimator for $T_n$ with $O(n\log n)$ complexity. Relative to the variance estimator proposed in \cite{Lin_Han_2025_CLT} and the $m$-out-of-$n$ bootstrap considered in \cite{Dette_Kroll_2025}, this provides a computationally more efficient inferential tool.




\subsection{Paper organization and notation}

\paragraph*{Paper organization.}

The rest of the paper is organized as follows. Section~\ref{sec:Lin_Han} revisits the limit theorems established in \cite{Lin_Han_2025_CLT} for the Azadkia--Chatterjee \emph{unconditional} correlation coefficient, and presents our new findings related to this unconditional dependence measure. Section~\ref{sec:SI} introduces the proposed inferential framework for the Azadkia--Chatterjee conditional correlation coefficient. Section~\ref{sec:theory} develops the corresponding theory. Section~\ref{sec:simu} reports numerical experiments illustrating the finite-sample performance of the proposed procedure. All proofs are deferred to the Appendix.

\paragraph*{Notation.}

For any integer $n \ge 1$, let $\llbracket n \rrbracket = \{1,2,\dots,n\}$. A set consisting of distinct elements $x_1,\dots,x_n$ is written either as $\{x_1,\dots,x_n\}$ or as $\{x_i\}_{i=1}^n$. For a real random vector ${\boldsymbol{W}}$, let ${\mathrm P}_{{\boldsymbol{W}}}$, $F_{{\boldsymbol{W}}}$, and $\mathrm{supp}({\boldsymbol{W}})$ denote its induced probability measure, cumulative distribution function, and support, respectively. We write $\mathbf{1}(\cdot)$ for the indicator function. For a vector ${\boldsymbol{v}} \in \mathbb{R}^d$, let $\|{\boldsymbol{v}}\|$ denote its Euclidean norm. For any $a,b \in \mathbb{R}$, define $a \vee b = \max\{a,b\}$ and $a \wedge b = \min\{a,b\}$. For a finite set $A$, let $\#A$ or $|A|$ denote its cardinality. The symbols $\lfloor \cdot \rfloor$ and $\lceil \cdot \rceil$ denote the floor and ceiling functions.
For a random variable $U$ and a random vector ${\boldsymbol{V}}$, let $\mu_U$ denote the law of $U$ and let $\mu_{U \mid {\boldsymbol{V}}}$ denote the conditional law of $U$ given ${\boldsymbol{V}}$.
Finally, $\stackrel{\mathrm{a.s.}} \to$, $\stackrel{\mathrm{P}} {\to}$, and $\stackrel{\mathcal D} {\to}$ denote convergence almost surely, in probability, and in distribution, respectively. Unless stated otherwise, the terms ``absolutely continuous'' and ``almost everywhere'' are understood with respect to Lebesgue measure.


\section{Revisiting \cite{Lin_Han_2025_CLT}: CLT of the Azadkia-Chatterjee's unconditional correlation coefficient} \label{sec:Lin_Han}

This section reviews and refines the CLT for Azadkia–Chatterjee’s unconditional correlation coefficient $\xi_n(Y,{\boldsymbol{X}})$.
To facilitate comparison with $T_n$ in later sections and to maintain notational consistency throughout the paper, in this section we use ${\boldsymbol{Z}}$ in place of ${\boldsymbol{X}}$ and write $\xi_n(Y,{\boldsymbol{Z}})$ instead of $\xi_n(Y,{\boldsymbol{X}})$.

Specifically, in this section let $Y$ be a real-valued random variable and ${\boldsymbol{Z}}$ be a random vector in $\mathbb{R}^q$, both defined on the same probability space. Let $(Y_1,{\boldsymbol{Z}}_1),\dots,(Y_n,{\boldsymbol{Z}}_n)$ be $n$ independent copies of $(Y,{\boldsymbol{Z}})$. Under the assumption that $(Y,{\boldsymbol{Z}})$ is continuously distributed, Azadkia-Chatterjee's unconditional correlation coefficient proposed by \cite{azadkia2019simple} is defined as
\begin{eqnarray}
\xi_n(Y,{\boldsymbol{Z}}) := \frac{6}{n^2-1} \sum_{i=1}^n \min\{R_i, R_{N(i)}\}-\frac{2n+1}{n-1},
\label{eq:xi_n}
\end{eqnarray}
where, as before, $R_i$ denotes the rank of $Y_i$ among $\{Y_j\}_{j=1}^n$, and $N(i)$ denotes the index of the first NN of ${\boldsymbol{Z}}_i$ among $\{{\boldsymbol{Z}}_j\}_{j=1}^n$ under the Euclidean metric. Note that $\xi_n$ extends Chatterjee's original correlation coefficient \citep{chatterjee2020new} from the univariate setting $p=1$ to the multivariate setting $p\ge 1$. As $n\to\infty$, $\xi_n$ converges almost surely to the population quantity
\begin{eqnarray}
\xi(Y,{\boldsymbol{Z}}) := \frac{\int \mathrm{Var}\Big[{\mathrm E} \big\{ \mathbf{1}(Y\ge y) \mid {\boldsymbol{Z}} \big\}\Big] {\, \mathrm{d}} {\mathrm P}_Y(y)}
{\int \mathrm{Var}\big\{ \mathbf{1}(Y\ge y)\big\} {\, \mathrm{d}} {\mathrm P}_Y(y)},
\label{eq:xi}
\end{eqnarray}
which is also known as the Dette--Siburg--Stoimenov dependence measure \citep{MR3024030}.

In the special case where ${\boldsymbol{Z}}$ and $Y$ are independent, a CLT for $\xi_n$ was first established in \citet[Theorem~3.1(ii)]{Shi_Drton_Han_2024_Bernoulli}. The subsequent work of \cite{Lin_Han_2025_CLT} extended this result to the general setting in which ${\boldsymbol{Z}}$ and $Y$ may be arbitrarily dependent. We summarize the main conclusions of \cite{Lin_Han_2025_CLT} in Proposition~\ref{prop:Lin_Han} below.

\begin{proposition}[Summary of results in \cite{Lin_Han_2025_CLT}]
\label{prop:Lin_Han}
Assume that $F_{Y,{\boldsymbol{Z}}}$ is fixed and continuous.
\begin{enumerate}[label=(\roman*)]
\item The limiting variance $\sigma^2_{\xi(Y,{\boldsymbol{Z}})} := \lim_{n \to \infty} n \mathrm{Var}(\xi_n)$ exists. Moreover, $\sigma^2_{\xi(Y,{\boldsymbol{Z}})} > 0$ if and only if $Y$ is not almost surely a measurable function of ${\boldsymbol{Z}}$.
\item There exists a consistent estimator $\widetilde{\sigma}^2$ of $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$. Moreover, $\widetilde{\sigma}^2$ can be computed in $O(n^2)$ time.
\item If $Y$ is not almost surely a measurable function of ${\boldsymbol{Z}}$, then, as $n \to \infty$,
\begin{eqnarray*}
\frac{\xi_n - {\mathrm E}(\xi_n)}{\sqrt{\mathrm{Var}(\xi_n)}} \stackrel{\mathcal D} {\to} N(0,1).
\end{eqnarray*}
\end{enumerate}
\end{proposition}

\vspace{0.2cm}

Despite these results, two notable gaps remain in \cite{Lin_Han_2025_CLT}: (1) no closed-form expression for $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$ is available; and (2) the estimator $\widetilde{\sigma}^2$ in \cite{Lin_Han_2025_CLT} requires $O(n^2)$ computational time, which is substantially more demanding than the $O(n\log n)$ complexity typically associated with rank- and graph-based statistics, and thus limits its practical usefulness.

As a byproduct of developing our general inferential theory for $T_n$, we resolve both of these issues in a definitive manner. We present these improvements to \cite{Lin_Han_2025_CLT} first, before turning to the general theory of $T_n$. We hope that this presentation makes it clearer that the present paper is not merely an extension of \cite{Lin_Han_2025_CLT} to the setting of conditional dependence.




\subsection{New probabilistic results on NNGs}

To derive the closed-form expression for $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$, we begin by reviewing and establishing several probabilistic asymptotic results for nearest neighbor graphs (NNGs), which form the foundation for our subsequent analysis.

Our first result in this section concerns a sample $\{{\boldsymbol{W}}_i\}_{i=1}^n$ consisting of $n$ independent copies of a random vector ${\boldsymbol{W}} \in \mathbb{R}^d$. Let $\mathcal{G}_n$ denote the associated directed nearest-neighbor graph (NNG) with vertex set $\llbracket n \rrbracket$. A directed edge $i \to j$ is drawn between two distinct vertices $i$ and $j$ whenever ${\boldsymbol{W}}_j$ is the nearest neighbor of ${\boldsymbol{W}}_i$. Denote by $\mathcal{E}(\mathcal{G}_n)$ the edge set of $\mathcal{G}_n$.

We begin by recalling Theorem 1 of \cite{MR937563} on the expected number of mutual NN pairs.

\begin{lemma}[\cite{MR937563}, expected number of mutual NN pairs]
\label{lemma:q_d}
Assume that ${\boldsymbol{W}}$ is Lebesgue absolutely continuous. Then, for any fixed $i$, we have
\begin{eqnarray*}
{\mathrm E}\Big(\#\big\{j \in \llbracket n \rrbracket: i\to j,\ j \to i \in \mathcal{E}(\mathcal{G}_n)\big\} \ \Big | \ {\boldsymbol{W}}_i \Big)  \stackrel{\mathrm{P}} {\to} \mathfrak{q}_d,
\end{eqnarray*}
where $\mathfrak{q}_d$ is a positive constant depending only on $d$, with explicit expression
\begin{eqnarray}
\mathfrak{q}_d := \Big\{2- I_{3/4}\Big(\frac{d+1}{2}, \frac{1}{2}\Big)\Big\}^{-1},
\qquad
I_{x}(a,b) = \frac{\int_{0}^x t^{a-1} (1-t)^{b-1} {\, \mathrm{d}} t}{\int_{0}^1 t^{a-1} (1-t)^{b-1} {\, \mathrm{d}} t}.
\label{eq:q_d}
\end{eqnarray}
Consequently,
\begin{eqnarray*}
{\mathrm E}\Big( \frac{1}{n} \#\big\{(i,j) \text{ distinct}: i\to j,\ j \to i \in \mathcal{E}(\mathcal{G}_n)\big\}  \Big) \to \mathfrak{q}_d.
\end{eqnarray*}
\end{lemma}

Note that the first statement may equivalently be written as ${\mathrm P}\{N(N(i))=i \mid {\boldsymbol{W}}_i\} \stackrel{\mathrm{P}} {\to} \mathfrak{q}_d$, that is, the conditional probability that the NN of ${\boldsymbol{W}}_i$ has ${\boldsymbol{W}}_i$ itself as its NN converges to $\mathfrak{q}_d$, regardless of the specific value of ${\boldsymbol{W}}_i$.

We are now ready to present our first new result, which extends the argument of Lemma~\ref{lemma:q_d} to the conditional expectation of the number of shared nearest-neighbor triplets.

\begin{lemma}[Conditional expected number of shared-NN triplets]
\label{lemma:o_d}
Assume that ${\boldsymbol{W}}$ is Lebesgue absolutely continuous and admits a continuous density on its support. Then, for any fixed $i$, as $n \to \infty$,
\begin{eqnarray*}
{\mathrm E}\Big(\#\big\{j\in \llbracket n \rrbracket: j \to k,\ i \to k \in \mathcal{E}(\mathcal{G}_n) \big\} \ \Big | \ {\boldsymbol{W}}_i \Big) \stackrel{\mathrm{P}} {\to} \mathfrak{o}_d,
\end{eqnarray*}
where $\mathfrak{o}_d$ is a positive constant depending only on $d$, with explicit expression
\begin{eqnarray}
\mathfrak{o}_d := \int_{\Gamma_{d;2}}
\exp\Big[ - \lambda \big\{ \mathcal{B}({\boldsymbol{w}}_1, \|{\boldsymbol{w}}_1\|) \cup \mathcal{B}({\boldsymbol{w}}_2, \|{\boldsymbol{w}}_2\|)\big\}\Big]{\, \mathrm{d}}({\boldsymbol{w}}_1, {\boldsymbol{w}}_2),
\label{eq:o_d}
\\
\Gamma_{d;2} :=
\Big\{({\boldsymbol{w}}_1, {\boldsymbol{w}}_2) \in (\mathbb{R}^d)^2: \max(\|{\boldsymbol{w}}_1\|, \|{\boldsymbol{w}}_2\|) < \|{\boldsymbol{w}}_1-{\boldsymbol{w}}_2\|\Big\},
\nonumber
\end{eqnarray}
with $\mathcal{B}({\boldsymbol{w}},r)$ denoting the ball of radius $r$ centered at ${\boldsymbol{w}}$, and $\lambda(\cdot)$ denoting Lebesgue measure.
\end{lemma}

Note that, by the bounded convergence theorem, together with the well-known fact that the maximum degree of an NNG is bounded \citep{MR682809}, Lemma~\ref{lemma:o_d} immediately recovers the existing result on the unconditional expectation from \cite{MR914597}:
\begin{eqnarray}
{\mathrm E}\Big(\frac{1}{n}\#\big\{(i,j,k) \text{ distinct}: j \to k,\ i \to k \in \mathcal{E}(\mathcal{G}_n) \big\} \Big) \to \mathfrak{o}_d.
\label{eq:E_to_od}
\end{eqnarray}

Table~\ref{table1} reports the values of $\mathfrak{q}_d$ and $\mathfrak{o}_d$ for the first ten dimensions, updating the calculations reported in \citet[Table~1]{han2024azadkia}.

\begin{table}[htbp]
\centering
\caption{The first $10$ values of $\mathfrak{q}_d$ and $\mathfrak{o}_d$.
Specifically, $\mathfrak{q}_d$ is computed by numerical integration according to \eqref{eq:q_d}, whereas $\mathfrak{o}_d$ is estimated by Monte Carlo simulation based on \eqref{eq:E_to_od} with $n = 10^7$.}
\label{table1}
\begin{tabular}{ccccccccccc}
\hline
$d$ & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\
\hline
$\mathfrak{q}_d$ & 0.667 & 0.622 & 0.593 & 0.573 & 0.558 & 0.547 & 0.538 & 0.531 & 0.528 & 0.521 \\
$\mathfrak{o}_d$ & 0.500 & 0.633 & 0.709 & 0.763 & 0.805 & 0.840 & 0.871 & 0.898 & 0.923 & 0.946 \\
\hline
\end{tabular}
\end{table}

\vspace{0.5cm}

Our second result concerns a setting involving two NNGs, generated from the full sample and from a subsample, respectively. Such configurations arise repeatedly in the analysis of the conditional correlation coefficient $T_n$.

To describe this setting, consider a sample $\{{\boldsymbol{W}}_i\}_{i=1}^n$, where each ${\boldsymbol{W}}_i = ({\boldsymbol{U}}_i,{\boldsymbol{V}}_i)$ is independently drawn from the random vector ${\boldsymbol{W}}=({\boldsymbol{U}},{\boldsymbol{V}})$, with ${\boldsymbol{U}} \in \mathbb{R}^{d_1}$ and ${\boldsymbol{V}} \in \mathbb{R}^{d_2}$. Let $\mathcal{G}^{{\boldsymbol{W}}}_n$ denote the NNG associated with the full sample $\{{\boldsymbol{W}}_i\}_{i=1}^n$, and let $\mathcal{G}^{{\boldsymbol{U}}}_n$ denote the NNG associated with the subsample $\{{\boldsymbol{U}}_i\}_{i=1}^n$. Lemma~\ref{lemma:two_NNGs} below establishes the convergence of the conditional expected number of shared-NN triplets across the two graphs $\mathcal{G}^{{\boldsymbol{W}}}_n$ and $\mathcal{G}^{{\boldsymbol{U}}}_n$.

\begin{lemma}[Shared-NN triplets across two NNGs]
\label{lemma:two_NNGs}
Assume that ${\boldsymbol{W}}$ is Lebesgue absolutely continuous and admits a continuous density on its support. Then, for each fixed $i$, as $n \to \infty$,
\begin{eqnarray*}
{\mathrm E}\Big(\#\big\{j\in \llbracket n \rrbracket: j \to k \in \mathcal{E}(\mathcal{G}^{\boldsymbol{U}}_n),\ i \to k \in \mathcal{E}(\mathcal{G}^{\boldsymbol{W}}_n) \big\} \ \Big | \ {\boldsymbol{W}}_i \Big) \stackrel{\mathrm{P}} {\to} 1.
\end{eqnarray*}
\end{lemma}


Of note, the unconditional version
\begin{eqnarray*}
{\mathrm E}\Big(\frac{1}{n}\#\big\{(i,j,k) \text{ distinct}: j \to k \in \mathcal{E}(\mathcal{G}^{\boldsymbol{U}}_n),\ i \to k \in \mathcal{E}(\mathcal{G}^{\boldsymbol{W}}_n) \big\}\Big) \to 1,
\end{eqnarray*}
was previously established in \citet[Lemma 7.4]{Shi_Drton_Han_2024_Bernoulli}. As in Lemmas~\ref{lemma:q_d} and \ref{lemma:o_d}, here we show that the asymptotic behavior of the corresponding conditional expectation is invariant with respect to the specific value of ${\boldsymbol{W}}_i$.


Lemmas~\ref{lemma:q_d} and \ref{lemma:o_d} will be used to derive the closed-form expression for the limiting variance of $\xi_n$, whereas Lemma~\ref{lemma:two_NNGs} will be used to derive the closed-form expression for the limiting variance of $T_n$ in Section~\ref{sec:theory}.


\subsection{Closed-form expression for the limiting variance of $\xi_n$}

In what follows, let ${\widetilde Y}_i$, ${\widetilde Y}'_i$, and ${\widetilde Y}''_i$ denote copies of $Y_i$ such that, conditional on ${\boldsymbol{Z}}_i$, they are independently and identically distributed ($\mathrm{i.i.d.}$) according to the conditional distribution of $Y_i$ given ${\boldsymbol{Z}}_i$.


Theorem~\ref{thm:var_xi} below constitutes the main theoretical result of this section. It provides an explicit expression for the limiting variance of $\xi_n$ when ${\boldsymbol{Z}}$ and $Y$ are possibly dependent. In this way, it complements \cite{Lin_Han_2025_CLT} and further extends the corresponding results of \cite{Shi_Drton_Han_2024_Bernoulli} and \cite{chhaibi2026martingaleapproachfluctuationsrank} to the settings of dependent pairs and multivariate ${\boldsymbol{Z}}$, respectively.

\begin{theorem}[Asymptotic variance of $\xi_n$ under dependence] \label{thm:var_xi}
Assume that $F_{Y,{\boldsymbol{Z}}}$ is fixed and continuous. Assume further that ${\boldsymbol{Z}}$ is Lebesgue absolutely continuous and admits a continuous density function on its support. We then have
\begin{align}
\hspace{-0.9cm}\sigma^2_{\xi(Y,{\boldsymbol{Z}})} &:=\lim_{n \to \infty} n \mathrm{Var}\big\{\xi_n(Y,{\boldsymbol{Z}})\big\} \cr
=& \ 36 \, \Big\{	(1+\mathfrak{q}_q) \,T_1 + (2-2\mathfrak{q}_q + \mathfrak{o}_q) \, T_2 - (2-\mathfrak{q}_q + \mathfrak{o}_q) \, T_3 +4 \, T_4 - 2 \, T_5 + T_6 - 4 \, T_7 \Big\}, \label{eq:sigma2_xi}
\end{align}
where
\begin{align}
T_1 &=  {\mathrm E}\big\{ F_Y^2(Y \wedge {\widetilde Y})\big\}, \quad & T_2 &= {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y}) \cdot F_Y(Y \wedge{\widetilde Y}')\big\},  \cr
T_3 &= {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y}) \cdot F_Y({\widetilde Y}' \wedge{\widetilde Y}'')\big\}, \quad & T_4 &= {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge {\widetilde Y}_2)\cdot F_Y(Y_1 \wedge {\widetilde Y}_1)\big\},  \cr
T_5 &=  {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge {\widetilde Y}_2)\cdot F_Y({\widetilde Y}_1 \wedge {\widetilde Y}_1')\big\}, \quad & T_6 &= {\mathrm E}\big\{F_Y(Y_1 \wedge {\widetilde Y}_1 \wedge Y_2 \wedge {\widetilde Y}_2)\big\}, \cr
T_7 &= {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y})\big\}^2, & \label{eq:T1-T7}
\end{align}
with $\mathfrak{q}_q$ and $\mathfrak{o}_q$ defined in \eqref{eq:q_d} and \eqref{eq:o_d}, respectively.
\end{theorem}

Notably, when $Y$ and ${\boldsymbol{Z}}$ are further assumed to be independent, all terms $T_1$ through $T_7$ in \eqref{eq:T1-T7} reduce to distribution-free constants, and \eqref{eq:sigma2_xi} further simplifies to the corresponding expression in \cite{Shi_Drton_Han_2024_Bernoulli}, summarized in Proposition~\ref{prop:var_xi_indep} below.

\begin{proposition}[Asymptotic variance under independence] \label{prop:var_xi_indep}
If $Y$ is independent of ${\boldsymbol{Z}}$, then the terms $T_1$ through $T_7$ in \eqref{eq:T1-T7} all reduce to constants, with values $T_1 = 1/6$, $T_2 = 2/15$, $T_3=1/9$, $T_4 = 1/15$, $T_5 = 1/9$, $T_6 = 1/5$, and $T_7 = 1/9$. Consequently, the limiting variance reduces to
\begin{eqnarray*}
\frac{2}{5} +\frac{2}{5}\mathfrak{q}_q + \frac{4}{5}\mathfrak{o}_q,
\end{eqnarray*}
which agrees with the expression in \citet[Theorem 3.1(ii)]{Shi_Drton_Han_2024_Bernoulli}.
\end{proposition}

\subsection{Variance estimation for $\xi_n$}

Recall that the variance estimator proposed in \cite{Lin_Han_2025_CLT} requires $O(n^2)$ computational time. In contrast, our new Theorem~\ref{thm:var_xi} yields, as a byproduct, a new estimator of $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$ that can be computed in $O(n\log n)$ time. We first present the form of this estimator,  along with its theoretical properties, in Theorem~\ref{thm:est_var_xi} below.

\begin{theorem}[Consistent estimator of the limiting variance]\label{thm:est_var_xi}
Assume the conditions of Theorem~\ref{thm:var_xi}. Then the following statistic converges in probability to $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$:
\begin{eqnarray}
\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})} &:=&  36\Big\{	(1+\mathfrak{q}_q) \,\widehat{T}_1 + (2-2\mathfrak{q}_q + \mathfrak{o}_q) \, \widehat{T}_2 - (2-\mathfrak{q}_q + \mathfrak{o}_q) \, \widehat{T}_3 \cr
&& \qquad +4 \, \widehat{T}_4 - 2 \, \widehat{T}_5 + \widehat{T}_6 - 4 \, \widehat{T}_7 \Big\},
\label{eq:est_var_cha}
\end{eqnarray}
where
\begin{align}
\widehat{T}_1 &=  \frac{1}{n^3} \sum_{i=1}^n  \big(R_i \wedge R_{N(i)} \big)^2,  &
\widehat{T}_2 &= \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{N(i)} \big) \big( R_i \wedge R_{N_2(i)}\big),  \cr
\widehat{T}_3 &= \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{N(i)} \big) \big( R_{N_2(i)} \wedge R_{N_3(i)}\big), \qquad &
\widehat{T}_4 &= \frac{1}{n^3}
\sum_{1 \leq i \neq j \leq n}
\mathbf{1}\big(R_i \leq R_j \wedge R_{N(j)}\big) \big(R_i \wedge R_{N(i)}\big),  \cr
\widehat{T}_5 &=   \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n} \mathbf{1}\big(R_i \leq R_j \wedge R_{N(j)}\big) \big(R_{N(i)} \wedge R_{N_2(i)}\big), \hspace{-4.5em} \cr
\widehat{T}_6 &= \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n}  R_i \wedge R_{N(i)} \wedge R_j \wedge R_{N(j)},  &
\widehat{T}_7 &= \Big(\frac{1}{n^2}\sum_{i=1}^n R_i \wedge R_{N(i)}\Big)^2, \label{eq:hatT1-hatT7}
\end{align}
with $N_2(i)$ and $N_3(i)$ denoting the indices of the second and third NNs of ${\boldsymbol{Z}}_i$, respectively.
\end{theorem}

Compared with the original estimator in \cite{Lin_Han_2025_CLT} (Theorem 1.1), the new estimator $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ is notably simpler, owing to the explicit closed-form expression of $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ established in Theorem~\ref{thm:var_xi} through the incorporation of the constants $\mathfrak{q}_q$ and $\mathfrak{o}_q$. Moreover, unlike $\widetilde\sigma^2$ in \cite{Lin_Han_2025_CLT}, the following proposition shows that $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ can be computed in $O(n\log n)$ time.

\begin{proposition} \label{prop:nlogn}
The estimator $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ in \eqref{eq:est_var_cha} can be computed in $O(n \log n)$ time. In particular, the terms $\widehat{T}_4$, $\widehat{T}_5$, and $\widehat{T}_6$ admit $O(n \log n)$ implementations via Algorithm~\ref{alg1}.
\end{proposition}

\begin{algorithm}[htbp]
\caption{Fast computations of $\widehat{T}_4$, $\widehat{T}_5$, and $\widehat{T}_6$ in $O(n \log n)$ time}
\label{alg1}
\begin{algorithmic}[1]

\Require Sample $\{({\boldsymbol{Z}}_i,Y_i)\}_{i=1}^n$.


\State Perform $k$-NN search on $\{{\boldsymbol{Z}}_i\}_{i=1}^n$, for $k=1$ and $2$. Obtain $N(i), N_2(i)$, for $i=1,\dots, n$.
\State Sort $\{Y_i\}_{i=1}^n$ and compute ranks $R_i$. Obtain
$R_i, R_{N(i)}, R_{N_2(i)}$, for $i=1,\dots, n$.
\State For $i=1,\dots, n$, compute $U_i = R_i \wedge R_{N(i)}$ and $V_i=  R_{N(i)} \wedge R_{N_2(i)}$. The remaining objective is to compute $\widehat{T}_4
=\sum_{i=1}^n U_i \cdot \sum_{j \in
\llbracket n \rrbracket, j\neq i}  \mathbf{1}(R_i \leq U_j)
$, $\widehat{T}_5= \sum_{i=1}^n V_i \cdot \sum_{j \in
\llbracket n \rrbracket, j\neq i}  \mathbf{1}(R_i \leq U_j)$, and
$\widehat{T}_6 = \sum_{i\neq j} U_i \wedge U_j = 2\cdot \sum_{i=1}^n U_i \cdot \big\{\sum_{j=1}^n \mathbf{1}(U_i \leq U_j) - 1\big\}$.
\State For $i=1,\dots,n$, use binary search to find the rank of $R_i$ among $\{U_j\}_{j=1}^n$, i.e., compute $R^*_i = \#\{j \in \llbracket n \rrbracket: R_i > U_j\}$; also, find the rank of $U_i$ among $\{U_j\}_{j=1}^n$, i.e., compute $R^\#_i = \#\{j \in \llbracket n \rrbracket: U_i > U_j\}$.
\State Compute $\widehat{T}_4 = \sum_{i=1}^n U_i \cdot \big\{n-R^*_i - \mathbf{1}(R_i \leq U_i)\big\}$, $\widehat{T}_5 = \sum_{i=1}^n V_i \cdot \big\{n-R^*_i - \mathbf{1}(R_i \leq U_i)\big\}$, and $\widehat{T}_6 =2\cdot\sum_{i=1}^n U_i\cdot(n- R^\#_i -1)$.
\Ensure Estimators $\widehat{T}_4$, $\widehat{T}_5$, and $\widehat{T}_6$.

\end{algorithmic}
\end{algorithm}




\section{Statistical inference of $T_n$} \label{sec:SI}

This section introduces inferential procedures for constructing confidence intervals for $T$ in \eqref{eq:T}, as well as for testing $H_0$ in \eqref{eq:null}, based on Azadkia--Chatterjee's conditional correlation coefficient $T_n$ in \eqref{eq:T_n}. Before proceeding, we first introduce some notation and preliminary observations.

Let $({\boldsymbol{X}}_1,Y_1,{\boldsymbol{Z}}_1),\dots,({\boldsymbol{X}}_n,Y_n,{\boldsymbol{Z}}_n)$ be $n$ $\mathrm{i.i.d.}$ copies of the random triplet $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$, where $Y \in \mathbb{R}$, ${\boldsymbol{X}} \in \mathbb{R}^p$, and ${\boldsymbol{Z}} \in \mathbb{R}^q$, with $p,q \ge 1$. Recall that $T_n$ in \eqref{eq:T_n} takes the form
\begin{eqnarray}
T_n=T_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) &=&  \frac{\tau_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}})}{\kappa_n(Y, {\boldsymbol{Z}})}, \cr
\text{with}	\hspace{2cm} \tau_n=\tau_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) &:=& \frac{1}{n^2}\sum_{i=1}^n \big(\min\{R_i, R_{M(i)}\}- \min\{R_i, R_{N(i)}\}\big), \label{eq:tau_n} \\
\kappa_n=\kappa_n(Y, {\boldsymbol{Z}}) &:=& \frac{1}{n^2}\sum_{i=1}^n\big(R_i-\min\{R_i,R_{N(i)}\}\big), \nonumber
\end{eqnarray}
where $\tau_n$ and $\kappa_n$ denote, respectively, the numerator and denominator of $T_n$ after scaling by $n^{-2}$, and $M(i)$ and $N(i)$ index the NNs of $({\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)$ and ${\boldsymbol{Z}}_i$, respectively. Whenever $Y$ is not almost surely a function of ${\boldsymbol{Z}}$, $T_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}})$ converges almost surely to the conditional dependence measure in \eqref{eq:T}, expressed as
\begin{eqnarray}
T=T(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) &=&  \tau(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) / \kappa(Y, {\boldsymbol{Z}}), \cr
\text{with}	\hspace{2cm}  \tau=\tau(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) &:=& \int {\mathrm E}\big[\mathrm{Var}\big\{{\mathrm P}(Y \geq y \mid {\boldsymbol{X}}, {\boldsymbol{Z}}) \ \big | \  {\boldsymbol{Z}} \big\}\big] {\, \mathrm{d}} {\mathrm P}_Y(y), \label{eq:tau} \\
\kappa=\kappa(Y, {\boldsymbol{Z}}) &:=& \int {\mathrm E}\big\{ \mathrm{Var} \big( \mathbf{1}(Y \geq y) \mid {\boldsymbol{Z}} \big)\big\} {\, \mathrm{d}} {\mathrm P}_Y(y), \label{eq:kappa}
\end{eqnarray}
where $\tau$ and $\kappa$ denote, respectively, the numerator and denominator of $T$.

\subsection{Confidence intervals} \label{sec:CI}

Constructing confidence intervals for $T$ using $T_n$ hinges on deriving the limiting distribution of $T_n - T$, where
\begin{eqnarray}
T_n - T \ = \
\frac{\tau_n- T\cdot \kappa_n}{\kappa_n}
\ =: \  \frac{\widetilde{T}_n}{\kappa_n}.  \label{eq:T_n-T}
\end{eqnarray}
It was shown in \cite{azadkia2019simple} that
\[
\kappa_n \stackrel{\mathrm{a.s.}} \to \kappa >0
\qquad\text{and}\qquad
\widetilde{T}_n \stackrel{\mathrm{a.s.}} \to 0.
\]
Accordingly, the main challenge is to infer the limiting distribution of the numerator term $\widetilde{T}_n$.

The construction of confidence intervals for $T$ proceeds in the following three steps.

\begin{description}
\item[Step 1] (CLT). Establish a CLT for $\widetilde{T}_n$:
\begin{eqnarray}
\sqrt{n} \big\{\widetilde{T}_n - {\mathrm E}(\widetilde{T}_n) \big\}\stackrel{\mathcal D} {\to} N(0, \sigma^2), \quad \text{as } n \to \infty, \label{eq:step1-CLT}
\end{eqnarray}
and construct a consistent estimator $\widehat{\sigma}^2$ of the limiting variance $\sigma^2$.

\item[Step 2] (Bias correction, if necessary).
Let $L_n ={\mathrm E}(\widetilde{T}_n)$ be the (asymptotic) bias in \eqref{eq:step1-CLT}. If $\sqrt{n} L_n \to 0$, then it is asymptotically negligible. In that case,
\begin{eqnarray*}
\sqrt{n} \widetilde{T}_n \stackrel{\mathcal D} {\to} N(0, \sigma^2),
\end{eqnarray*}
and, by Slutsky's theorem,
\begin{eqnarray*}
\sqrt{n} \big(T_n - T \big)\stackrel{\mathcal D} {\to} N(0, \sigma^2/\kappa^2), \quad \text{as } n \to \infty.
\end{eqnarray*}

Otherwise, the bias is not negligible. In such cases, let
\[
L^{(\tau)}_n = {\mathrm E}(\tau_n) - \tau~~~ {\rm and}~~~ L^{(\kappa)}_n = {\mathrm E}(\kappa_n) - \kappa
\]
denote the biases of $\tau_n$ and $\kappa_n$, respectively.
Construct consistent bias estimators $\widehat{L}^{(\tau)}_n$ and $\widehat{L}^{(\kappa)}_n$ such that
\begin{eqnarray}
\widehat{L}^{(\tau)}_n - L^{(\tau)}_n = o_{\mathrm P}(n^{-1/2}), \qquad  \widehat{L}^{(\kappa)}_n - L^{(\kappa)}_n = o_{\mathrm P}(n^{-1/2}).\label{eq:step2-bias}
\end{eqnarray}
Then the bias-corrected conditional correlation coefficient
\begin{eqnarray*}
T_n^{\mathrm{bc}} \ : = \
\big(\tau_n - 	\widehat{L}^{(\tau)}_n\big)/\big(\kappa_n - 	\widehat{L}^{(\kappa)}_n\big)
\end{eqnarray*}
satisfies
\begin{eqnarray*}
\sqrt{n} \big(	T_n^{\mathrm{bc}} - T \big)\stackrel{\mathcal D} {\to} N(0, \sigma^2/\kappa^2), \quad \text{as } n \to \infty.
\end{eqnarray*}

\item[Step 3] (Confidence interval). A $(1-\alpha)$ confidence interval for $T$ is given by
\begin{eqnarray}
&&\mathsf{CI}_{\alpha}(T) = \Big(T_n -\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n}} \ , \ T_n +\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n} } \Big), \quad \text{if bias correction is unnecessary,} \cr
\text{or}&&\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T) = \Big(T_n^{\mathrm{bc}} -\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n}} \ , \ T_n^{\mathrm{bc}} +\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n} } \Big), \quad \text{if bias correction is necessary,}  \qquad \label{eq:step3-CI}
\end{eqnarray}
where $z_{\alpha/2}$ denotes the $(1-\alpha/2)$-quantile of the standard normal distribution.
\end{description}

\subsection{Conditional independence testing} \label{sec:test}

Further simplifications arise when the goal is to test $H_0$ in \eqref{eq:null}. Specifically, under $H_0$, we have $T=\tau=0$, so that
\[
T_n - T = \tau_n/\kappa_n.
\]
Accordingly, testing $H_0$ is equivalent to testing $\tau=0$, which can be carried out using $\tau_n$ alone.

The construction of a test of $H_0$ based on $\tau_n$ proceeds as follows.

\begin{description}
\item[Step $\mathbf{1'}$] (CLT under $H_0$). Establish a CLT for $\tau_n$ under $H_0$:
\begin{eqnarray}\label{eq:han-CI-test}
\sqrt{n} \big\{\tau_n - {\mathrm E}(\tau_n) \big\}\stackrel{\mathcal D} {\to} N(0, \sigma_0^2), \quad \text{as } n \to \infty.
\end{eqnarray}
Construct a consistent estimator $\widehat{\sigma}_0^2$ of the limiting variance $\sigma_0^2$. \footnote{Note that under $H_0$, $T=0$, so $\tau_n$ coincides with $\widetilde{T}_n$ in \eqref{eq:T_n-T}; hence \eqref{eq:han-CI-test} is a special case of \eqref{eq:step1-CLT}.}

\item[Step $\mathbf{2'}$] (Bias correction, if necessary). Let $L^{(\tau)}_n = {\mathrm E}(\tau_n) - \tau$ denote the bias of $\tau_n$.
Whenever $\sqrt{n} L_n \nrightarrow 0$, the bias is not negligible.
In such cases, construct a bias estimator $\widehat{L}^{(\tau)}_n$ such that
\begin{eqnarray}
\widehat{L}^{(\tau)}_n - L^{(\tau)}_n = o_{\mathrm P}(n^{-1/2}), \label{eq:step2'-bias}
\end{eqnarray}
analogously to \eqref{eq:step2-bias} in Step 2 above.

\item[Step $\mathbf{3'}$] (A test of $H_0$).
The resulting level-$\alpha$ test is given by
\begin{align*}
&\mathsf{T}_{\alpha} = \mathbf{1}\big(\sqrt{n} \tau_n / \widehat{\sigma} > z_\alpha\big),  && \hspace{-1.2cm} \text{if bias correction is unnecessary,} \cr
\text{or } \quad &\mathsf{T}^{\mathrm{bc}}_\alpha = \mathbf{1}\big(	 \sqrt{n} (\tau_n - 	\widehat{L}^{(\tau)}_n ) / \widehat{\sigma} > z_\alpha\big),  && \hspace{-1.2cm} \text{if bias correction is necessary.}
\end{align*}
\end{description}




\section{Theory} \label{sec:theory}

This section provides the theoretical foundation for the inferential procedures described in Section~\ref{sec:SI}. In particular,
\begin{enumerate}[label=(\roman*)]
\item For Steps 1 and $1'$, we establish a CLT, derive a closed-form expression for the limiting variance (Section~\ref{sec:CLT}), and provide consistent estimators thereof (Section~\ref{sec:est_var});
\item For Steps 2 and $2'$, we show that bias correction could be unnecessary when the combined dimension of ${\boldsymbol{X}}$ and ${\boldsymbol{Z}}$ satisfies $p+q \leq 3$; otherwise, a bias-correction procedure is justified (Section~\ref{sec:bias_correct});
\item Finally, for Steps 3 and $3'$, we establish the validity of the proposed confidence intervals and tests (Section~\ref{sec:infer_theory}).
\end{enumerate}

\subsection{CLT} \label{sec:CLT}

Before presenting the main theorems in this section, we first introduce the following assumptions.

\begin{assumption} \label{assump_4.1}
Assume that $\{({\boldsymbol{X}}_i,Y_i,{\boldsymbol{Z}}_i): i\in\llbracket n\rrbracket\}$ are $n$ independent copies of $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$.
\end{assumption}

\begin{assumption}\label{assump_4.2}
The joint cumulative distribution function $F_{{\boldsymbol{X}},Y,{\boldsymbol{Z}}}$ of $({\boldsymbol{X}},Y,{\boldsymbol{Z}})$ is continuous.
\end{assumption}

\begin{assumption} \label{assump_4.3}
$({\boldsymbol{X}},{\boldsymbol{Z}})$ is absolutely continuous and admits a density function $f_{{\boldsymbol{X}},{\boldsymbol{Z}}}({\boldsymbol{x}},{\boldsymbol{z}})$ that is continuous on its support.
\end{assumption}

\begin{assumption}\label{assump_4.4}
$Y$ is not almost surely equal to a function of ${\boldsymbol{Z}}$.
\end{assumption}

\begin{assumption}\label{assump_4.5}
Define
$G_{\boldsymbol{z}}(t) = {\mathrm E}[\mathbf{1}(Y \geq t) \mid {\boldsymbol{Z}} = {\boldsymbol{z}}]$ and $G_{{\boldsymbol{x}},{\boldsymbol{z}}}(t) = {\mathrm E}[\mathbf{1}(Y \geq t) \mid ({\boldsymbol{X}}, {\boldsymbol{Z}}) = ({\boldsymbol{x}},{\boldsymbol{z}})]$.
For any fixed $t$, assume that the mapping ${\boldsymbol{z}} \mapsto G_{{\boldsymbol{z}}}(t)$ is continuous almost everywhere on $\mathrm{supp}({\boldsymbol{Z}})$, and that the mapping $({\boldsymbol{x}},{\boldsymbol{z}}) \mapsto G_{{\boldsymbol{x}},{\boldsymbol{z}}}(t)$ is continuous almost everywhere on $\mathrm{supp}(({\boldsymbol{X}},{\boldsymbol{Z}}))$.
\end{assumption}

Recall that $T_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) = \tau_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) / \kappa_n(Y, {\boldsymbol{Z}})$ in \eqref{eq:T_n}, with $\tau_n$ and $\kappa_n$ given by
\begin{align}
&	\tau_n(Y, {\boldsymbol{X}} \mid {\boldsymbol{Z}}) = \widetilde{\xi}_{1,n} -  \widetilde{\xi}_{2,n},     &\kappa_n(Y, {\boldsymbol{Z}}) 	= (n+1)/(2n)- 3^{-1} -  \widetilde{\xi}_{2,n}, \cr
\text{where} \hspace{1cm} & \widetilde{\xi}_{1,n} := \frac{1}{n^2}\sum_{i=1}^n \min\{R_i, R_{M(i)}\} -\frac{1}{3},
& \widetilde{\xi}_{2,n} := \frac{1}{n^2} \sum_{i=1}^n \min\{R_i, R_{N(i)}\} - \frac{1}{3}. \label{eq:txi}
\end{align}
It is worth noting that $\widetilde{\xi}_{1,n}$ and $\widetilde{\xi}_{2,n}$ can be viewed as unnormalized versions of the Azadkia--Chatterjee unconditional correlation coefficient $\xi_n$ in \eqref{eq:xi_n}, in the sense that
\begin{eqnarray}
\widetilde{\xi}_{1,n} &=& 6^{-1}\cdot \xi_n(Y, ({\boldsymbol{X}}, {\boldsymbol{Z}})) + O_{\mathrm P}(n^{-2}) +O(n^{-1}), \cr \widetilde{\xi}_{2,n} &=&  6^{-1} \cdot \xi_n(Y,{\boldsymbol{Z}}) + O_{\mathrm P}(n^{-2}) + O(n^{-1}). \label{eq:txi_2}
\end{eqnarray}
Using this notation, $\widetilde{T}_n$ in \eqref{eq:T_n-T} admits the decomposition
\begin{eqnarray}
\widetilde{T}_n \ = \kappa_n\cdot (T_n -T) \ = \ \widetilde{\xi}_{1,n} -(1 -T) \cdot \widetilde{\xi}_{2,n} - \big\{(n+1)/(2n)-3^{-1}\big\}\cdot T.  \label{eq:tTn}
\end{eqnarray}

We first derive the general CLT. As noted earlier in Section~\ref{sec:CI}, the CLT for $T_n$ relies on that for $\widetilde{T}_n$. We therefore begin by establishing a CLT for $\widetilde{T}_n$ in the general case where $Y$ may depend on ${\boldsymbol{X}}$ conditionally on ${\boldsymbol{Z}}$. Throughout Section~\ref{sec:theory}, $\overline{Y}$, $\overline{Y}'$, ${\widetilde Y}$, and ${\widetilde Y}'$ denote copies of $Y$ such that, conditional on $({\boldsymbol{X}},{\boldsymbol{Z}})$, (i) they are mutually independent, (ii) $\overline{Y},\overline{Y}' \sim \mu_{Y \mid {\boldsymbol{X}},{\boldsymbol{Z}}}$, and (iii) ${\widetilde Y},{\widetilde Y}' \sim \mu_{Y \mid {\boldsymbol{Z}}}$.

\begin{theorem}[CLT of $\widetilde{T}_n$] \label{thm:CLT-main}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Then, as $n \to \infty$, $\widetilde{T}_n \ = \kappa_n\cdot (T_n -T)$ satisfies the CLT
\begin{eqnarray*}
\sqrt{n} \big\{\widetilde{T}_n - {\mathrm E}(\widetilde{T}_n) \big\}\stackrel{\mathcal D} {\to} N(0, \sigma^2),
\end{eqnarray*}
where
\begin{eqnarray}
\sigma^2 = \lim_{n \to \infty} n \mathrm{Var}(\widetilde{T}_n) = \sigma_1^2 + (1-T)^2 \cdot \sigma_2^2 - 2\cdot (1-T) \cdot \sigma_{1,2}.  \label{eq:sigma2_tTn}
\end{eqnarray}
The explicit expressions of $\sigma_1^2$, $\sigma_2^2$, and $\sigma_{1,2}$ are given as:
\begin{eqnarray*}
&&	\sigma_1^2 = \lim_{n \to \infty} n \mathrm{Var}(\widetilde{\xi}_{1,n}) = \lim_{n \to \infty} n \mathrm{Var}\big(6^{-1} \cdot \xi_n(Y,{\boldsymbol{Z}}) \big) = 36^{-1}\cdot \sigma^2_{\xi(Y,{\boldsymbol{Z}})}, \cr
&&	\sigma_2^2 = \lim_{n \to \infty} n \mathrm{Var}(\widetilde{\xi}_{2,n}) = \lim_{n \to \infty} n \mathrm{Var}\big(6^{-1} \cdot \xi_n(Y,({\boldsymbol{X}},{\boldsymbol{Z}})) \big) = 36^{-1}\cdot \sigma^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))},
\end{eqnarray*}
where $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$ and $\sigma^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))}$ are the limiting variances of the Azadkia--Chatterjee unconditional coefficients $\xi(Y,{\boldsymbol{Z}})$ and $\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))$, respectively.
\footnote{The closed form of $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$ is provided in \eqref{eq:sigma2_xi}–\eqref{eq:T1-T7}, and $\sigma^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))}$ is defined in the same manner as $\sigma^2_{\xi(Y,{\boldsymbol{Z}})}$, with ${\boldsymbol{Z}}$ replaced by $({\boldsymbol{X}},{\boldsymbol{Z}})$.}  Furthermore,
\begin{eqnarray}
\sigma_{1,2} &=& \lim_{n \to \infty} n \cdot {\rm Cov}\big(\widetilde{\xi}_{1,n} \, , \, \widetilde{\xi}_{2,n}\big) \cr
&=& 4 \, U_1 -2\, U_2 - U_3 +2 \, U_4 - U_5 + 2\, U_6 - U_7 + U_8 - 4 \, U_9,
\label{eq:sigma12}
\end{eqnarray}
with
\begin{align}
U_1 &=  {\mathrm E}\big\{F_Y(Y \wedge \overline{Y}) \cdot F_Y(Y \wedge{\widetilde Y})\big\}, \quad & U_2 &= {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y}) \cdot F_Y(\overline{Y} \wedge\overline{Y}')\big\},  \cr
U_3 &= {\mathrm E}\big\{F_Y(Y \wedge \overline{Y}) \cdot F_Y({\widetilde Y} \wedge{\widetilde Y}')\big\}, \quad & U_4 &= {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge {\widetilde Y}_2)\cdot F_Y(Y_1 \wedge \overline{Y}_1)\big\},  \cr
U_5 &=  {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge {\widetilde Y}_2)\cdot F_Y(\overline{Y}_1 \wedge \overline{Y}_1')\big\}, \quad & U_6 &=  {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge \overline{Y}_2)\cdot F_Y(Y_1 \wedge {\widetilde Y}_1)\big\}, \cr
U_7 &=  {\mathrm E}\big\{\mathbf{1}(Y_1 \leq Y_2 \wedge \overline{Y}_2)\cdot F_Y({\widetilde Y}_1 \wedge {\widetilde Y}_1')\big\}, \quad & U_8 &= {\mathrm E}\big\{F_Y(Y_1 \wedge \overline{Y}_1 \wedge Y_2 \wedge {\widetilde Y}_2)\big\}, \cr
U_9 &= {\mathrm E}\big\{F_Y(Y \wedge \overline{Y})\big\}\cdot {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y})\big\}. & \label{eq:U1-U9}
\end{align}
\end{theorem}

By Slutsky's theorem, Theorem~\ref{thm:CLT-main} directly yields the CLT for $T_n$ and $T_n^{\mathrm{bc}}$, stated in Corollary~\ref{cor: CLT-Tn} below.

\begin{corollary}[CLT of $T_n$ and $T_n^{\mathrm{bc}}$] \label{cor: CLT-Tn}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Let $L_n ={\mathrm E}(\widetilde{T}_n)$ denote the bias of $\widetilde{T}_n$.
\begin{enumerate}[label=(\roman*)]
\item If $\sqrt{n} L_n \to 0$ as $n \to \infty$, then
\begin{eqnarray}\label{eq:han-Tn-var}
\sqrt{n} \big(T_n - T \big)\stackrel{\mathcal D} {\to} N(0, \sigma^2/\kappa^2).
\end{eqnarray}
\item Assume that the bias-corrected correlation coefficient $T_n^{\mathrm{bc}}  =
\big(\tau_n - 	\widehat{L}^{(\tau)}_n\big)/\big(\kappa_n - 	\widehat{L}^{(\kappa)}_n\big)$
satisfies \eqref{eq:step2-bias}. Then, as $n \to \infty$,
\begin{eqnarray}\label{eq:han-Tn-var2}
\sqrt{n} \big(T_n^{\mathrm{bc}}   - T \big)\stackrel{\mathcal D} {\to} N(0, \sigma^2/\kappa^2).
\end{eqnarray}
\end{enumerate}
\end{corollary}

\begin{remark}
As will be shown in Theorem~\ref{thm:bias_correct} below, the condition $p+q \leq 3$ could imply $\sqrt{n}L_n \to 0$. In this case, $\sqrt{n}(T_n-T)$ is asymptotically normal, so that $T_n$ converges to $T$ at the parametric rate $n^{-1/2}$ without the need for bias correction.
\end{remark}

We next derive the CLT under $H_0$. To this end, only the limiting distribution of $\tau_n$ is needed.

\begin{corollary}[CLT of $\tau_n$ under conditional independence]
\label{cor: CLT-H0}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Assume further that $Y$ is conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$.
Then, as $n \to \infty$, $\tau_n = \widetilde{T}_n$ satisfies the CLT
\begin{eqnarray*}
\sqrt{n} \big\{\tau_n - {\mathrm E}(\tau_n) \big\}\stackrel{\mathcal D} {\to} N(0, \sigma_0^2),
\end{eqnarray*}
where $\sigma_0^2$ is strictly positive and admits the simplified expression
\begin{eqnarray}
\sigma_0^2 &=& (2+\mathfrak{q}_q+ \mathfrak{q}_{p+q}) \cdot{\mathrm E}\big\{ F_Y^2(Y \wedge {\widetilde Y})\big\}
\cr
&&+	(\mathfrak{o}_q+ \mathfrak{o}_{p+q} - 2\mathfrak{q}_q- 2\mathfrak{q}_{p+q} -4) \cdot {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y}) \cdot F_Y(Y \wedge{\widetilde Y}')\big\}
\cr
&&+(2+\mathfrak{q}_q
+ \mathfrak{q}_{p+q}-\mathfrak{o}_q- \mathfrak{o}_{p+q}) \cdot {\mathrm E}\big\{F_Y(Y \wedge {\widetilde Y}) \cdot F_Y({\widetilde Y}' \wedge{\widetilde Y}'')\big\}, \label{eq:sigma2_H0}
\end{eqnarray}
where, as before, $\mathfrak{q}_d$ and $\mathfrak{o}_d$ are defined in \eqref{eq:q_d} and \eqref{eq:o_d}, respectively.
\end{corollary}


\subsection{Estimation of the limiting variance} \label{sec:est_var}

We begin with the general case. In view of Theorem~\ref{thm:CLT-main}, we can construct a consistent estimator of $\sigma^2$ in a manner analogous to that of Theorem~\ref{thm:est_var_xi}.

\begin{theorem}[Consistent estimator of limiting variance $\sigma^2$] \label{thm: est_var-main}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Then $\sigma^2$ in Theorem~\ref{thm:CLT-main} admits the consistent estimator
\begin{eqnarray}
\widehat{\sigma}^2 = \widehat{\sigma}_1^2 + (1-T_n)^2 \cdot \widehat{\sigma}_2^2 - 2\cdot (1-T_n) \cdot \widehat{\sigma}_{1,2}, \label{eq:sigma2_est}
\end{eqnarray}
where the explicit expressions of $\widehat{\sigma}_1^2$, $\widehat{\sigma}_2^2$, and $\widehat{\sigma}_{1,2}$ are given as follows:
\begin{eqnarray*}
\widehat{\sigma}_1^2 = 36^{-1}\cdot\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})},  \quad \text{and} \quad \widehat{\sigma}_2^2  = 36^{-1}\cdot\widehat{\sigma}^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))}.
\end{eqnarray*}
Here $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ and $\widehat{\sigma}^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))}$ are consistent estimators of the limiting variances of the Azadkia--Chatterjee unconditional correlation coefficients $\xi(Y,{\boldsymbol{Z}})$ and $\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))$, respectively.
\footnote{Here the form of $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$ is provided in \eqref{eq:est_var_cha} in Theorem~\ref{thm:est_var_xi}, and $\widehat{\sigma}^2_{\xi(Y,({\boldsymbol{X}},{\boldsymbol{Z}}))}$ is defined in the same manner as $\widehat{\sigma}^2_{\xi(Y,{\boldsymbol{Z}})}$, with ${\boldsymbol{Z}}$ replaced by $({\boldsymbol{X}},{\boldsymbol{Z}})$.}
Furthermore, we introduce
\begin{eqnarray*}
\widehat{\sigma}_{1,2} = 4 \, \widehat{U}_1 -2\, \widehat{U}_2 - \widehat{U}_3 +2 \, \widehat{U}_4 - \widehat{U}_5 + 2\, \widehat{U}_6 - \widehat{U}_7 + \widehat{U}_8 - 4 \, \widehat{U}_9,
\end{eqnarray*}
where
\begin{align*}
\widehat{U}_1 &=  \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{M(i)} \big) \big( R_i \wedge R_{N(i)}\big),  &
\widehat{U}_2 &= \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{N(i)} \big) \big( R_{M(i)} \wedge R_{M_2(i)}\big),  \cr
\widehat{U}_3 &= \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{M(i)} \big) \big( R_{N(i)} \wedge R_{N_2(i)}\big), \qquad &
\widehat{U}_4 &= \frac{1}{n^3}
\sum_{1 \leq i \neq j \leq n}
\mathbf{1}\big(R_i \leq R_j \wedge R_{N(j)}\big) \big(R_i \wedge R_{M(i)}\big),  \cr
\widehat{U}_5 &=   \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n} \mathbf{1}\big(R_i \leq R_j \wedge R_{N(j)}\big) \big(R_{M(i)} \wedge R_{M_2(i)}\big), \hspace{-4.5em} \cr
\widehat{U}_6 &=   \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n} \mathbf{1}\big(R_i \leq R_j \wedge R_{M(j)}\big) \big(R_i \wedge R_{N(i)}\big), \hspace{-4.5em} \cr
\widehat{U}_7 &=   \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n} \mathbf{1}\big(R_i \leq R_j \wedge R_{M(j)}\big) \big(R_{N(i)} \wedge R_{N_2(i)}\big), \hspace{-4.5em} \cr
\widehat{U}_8 &= \frac{1}{n^3} \sum_{1 \leq i \neq j \leq n}  R_i \wedge R_{M(i)} \wedge R_j \wedge R_{N(j)},  &
\widehat{U}_9 &= \Big(\frac{1}{n^2}\sum_{i=1}^n R_i \wedge R_{M(i)}\Big) \Big(\frac{1}{n^2}\sum_{i=1}^n R_i \wedge R_{N(i)}\Big), \label{eq:hatT1-hatT7+}
\end{align*}
with $M_2(i)$ and $N_2(i)$ denoting the indices of the second NNs of $({\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)$ and ${\boldsymbol{Z}}_i$, respectively.
\end{theorem}

According to Theorem~\ref{thm:est_var_xi} and Proposition~\ref{prop:nlogn}, the computational complexities of $\widehat{\sigma}_1^2$ and $\widehat{\sigma}_2^2$ are both of order $O(n \log n)$. For $\widehat{\sigma}_{1,2}$, an analysis similar to that in Proposition \ref{prop:nlogn} shows that its computational complexity is also of order $O(n \log n)$. Indeed, the terms $\widehat U_4$--$\widehat U_8$ involved in $\widehat{\sigma}_{1,2}$ can be computed via fast algorithms analogous to Algorithm~\ref{alg1}, so that each term can be evaluated in $O(n \log n)$ time. Consequently, the overall computational complexity of $\widehat{\sigma}^2$ is of order $O(n \log n)$.

\begin{proposition} \label{prop:nlogn_new}
The estimator $\widehat{\sigma}^2$ in \eqref{eq:sigma2_est} can be computed in $O(n \log n)$ time.
In particular, the terms $\widehat U_4$--$\widehat U_8$ can be computed in $O(n \log n)$ time via fast algorithms analogous to Algorithm~\ref{alg1}.
\end{proposition}

Next, for conditional independence testing, it suffices to estimate $\sigma_0^2$ in \eqref{eq:sigma2_H0}. To this end, we consider two alternative estimators: (1) the fast simplified estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$, and (2) the $m$-out-of-$n$ bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$. Compared with $\widehat{\sigma}^2$ in Theorem~\ref{thm: est_var-main}, both alternative estimators remain consistent under $H_0$, while offering simpler computation and improved estimation accuracy.

We first discuss the direct estimation approach based on Theorem~\ref{thm: est_var-main}, which has time complexity $O(n\log n)$.

\begin{corollary}[Consistency of fast simplified estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$]
\label{cor: est_var-H0}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Assume that $Y$ is conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$. Then $\sigma_0^2$ in Corollary~\ref{cor: CLT-H0} admits the simplified consistent estimator
\begin{eqnarray}
\widehat{\sigma}^2_{0,\mathrm{F}} &=& (2+\mathfrak{q}_q+ \mathfrak{q}_{p+q})\cdot  \frac{1}{n^3} \sum_{i=1}^n  \big(R_i \wedge R_{N(i)} \big)^2 \cr
&& + \	(\mathfrak{o}_q+ \mathfrak{o}_{p+q} - 2\mathfrak{q}_q- 2\mathfrak{q}_{p+q} -4) \cdot \frac{1}{n^3} \sum_{i=1}^n \big(R_i \wedge R_{N(i)} \big) \big( R_i \wedge R_{N_2(i)}\big)  \cr
&&+ \ (2+\mathfrak{q}_q
+ \mathfrak{q}_{p+q}-\mathfrak{o}_q- \mathfrak{o}_{p+q})\cdot \frac{1}{n^3}\sum_{i=1}^n \big(R_i \wedge R_{N(i)} \big) \big( R_{N_2(i)} \wedge R_{N_3(i)}\big),
\label{eq:sigma2_F}
\end{eqnarray}
with $N_2(i)$ and $N_3(i)$ denoting the indices of the second and third NNs of ${\boldsymbol{Z}}_i$, respectively.
\end{corollary}

We next consider the $m$-out-of-$n$ bootstrap procedure proposed in \cite{Dette_Kroll_2025}, which has computational complexity $O(B \, m \log m)$, where $B$ denotes the number of bootstrap replicates. The procedure is as follows. For each bootstrap iteration $b = 1,\dots,B$, draw $m<n$ observations without replacement from $\{(Y_i, {\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)\}_{i=1}^n$, denoted by $\{(Y_{b,j}^*, {\boldsymbol{X}}_{b,j}^*, {\boldsymbol{Z}}_{b,j}^*)\}_{j=1}^m$, and compute the statistic $\tau_m$ in \eqref{eq:tau_n} based on this bootstrap sample, denoted by $\tau_{m,b}^*$. The bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$ of $\sigma^2$ is then defined by
\begin{eqnarray}
\widehat{\sigma}^2_{0,\mathrm{B}} = \frac{m}{B}\sum_{b=1}^B \Big(\tau_{m,b}^* - \frac{1}{B}\sum_{j=1}^B \tau_{m,j}^*\Big)^2.
\label{eq:sigma2_B}
\end{eqnarray}

\begin{proposition}[Consistency of $m$-out-of-$n$ bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$]
\label{prop:est_var-B}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}.
Assume that $Y$ is conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$.
Assume $B\to \infty$, $m \to \infty$, and $m = o(n)$. Then under the null hypothesis, $\widehat{\sigma}^2_{0,\mathrm{B}}$ is consistent, in the sense that $\widehat{\sigma}^2_{0,\mathrm{B}} \stackrel{\mathrm{P}} {\to} \sigma_0^2$ as $n\to \infty$.
\end{proposition}


\subsection{Bias correction} \label{sec:bias_correct}

Recall that $L^{(\tau)}_n = {\mathrm E}(\tau_n) - \tau$ and $L^{(\kappa)}_n = {\mathrm E}(\kappa_n) - \kappa$ represent the biases of $\tau_n$ and $\kappa_n$, respectively.
The goal of this subsection is to introduce consistent estimators $\widehat{L}^{(\tau)}_n$ and $\widehat{L}^{(\kappa)}_n$ such that $\widehat{L}^{(\tau)}_n - L^{(\tau)}_n= o_{\mathrm P}(n^{-1/2})$ and $\widehat{L}^{(\kappa)}_n - L^{(\kappa)}_n= o_{\mathrm P}(n^{-1/2})$, as required in \eqref{eq:step2-bias} and \eqref{eq:step2'-bias} for our inferential procedures.

Recall from \eqref{eq:txi} that $\tau_n$ and $\kappa_n$ admit the decompositions
\begin{align*}
&	\tau_n = \widetilde{\xi}_{1,n} -  \widetilde{\xi}_{2,n},     &\kappa_n	= (n+1)/(2n)- 3^{-1} -  \widetilde{\xi}_{2,n}, \cr
\text{where} \hspace{2cm} & \widetilde{\xi}_{1,n} := \sum_{i=1}^n \min\{R_i, R_{M(i)}\} -\frac{1}{3},
& \widetilde{\xi}_{2,n} := \sum_{i=1}^n \min\{R_i, R_{N(i)}\} - \frac{1}{3}.
\end{align*}
As noted earlier in \eqref{eq:txi_2}, $\widetilde{\xi}_{1,n}$ and $\widetilde{\xi}_{2,n}$ can be viewed as unnormalized versions of the Azadkia--Chatterjee unconditional correlation coefficient $\xi_n$ in \eqref{eq:xi_n}.
Therefore, as $n\to \infty$, we have
\begin{eqnarray*}
&&	 \hspace{1cm}\widetilde{\xi}_{1,n} \stackrel{\mathrm{a.s.}} \to \widetilde{\xi}_1, \qquad \widetilde{\xi}_{2,n} \stackrel{\mathrm{a.s.}} \to \widetilde{\xi}_2, \cr
\text{where}  \hspace{1cm}  &&
\widetilde{\xi}_1  := \int \mathrm{Var}\big[{\mathrm E} \big\{ \mathbf{1}(Y\geq y) \mid {\boldsymbol{X}}, {\boldsymbol{Z}}\big\}\big] {\, \mathrm{d}} {\mathrm P}_Y(y), \cr
&& \widetilde{\xi}_2 := \int \mathrm{Var}\big[{\mathrm E} \big\{ \mathbf{1}(Y\geq y) \mid {\boldsymbol{Z}} \big\}\big] {\, \mathrm{d}} {\mathrm P}_Y(y),
\end{eqnarray*}
which correspond to the numerator of $\xi$ in \eqref{eq:xi}, with ${\boldsymbol{X}}$ therein replaced by $({\boldsymbol{X}},{\boldsymbol{Z}})$ and ${\boldsymbol{Z}}$, respectively.
It is also straightforward to verify that the population quantities $\tau$ and $\kappa$ in \eqref{eq:tau} and \eqref{eq:kappa} admit analogous decompositions, namely, $\tau = \widetilde{\xi}_1 - \widetilde{\xi}_2$ and $\kappa = 6^{-1} - \widetilde{\xi}_2$.

Denote the biases of $\widetilde{\xi}_{1,n}$ and $\widetilde{\xi}_{2,n}$ by
\begin{eqnarray*}
L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}: = {\mathrm E}(\widetilde{\xi}_{1,n}) - \widetilde{\xi}_1, \qquad L_n^{{\boldsymbol{Z}}}: = {\mathrm E}(\widetilde{\xi}_{2,n}) - \widetilde{\xi}_2.
\end{eqnarray*}
From the above decompositions, it is readily verified that
\begin{eqnarray*}
&&L^{(\tau)}_n  \ = \ {\mathrm E}(\tau_n) - \tau  \ = \ {\mathrm E}( \widetilde{\xi}_{1,n} -  \widetilde{\xi}_{2,n}) - (\widetilde{\xi}_1 - \widetilde{\xi}_2) \ = \ L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} - L_n^{{\boldsymbol{Z}}}, \cr
&&L^{(\kappa)}_n \ = \ {\mathrm E}(\kappa_n) - \kappa \ = \ - \big\{ {\mathrm E}(  \widetilde{\xi}_{2,n})   -\widetilde{\xi}_2 \big\} + 1/(2n) \ = \ - L_n^{{\boldsymbol{Z}}} + O(n^{-1}).
\end{eqnarray*}
Therefore, it suffices to perform bias correction separately for $L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}$ and $L_n^{{\boldsymbol{Z}}}$; that is, to construct estimators $\widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}$ and $\widehat{L}_n^{{\boldsymbol{Z}}}$ such that $\widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} =L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} + o_{\mathrm P}(n^{-1/2})$ and $\widehat{L}_n^{{\boldsymbol{Z}}} =L_n^{{\boldsymbol{Z}}} + o_{\mathrm P}(n^{-1/2})$. Then
\begin{eqnarray}
\widehat{L}^{(\tau)}_n = \widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} - \widehat{L}_n^{{\boldsymbol{Z}}}, \qquad \text{and} \quad \widehat{L}^{(\kappa)}_n = - \widehat{L}_n^{{\boldsymbol{Z}}},   \label{eq:bias_est}
\end{eqnarray}
serve as the desired bias estimators satisfying \eqref{eq:step2-bias} and \eqref{eq:step2'-bias}.

Note that the construction of $L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}$ is entirely analogous to that of $L_n^{{\boldsymbol{Z}}}$: one simply replaces the data $\{{\boldsymbol{Z}}_i\}_{i=1}^n$ by $\{({\boldsymbol{X}}_i,{\boldsymbol{Z}}_i)\}_{i=1}^n$.
We therefore focus on the construction of $L_n^{{\boldsymbol{Z}}}$ below to illustrate the bias-correction procedure and its theoretical justification, following the framework of \cite{azadkia2026biascorrection}.

Recall the function $G_{\boldsymbol{z}}(t) = {\mathrm E}\{\mathbf{1}(Y\geq t) \mid {\boldsymbol{Z}}={\boldsymbol{z}}\}$ defined in Assumption~\ref{assump_4.5}.
Section 3 of \cite{azadkia2026biascorrection} establishes the following alternative representation of the bias:
\begin{eqnarray}
L_n^{\boldsymbol{Z}} &=& \int {\mathrm E} \Big\{G_{{\boldsymbol{Z}}_1}(t) G_{{\boldsymbol{Z}}_{N(1)}}(t) - G^2_{{\boldsymbol{Z}}_1}(t)\Big\} {\, \mathrm{d}} {\mathrm P}_Y(t) \cr
&=&
{\mathrm E} \Big\{G_{{\boldsymbol{Z}}_1}(Y^*) G_{{\boldsymbol{Z}}_{N(1)}}(Y^*) - G^2_{{\boldsymbol{Z}}_1}(Y^*)\Big\}, \label{eq:bias}
\end{eqnarray}
where $Y^*$ is independent of ${\boldsymbol{Z}}$ and has the same marginal distribution as $Y$.
This identity naturally suggests a two-step procedure for constructing an estimator.
First, construct an appropriate estimator $\widehat{G}.(\cdot)$ of the bivariate regression function $G.(\cdot)$.
Second, approximate the expectation in \eqref{eq:bias} by replacing the population mean with the empirical distribution of $Y^*$ and ${\boldsymbol{Z}}$, leading to the estimator
\begin{eqnarray*}
\widehat{L}_n^{{\boldsymbol{Z}}} = \frac{1}{n(n-1)} \sum_{1 \leq i \neq j \leq n} \Big\{ \widehat{G}_{{\boldsymbol{Z}}_i}(Y_j) \widehat{G}_{{\boldsymbol{Z}}_{N(i)}}(Y_j) - \widehat{G}^2_{{\boldsymbol{Z}}_i}(Y_j) \Big\}.
\end{eqnarray*}
Note that, for each fixed $t$, $G_{\boldsymbol{Z}}(t) = {\mathrm E}\{\mathbf{1}(Y\geq t) \mid {\boldsymbol{Z}}\}$ is the regression mean function of $\mathbf{1}(Y\geq t)$ on ${\boldsymbol{Z}}$. This motivates estimating $G.(t)$ by regression techniques.
According to the results in Section 4.2.2 of \cite{azadkia2026biascorrection}, $G.(t)$ can be effectively estimated by ridge least squares \citep{tuo2024asymptotic}, which yields the estimator $\widehat{G}.(t)$ and hence $\widehat{L}_n^{{\boldsymbol{Z}}}$. The complete procedure for computing $\widehat{L}_n^{{\boldsymbol{Z}}}$ is summarized in Algorithm~\ref{alg2}.

\begin{algorithm}[htbp]
\caption{Compute bias estimator $\widehat{L}_n^{{\boldsymbol{Z}}}$}
\label{alg2}
\begin{algorithmic}[1]

\Require Sample $\{(Y_i, {\boldsymbol{Z}}_i)\}_{i=1}^n$; $K$ basis functions ${\boldsymbol{p}}({\boldsymbol{z}}) = (p_1({\boldsymbol{z}}),...,p_K({\boldsymbol{z}}))^\top$, defined for ${\boldsymbol{z}} \in \mathbb{R}^q$; regularization penalty parameter $\lambda_n>0$.
\State
For each $j = 1, \dots, n$, let $t= Y_j$. Solve the ridge estimator
\begin{eqnarray*}
\widehat{{\boldsymbol{\beta}}}_j = \arg \min_{{\boldsymbol{\beta}} \in \mathbb{R}^K} \Big\{\frac{1}{n}\sum_{i=1}^n  \big\{ \mathbf{1}(Y_i \geq t) - {\boldsymbol{p}}({\boldsymbol{Z}}_i)^\top {\boldsymbol{\beta}}\big\} + \lambda_n \|{\boldsymbol{\beta}}\|^2 \Big\}.
\end{eqnarray*}
Obtain $\widehat{G}_{\boldsymbol{z}}(Y_j) = \widehat{{\boldsymbol{\beta}}}_j^\top {\boldsymbol{p}}({\boldsymbol{z}})$, for any given ${\boldsymbol{z}} \in \mathbb{R}^q$.
\State Perform NN search on $\{{\boldsymbol{Z}}_i\}_{i=1}^n$. Obtain $N(i)$ for $i=1,\dots,n$.
\State Compute
$
\widehat{L}_n^{{\boldsymbol{Z}}} = n^{-1}(n-1)^{-1} \sum_{1 \leq  i \neq j \leq n} \Big\{ \widehat{G}_{{\boldsymbol{Z}}_i}(Y_j) \widehat{G}_{{\boldsymbol{Z}}_{N(i)}}(Y_j) - \widehat{G}^2_{{\boldsymbol{Z}}_i}(Y_j) \Big\}.
$
\Ensure Bias estimator $\widehat{L}_n^{{\boldsymbol{Z}}}$.
\end{algorithmic}
\noindent \textit{Note:} The bias estimator $\widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}$ can be obtained by repeating Algorithm~\ref{alg2} with ${\boldsymbol{Z}}$ replaced by $({\boldsymbol{X}},{\boldsymbol{Z}})$.
\end{algorithm}

Theorem~\ref{thm:bias_correct} below establishes the bias rate and the convergence rate of the bias estimators. For ease of exposition, Assumptions~\ref{assump_A.1}--\ref{assump_A.4}, which are needed in this subsection, are deferred to Appendix~\ref{secA:assump}.

\begin{theorem} \label{thm:bias_correct}
We establish the following results:
\begin{enumerate}[label=(\roman*)]
\item (Bias rate) Assume Assumptions~\ref{assump_4.1}, \ref{assump_4.2}, and \ref{assump_A.1}. Then, as $n\to \infty$,
\begin{eqnarray*}
L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} = O(n^{-2/(p+q)}),  \qquad  L_n^{{\boldsymbol{Z}}} = O(n^{-2/q} + n^{-1}).
\end{eqnarray*}
If $p+q \leq 3$, then $L^{(\tau)}_n  = L_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} - L_n^{{\boldsymbol{Z}}}  =  o(n^{-1/2})$ and $L^{(\kappa)}_n  = - L_n^{{\boldsymbol{Z}}} +O(n^{-1}) =o(n^{-1/2})$, that is, both $L^{(\tau)}_n$ and $L^{(\kappa)}_n$ are asymptotically negligible and bias correction is unnecessary.

\item (Bias-correction efficacy) Assume Assumptions~\ref{assump_4.1}, \ref{assump_4.2}, \ref{assump_A.2}, \ref{assump_A.3}, and \ref{assump_A.4}. Then, as $n \to \infty$, the bias estimators $\widehat{L}_n^{{\boldsymbol{Z}}}$ and $\widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}}$ produced by Algorithm~\ref{alg2} satisfy
\begin{eqnarray*}
\widehat{L}_n^{{\boldsymbol{Z}}} - L_n^{{\boldsymbol{Z}}} = o_{\mathrm P}(n^{-1/2}), \qquad \widehat{L}_n^{{\boldsymbol{X}},{\boldsymbol{Z}}} - L_n^{{\boldsymbol{X}},{\boldsymbol{Z}}} = o_{\mathrm P}(n^{-1/2}).
\end{eqnarray*}
\end{enumerate}
\end{theorem}



\subsection{Inferential validity} \label{sec:infer_theory}

Recall the confidence intervals proposed in \eqref{eq:step3-CI}:
\begin{eqnarray*}
&&\mathsf{CI}_{\alpha}(T) = \Big(T_n -\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n}} \ , \ T_n +\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n} } \Big), \qquad
\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T) = \Big(T_n^{\mathrm{bc}} -\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n}} \ , \ T_n^{\mathrm{bc}} +\frac{z_{\alpha/2}\cdot \widehat{\sigma}}{\kappa_n \cdot\sqrt{n} } \Big).
\end{eqnarray*}
Based on the explicit variance-estimation and bias-correction procedures developed in Section~\ref{sec:theory}, we can now specify the concrete forms of $\widehat{\sigma}$ and $T_n^{\mathrm{bc}}$. Specifically, $\widehat{\sigma}$ is a consistent estimator of $\sigma$, as given in \eqref{eq:sigma2_est} of Theorem~\ref{thm: est_var-main}, whereas
$
T_n^{\mathrm{bc}} =
\big(\tau_n - 	\widehat{L}^{(\tau)}_n\big)/\big(\kappa_n - 	\widehat{L}^{(\kappa)}_n\big),
$
with $\widehat{L}^{(\tau)}_n$ and $\widehat{L}^{(\kappa)}_n$ defined in \eqref{eq:bias_est}.

As a direct consequence of the CLT established in Corollary~\ref{cor: CLT-Tn}, we then obtain the validity of the confidence intervals in Theorem~\ref{thm:CI} below.

\begin{theorem}[Confidence interval validity] \label{thm:CI}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5}. Let $\alpha \in (0,1)$ be a prespecified significance level. Assume that $\sigma^2$ in \eqref{eq:sigma2_tTn} is strictly positive.
\begin{enumerate}[label=(\roman*)]
\item Assume Assumption~\ref{assump_A.1}. If $p+q\leq 3$, then the confidence interval $\mathsf{CI}_{\alpha}(T)$ is valid, in the sense that
\begin{eqnarray*}
\lim_{n \to \infty}{\mathrm P}\big( T \in \mathsf{CI}_{\alpha}(T)\big) = 1-\alpha.
\end{eqnarray*}
\item Assume Assumptions~\ref{assump_A.2}--\ref{assump_A.4}. For general $p,q \geq 1$, the bias-corrected confidence interval $\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T)$ is valid, in the sense that
\begin{eqnarray*}
\lim_{n \to \infty}{\mathrm P}\big( T \in \mathsf{CI}^{\mathrm{bc}}_{\alpha}(T)\big) = 1-\alpha.
\end{eqnarray*}
\end{enumerate}
Moreover, as $n \to \infty$, the lengths of both confidence intervals converge to $0$ at rate $O_{\mathrm P}(n^{-1/2})$.
\end{theorem}

For conditional independence testing, by combining the two limiting variance-estimation methods developed in the latter part of Section~\ref{sec:est_var} with the bias-correction procedure in Section~\ref{sec:bias_correct}, we obtain four distinct level-$\alpha$ tests:
\begin{eqnarray}
&&\mathsf{T}^{\mathrm{F}}_{\alpha} = \mathbf{1} \big(	\sqrt{n} \tau_n / \widehat{\sigma}_{0,\mathrm{F}}> z_\alpha\big), \qquad
\mathsf{T}^{\mathrm{F,bc}}_\alpha = \mathbf{1}\big(	\sqrt{n} (\tau_n - \widehat{L}_n^{(\tau)}) / \widehat{\sigma}_{0,\mathrm{F}} > z_\alpha\big), \cr
&&\mathsf{T}^{\mathrm{B}}_{\alpha} = \mathbf{1}\big(\sqrt{n} \tau_n / \widehat{\sigma}_{0,\mathrm{B}}> z_\alpha\big), \qquad
\mathsf{T}^{\mathrm{B,bc}}_\alpha = \mathbf{1}\big(	\sqrt{n} (\tau_n - \widehat{L}_n^{(\tau)}) / \widehat{\sigma}_{0,\mathrm{B}} > z_\alpha\big),
\label{eq:four_tests}
\end{eqnarray}
where $\widehat{\sigma}_{0,\mathrm{F}}$ and $\widehat{\sigma}_{0,\mathrm{B}}$ correspond to the fast estimator in \eqref{eq:sigma2_F} and the bootstrap estimator in \eqref{eq:sigma2_B}, respectively, and $\widehat{L}_n^{(\tau)} = \widehat{L}_n^{{\boldsymbol{X}}, {\boldsymbol{Z}}} - \widehat{L}_n^{{\boldsymbol{Z}}}$ is the bias estimator.
Let $H_1$ denote the alternative hypothesis consisting of all distributions of $({\boldsymbol{X}}, Y, {\boldsymbol{Z}})$ under which $Y$ is not conditionally independent of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$.
Theorem~\ref{thm:test} establishes the asymptotic size control and consistency of these four tests.

\begin{theorem}[Test validity and consistency] \label{thm:test}
Assume Assumptions~\ref{assump_4.1}--\ref{assump_4.5} and \ref{assump_A.1}--\ref{assump_A.4}. Let $\alpha \in (0,1)$ be a prespecified significance level. Assume $B\to \infty$, $m \to \infty$, and $m = o(n)$.
We have the following results:
\begin{enumerate}[label=(\roman*)]
\item The tests $\mathsf{T}^{\mathrm{F,bc}}_\alpha$ and $\mathsf{T}^{\mathrm{B,bc}}_\alpha$ are valid in the sense that, for any $({\boldsymbol{X}}, Y, {\boldsymbol{Z}})$ whose distribution is fixed and satisfies $H_0$ in \eqref{eq:null},
\begin{eqnarray*}
\lim_{n \to \infty} {\mathrm P}(\mathsf{T}^{\mathrm{F,bc}}_\alpha = 1) = \alpha, \qquad  \lim_{n \to \infty} {\mathrm P}(\mathsf{T}^{\mathrm{B,bc}}_\alpha = 1) = \alpha.
\end{eqnarray*}
Moreover, if $p+q \leq 3$, then the same conclusion also holds for the corresponding non-bias-corrected tests $\mathsf{T}^{\mathrm{F}}_{\alpha}$ and $\mathsf{T}^{\mathrm{B}}_{\alpha}$.
\item The four proposed tests are consistent in the sense that, for any $({\boldsymbol{X}}, Y, {\boldsymbol{Z}})$ whose distribution is fixed and does not satisfy $H_0$ in \eqref{eq:null},
\begin{eqnarray*}
\lim_{n \to \infty} {\mathrm P}(\mathsf{T} = 1) =1,
\qquad \text{for each }  \mathsf{T} \in \{\mathsf{T}^{\mathrm{F}}_{\alpha}, \mathsf{T}^{\mathrm{B}}_{\alpha}, \mathsf{T}^{\mathrm{F,bc}}_\alpha , \mathsf{T}^{\mathrm{B,bc}}_\alpha\}.
\end{eqnarray*}
\end{enumerate}
\end{theorem}

\begin{remark}
Unfortunately, a calculation of the joint limiting distribution of $\tau_n$ and the log-likelihood ratio, analogous to that in \cite{shi2020power} (see also \cite{shi2020rate}), shows that all the tests considered in Theorem~\ref{thm:test} still have Pitman efficiency zero. Nevertheless, by adopting ideas similar to those in \citet[Section~4]{bickel2022measures}, one may either use $T_n$ solely for measuring conditional dependence, or combine the conditional independence test based on $\tau_n$ with any conditional independence test that possesses positive Pitman efficiency. The resulting combined test can be size-adjusted either via a naive Bonferroni correction or through a more refined analysis of the joint distribution of the two test statistics. In this way, one obtains a combined procedure with nonzero local efficiency.
\end{remark}



\section{Simulations} \label{sec:simu}



We conduct simulation studies to investigate the finite-sample performance of the proposed confidence intervals and conditional independence tests.
To this end, we consider the following two data-generating models on $({\boldsymbol{X}}, Y, {\boldsymbol{Z}})$.

\begin{description}
\item[Model 1] (Uniform distribution). Let ${\boldsymbol{Z}} = (Z_1, \dots, Z_q) \sim \mathrm{Unif}[0,1]^q$ follow a $q$-dimensional uniform distribution.
Define the $p$-dimensional extension $\widetilde{\boldsymbol{Z}} = (\widetilde{Z}_1, \dots, \widetilde{Z}_p)$ of ${\boldsymbol{Z}}$, with
\[
\widetilde{Z}_j = Z_j \text{ for } j \leq q,\text{ and } \widetilde{Z}_j = Z_1 \text{ for } j > q \text{ (if $p>q$)} .
\]
Let $\epsilon_0, \epsilon_1, \dots, \epsilon_p \sim_\mathrm{i.i.d.} \mathrm{Unif}[0,1]$ be independent noise variables.
For a prespecified $\rho \in [0,1]$, define $Y$ and ${\boldsymbol{X}} = (X_1,\dots, X_p)$ as
\begin{eqnarray}
Y = Z_q + \rho \,\epsilon_1 + \sqrt{1-\rho^2} \, \epsilon_0, \qquad \text{and } \quad  X_j = \widetilde{Z}_j + \epsilon_j, \quad \text{for } j = 1\dots, p.  \label{eq:model1}
\end{eqnarray}
It is clear that as $\rho$ increases, the conditional dependence of $Y$ on ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$ becomes stronger. In particular, $\rho = 0$ corresponds to the null hypothesis \eqref{eq:null}, whereas $\rho = 1$ corresponds to the case where $Y$ is a deterministic function of ${\boldsymbol{X}}$ given ${\boldsymbol{Z}}$.

\item[Model 2] (Gaussian distribution). Model 2 is structurally similar to Model 1, except that all $\mathrm{Unif}[0,1]$ distributions are replaced with the standard normal $N(0,1)$ distributions.
\end{description}

For the settings of dimensions $p, q$, we consider the following five scenarios
\[
(p,q) \in \{(1,1), (1,2), (3,1), (3,3), (5,5)\}.
\]
In each simulation run, we generate $\mathrm{i.i.d.}$  $\{({\boldsymbol{X}}_i, Y_i, {\boldsymbol{Z}}_i)\}_{i=1}^n$ from the selected model, with sample size
\[
n \in \{1000,5000,10000\}.
\]
For the the confidence interval $\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T) $ and testing
 methods $\mathsf{T}^{\mathrm{F,bc}}$, $\mathsf{T}^{\mathrm{B,bc}}$ that involve bias correction, we set the penalty parameter $\lambda_n = \lceil n^{-0.85} \rceil$ in Algorithm~\ref{alg2}.
 For the
methods $\mathsf{T}^{\mathrm{B}}$ and $\mathsf{T}^{\mathrm{B,bc}}$ that employ the $m$-out-of-$n$ bootstrap, we set the number of bootstrap replications to $B = 200$ and the subsample size to $m = \lceil n^{0.5} \rceil$. In each scenario, the true values of $T$ and $\sigma^2$ (or $\sigma_0^2$) are approximated by Monte Carlo simulation by evaluating $T_n$ and $\widehat{\sigma}^2$ with $n=10^7$.

Note that {\bf Model 2} represents a challenging setting in which the observed random variables are unbounded, making bias correction---which essentially amounts to a nonparametric regression adjustment---substantially more difficult. We include this setting to contrast it with {\bf Model 1}, which represents the most idealized case, and to assess the robustness of our methods in scenarios that fall outside the scope of the available theoretical guarantees.

\vspace{0.3cm}

All code required to reproduce the results in this paper is publicly available at: \\\url{https://github.com/MuhongGao/Conditional_Independence}.


\subsection{Empirical coverage probabilities}
Tables \ref{Table_2} and \ref{Table_3} report the empirical coverage probabilities (ECPs) of the proposed confidence intervals
$\mathsf{CI}_{\alpha}(T) $ and $\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T) $ under Models 1 and 2, respectively. In addition, the tables report the relative empirical root mean squared error (rRMSE), defined below, to assess the accuracy of the limiting variance estimator $\widehat{\sigma}^2$ in \eqref{eq:sigma2_est}:
\begin{eqnarray}
	\text{rRMSE} = \sqrt{{\mathrm E} \Big\{\Big(\frac{\widehat{\sigma}^2 -\sigma^2}{\sigma^2}\Big)^2\Big\}}.   \label{eq:rRMSE}
\end{eqnarray}
From the tables, we observe that when $(p,q)=(1,1)$ and $(1,2)$, so that $p+q\leq 3$, both $\mathsf{CI}$ and $\mathsf{CI}^{\mathrm{bc}}$ achieve ECPs close to the nominal level $1-\alpha=0.9$. By contrast, when $(p,q)=(3,3)$, so that $p+q>3$, the ECPs of $\mathsf{CI}$ deteriorate substantially, whereas those of $\mathsf{CI}^{\mathrm{bc}}$ remain close to $0.9$. This pattern is especially pronounced under {\bf Model 1} when $\rho$ is not too close to $1$. These findings are consistent with Theorem~\ref{thm:bias_correct}, which indicates that bias correction is necessary when $p+q>3$, and they further demonstrate the effectiveness of the proposed bias-correction procedure. The case $(p,q)=(5,5)$, on the other hand, illustrates the curse of dimensionality, as one would expect.

The situation is also of interest under {\bf Model 2}, where the assumptions required for bias correction are violated because $({\boldsymbol{X}},{\boldsymbol{Z}})$ is supported on an unbounded domain. In this case, as shown in Table~\ref{Table_3}, both $\mathsf{CI}$ and $\mathsf{CI}^{\mathrm{bc}}$ continue to perform well when $p+q\leq 3$, suggesting that bias correction may serve as a safe alternative to $\mathsf{CI}$, albeit at a higher computational cost. On the other hand, when $p+q>3$, $\mathsf{CI}^{\mathrm{bc}}$ no longer performs as well as it does in Table~\ref{Table_2}, although it remains substantially superior to $\mathsf{CI}$. This suggests that the bias-correction step is indeed sensitive to tail observations, and points to the potential value of applying a rank transformation, as in \citet[Section~4]{cattaneo2025rosenbaum}, to stabilize the bias correction. Given the already broad scope of the present paper, we do not pursue this direction further.

\begin{table}[tbp]
	\caption{
		\textbf{(Model 1: ECPs).} $\mathsf{CI}$ and $\mathsf{CI}^{\mathrm{bc}}$ denote the ECPs of the proposed $(1-\alpha)$-level confidence intervals $\mathsf{CI}_{\alpha}(T)$ and $\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T)$, respectively, with $\alpha = 0.1$. rRMSE represents the relative empirical root mean squared error of $\widehat{\sigma}^2$.
		The results are based on 1000 simulation replicates.
	}
	\label{Table_2}
	\resizebox{\linewidth}{!}{
		\begin{tabular}{cr@{\hspace{6mm}}cccccccccccc}
			\toprule
			\multicolumn{2}{c}{} & \multicolumn{3}{c}{$(p,q)=(1,1)$} & \multicolumn{3}{c}{$(p,q)=(1,2)$} & \multicolumn{3}{c}{$(p,q)=(3,3)$} & \multicolumn{3}{c}{$(p,q)=(5,5)$} \\[1mm]
			\cmidrule(lr){3-5}  \cmidrule(lr){6-8}
			\cmidrule(lr){9-11}  \cmidrule(lr){12-14}
			$\rho$ & \multicolumn{1}{c}{$n$} & $\mathsf{CI}$     & $\mathsf{CI}^{\mathrm{bc}}$  & \footnotesize rRMSE  & $\mathsf{CI}$    & $\mathsf{CI}^{\mathrm{bc}}$ & \footnotesize rRMSE  & $\mathsf{CI}$    & $\mathsf{CI}^{\mathrm{bc}}$ & \footnotesize rRMSE  & $\mathsf{CI}$    & $\mathsf{CI}^{\mathrm{bc}}$ & \footnotesize rRMSE \\[1mm] \midrule
			\multirow{3}[0]{*}{0} & 1000  & 0.88  & 0.88  & 0.40  & 0.85  & 0.86  & 0.39  & 0.81  & 0.85  & 0.40  & 0.72  & 0.85  & 0.39 \\      & 5000  & 0.88  & 0.88  & 0.18  & 0.90  & 0.90  & 0.18  & 0.79  & 0.90  & 0.17  & 0.44  & 0.89  & 0.18 \\      & 10000 & 0.89  & 0.89  & 0.13  & 0.89  & 0.89  & 0.13  & 0.78  & 0.88  & 0.12  & 0.27  & 0.89  & 0.12 \\[2mm]
			\multirow{3}[0]{*}{0.3} & 1000  & 0.89  & 0.88  & 0.40  & 0.88  & 0.88  & 0.38  & 0.78  & 0.86  & 0.40  & 0.70  & 0.84  & 0.39 \\      & 5000  & 0.89  & 0.89  & 0.18  & 0.91  & 0.91  & 0.18  & 0.74  & 0.89  & 0.17  & 0.44  & 0.82  & 0.18 \\      & 10000 & 0.90  & 0.90  & 0.12  & 0.89  & 0.88  & 0.13  & 0.74  & 0.88  & 0.12  & 0.30  & 0.78  & 0.12 \\[2mm]
			\multirow{3}[0]{*}{0.5} & 1000  & 0.88  & 0.88  & 0.39  & 0.88  & 0.88  & 0.37  & 0.70  & 0.87  & 0.40  & 0.48  & 0.86  & 0.41 \\      & 5000  & 0.88  & 0.88  & 0.18  & 0.89  & 0.90  & 0.17  & 0.64  & 0.90  & 0.17  & 0.15  & 0.81  & 0.20 \\      & 10000 & 0.91  & 0.91  & 0.12  & 0.89  & 0.89  & 0.12  & 0.59  & 0.88  & 0.12  & 0.06  & 0.72  & 0.14 \\[2mm]
			\multirow{3}[0]{*}{0.7} & 1000  & 0.89  & 0.89  & 0.40  & 0.87  & 0.88  & 0.38  & 0.48  & 0.89  & 0.52  & 0.08  & 0.91  & 0.52 \\      & 5000  & 0.90  & 0.90  & 0.18  & 0.88  & 0.90  & 0.18  & 0.28  & 0.93  & 0.26  & 0.00  & 0.87  & 0.29 \\      & 10000 & 0.89  & 0.89  & 0.13  & 0.89  & 0.89  & 0.12  & 0.19  & 0.91  & 0.20  & 0.00  & 0.71  & 0.23 \\[2mm]
			\multirow{3}[0]{*}{0.9} & 1000  & 0.85  & 0.85  & 0.54  & 0.84  & 0.92  & 0.89  & 0.01  & 0.76  & 1.93  & 0.00  & 0.72  & 0.97 \\      & 5000  & 0.89  & 0.89  & 0.26  & 0.84  & 0.91  & 0.37  & 0.00  & 0.72  & 1.35  & 0.00  & 0.80  & 0.78 \\      & 10000 & 0.88  & 0.88  & 0.18  & 0.86  & 0.92  & 0.24  & 0.00  & 0.69  & 1.11  & 0.00  & 0.92  & 0.70 \\
			\bottomrule
		\end{tabular}
	}
\end{table}


Regarding the rRMSE, it decreases uniformly with $n$ across all scenarios, providing empirical evidence for the consistency of $\widehat{\sigma}^2$ established in Theorem \ref{thm: est_var-main}. Moreover, the rRMSE increases noticeably with both $p+q$ and $\rho$, indicating that the convergence rate of $\widehat{\sigma}^2$ deteriorates in higher-dimensional settings and under stronger dependence. Such behavior is in line with theoretical intuition.



\subsection{Empirical powers of tests of $H_0$} \label{sec:powers}

The empirical powers under {\bf Models 1} and {\bf 2} are reported in Figures~\ref{Figure_1} and \ref{Figure_2}, respectively.
We begin with the results under {\bf Model 1} (Figure~\ref{Figure_1}). When $(p,q) = (1,1), (1,2)$, or $(3,1)$, all four methods perform similarly well, with power curves that nearly overlap. By contrast, when $(p,q) = (3,3)$ or $(5,5)$, the bias-corrected methods $\mathsf{T}^{\mathrm{F,bc}}$ and $\mathsf{T}^{\mathrm{B,bc}}$ clearly outperform their non-bias-corrected counterparts $\mathsf{T}^{\mathrm{F}}$ and $\mathsf{T}^{\mathrm{B}}$, and this advantage becomes more pronounced as $(p,q)$ increases. This pattern is in line with Theorem~\ref{thm:bias_correct}, which shows that bias correction is unnecessary when $p+q < 4$, whereas its benefit becomes increasingly substantial as the dimension grows.
Furthermore, when $p$ and $q$ are large, although the non-bias-corrected methods $\mathsf{T}^{\mathrm{F}}$ and $\mathsf{T}^{\mathrm{B}}$ control size under $H_0$ (that is, when $\rho = 0$) well below the nominal level $\alpha = 0.05$, their power increases much more slowly as $\rho$ grows. This suggests that, under {\bf Model 1}, the bias $L_n = {\mathrm E}(\tau_n) - \tau$ is positive, rendering the tests more conservative and thereby lowering their rejection probabilities.
Finally, comparing the two variance-estimation methods, we observe little difference in either size or empirical power: $\mathsf{T}^{\mathrm{F,bc}}$ and $\mathsf{T}^{\mathrm{B,bc}}$ behave almost identically, and the same is true for $\mathsf{T}^{\mathrm{F}}$ and $\mathsf{T}^{\mathrm{B}}$.

\begin{table}[tbp]
	\caption{
		\textbf{(Model 2: ECPs).} $\mathsf{CI}$ and $\mathsf{CI}^{\mathrm{bc}}$ denote the ECPs of the proposed $(1-\alpha)$-level confidence intervals $\mathsf{CI}_{\alpha}(T)$ and $\mathsf{CI}^{\mathrm{bc}}_{\alpha}(T)$, respectively, with $\alpha = 0.1$. rRMSE represents the relative empirical root mean squared error of $\widehat{\sigma}^2$. The results are based on 1000 simulation replicates.
	}
	\label{Table_3}
	\resizebox{\linewidth}{!}{
		\begin{tabular}{cr@{\hspace{6mm}}cccccccccccc}\toprule\multicolumn{2}{c}{} & \multicolumn{3}{c}{$(p,q)=(1,1)$} & \multicolumn{3}{c}{$(p,q)=(1,2)$} & \multicolumn{3}{c}{$(p,q)=(3,3)$} & \multicolumn{3}{c}{$(p,q)=(5,5)$} \\[1mm]
			\cmidrule(lr){3-5}  \cmidrule(lr){6-8}
			\cmidrule(lr){9-11}  \cmidrule(lr){12-14}
			$\rho$ & \multicolumn{1}{c}{$n$} & \multicolumn{1}{c}{$\mathsf{CI}$} & \multicolumn{1}{c}{$\mathsf{CI}^{\mathrm{bc}}$ } & \multicolumn{1}{c}{\footnotesize rRMSE} & \multicolumn{1}{c}{$\mathsf{CI}$} & \multicolumn{1}{c}{$\mathsf{CI}^{\mathrm{bc}}$ } & \multicolumn{1}{c}{\footnotesize rRMSE} & \multicolumn{1}{c}{$\mathsf{CI}$} & \multicolumn{1}{c}{$\mathsf{CI}^{\mathrm{bc}}$ } & \multicolumn{1}{c}{\footnotesize rRMSE} & \multicolumn{1}{c}{$\mathsf{CI}$} & \multicolumn{1}{c}{$\mathsf{CI}^{\mathrm{bc}}$ } & \multicolumn{1}{c}{\footnotesize rRMSE} \\[1mm]
			\midrule\multirow{3}[1]{*}{0} & 1000  & 0.84  & 0.84  & 0.42  & 0.86  & 0.84  & 0.41  & 0.75  & 0.83  & 0.42  & 0.67  & 0.80  & 0.44 \\      & 5000  & 0.90  & 0.89  & 0.19  & 0.89  & 0.88  & 0.19  & 0.64  & 0.86  & 0.18  & 0.28  & 0.87  & 0.19 \\      & 10000 & 0.90  & 0.90  & 0.13  & 0.90  & 0.90  & 0.13  & 0.61  & 0.88  & 0.13  & 0.12  & 0.87  & 0.13 \\[2mm]
			\multirow{3}[0]{*}{0.3} & 1000  & 0.84  & 0.83  & 0.41  & 0.86  & 0.83  & 0.40  & 0.71  & 0.83  & 0.42  & 0.65  & 0.71  & 0.44 \\      & 5000  & 0.91  & 0.90  & 0.18  & 0.89  & 0.88  & 0.18  & 0.59  & 0.85  & 0.18  & 0.36  & 0.67  & 0.19 \\      & 10000 & 0.91  & 0.91  & 0.13  & 0.91  & 0.91  & 0.12  & 0.53  & 0.86  & 0.12  & 0.19  & 0.55  & 0.13 \\[2mm]
			\multirow{3}[0]{*}{0.5} & 1000  & 0.85  & 0.85  & 0.39  & 0.85  & 0.86  & 0.38  & 0.60  & 0.83  & 0.42  & 0.43  & 0.70  & 0.45 \\      & 5000  & 0.91  & 0.91  & 0.17  & 0.89  & 0.89  & 0.17  & 0.40  & 0.83  & 0.18  & 0.09  & 0.57  & 0.20 \\      & 10000 & 0.91  & 0.91  & 0.12  & 0.92  & 0.90  & 0.12  & 0.29  & 0.79  & 0.12  & 0.02  & 0.41  & 0.14 \\[2mm]
			\multirow{3}[0]{*}{0.7} & 1000  & 0.84  & 0.84  & 0.38  & 0.84  & 0.87  & 0.39  & 0.30  & 0.88  & 0.51  & 0.06  & 0.72  & 0.52 \\      & 5000  & 0.90  & 0.90  & 0.17  & 0.88  & 0.88  & 0.17  & 0.07  & 0.84  & 0.27  & 0.00  & 0.51  & 0.26 \\      & 10000 & 0.90  & 0.90  & 0.12  & 0.91  & 0.89  & 0.12  & 0.02  & 0.76  & 0.21  & 0.00  & 0.25  & 0.21 \\[2mm]
			\multirow{3}[1]{*}{0.9} & 1000  & 0.84  & 0.84  & 0.49  & 0.78  & 0.95  & 1.03  & 0.00  & 0.94  & 1.44  & 0.00  & 0.97  & 0.68 \\      & 5000  & 0.90  & 0.90  & 0.23  & 0.78  & 0.86  & 0.46  & 0.00  & 0.97  & 1.10  & 0.00  & 0.95  & 0.45 \\      & 10000 & 0.90  & 0.89  & 0.16  & 0.81  & 0.82  & 0.32  & 0.00  & 0.98  & 0.96  & 0.00  & 0.76  & 0.42 \\\bottomrule\end{tabular}
	}
\end{table}

We next examine the results under {\bf Model 2} (Figure~\ref{Figure_2}). Overall, the patterns are highly consistent with those observed under {\bf Model 1}, exhibiting similar trends and leading to the same qualitative conclusions. A closer inspection reveals that the discrepancies are, if anything, slightly more pronounced under {\bf Model 2}. In particular, (i) when $(p,q) = (1,2)$ or $(3,1)$, the four methods exhibit less overlap in their power curves; and (ii) when $(p,q) = (5,5)$ and the sample size is relatively small (e.g., $n = 1000$), the bias-corrected methods $\mathsf{T}^{\mathrm{F,bc}}$ and $\mathsf{T}^{\mathrm{B,bc}}$ exhibit size inflation under $H_0$ (that is, when $\rho = 0$), with rejection probabilities around $0.2$, substantially above the nominal level $0.05$. However, this issue is quickly alleviated as the sample size increases.

Overall, these results show that, unlike the confidence-interval counterpart, the proposed testing procedures exhibit consistently stable and favorable performance across different models, even in the presence of unbounded distributions. This is in line with the general intuition that testing is often statistically easier than estimation.

\subsection{Comparison of variance estimators under $H_0$}

We next take a closer look at the performance of the two asymptotic variance estimators under $H_0$: the fast $k$NN-based estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$ in \eqref{eq:sigma2_F}$,$ and the $m$-out-of-$n$ bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$ in \eqref{eq:sigma2_B}.

We focus on the following three aspects:
\begin{enumerate}[label=(\roman*)]
\item the rRMSE for estimating $\sigma_0^2$ (defined analogously to \eqref{eq:rRMSE});
\item the ECP of $\tau$, defined as the empirical probability that the $(1-\alpha)$-level confidence interval
\[
\Big( \tau_n -z_{1-\alpha/2} \widehat{\sigma} /\sqrt{n} \ , \  \tau_n +z_{1-\alpha/2} \widehat{\sigma} /\sqrt{n}\Big)
\]
covers the true value $\tau = 0$ under $H_0$, where $\widehat{\sigma}$ denotes either $\widehat{\sigma}_{0,\mathrm{F}}$ or $\widehat{\sigma}_{0,\mathrm{B}}$, and where we set $\alpha = 0.05$;
\item the CPU time required to compute each estimator. All numerical experiments were implemented in MATLAB R2023b on a Windows desktop equipped with an Intel Xeon Platinum 8370C 64-core processor.
\end{enumerate}

Figures~\ref{Figure_4} and \ref{Figure_5} report the results for {\bf Models 1} and {\bf 2}, respectively. The overall patterns in the two figures are nearly identical. We therefore focus on Figure~\ref{Figure_4} under {\bf Model 1}.
For the rRMSE, the bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$ consistently maintains a relatively low level, whereas the fast estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$ exhibits larger values when the sample size is small (e.g., $n = 1000$). As $n$ increases, however, the rRMSE of $\widehat{\sigma}^2_{0,\mathrm{F}}$ decreases substantially, eventually becoming comparable to, and even slightly smaller than, that of $\widehat{\sigma}^2_{0,\mathrm{B}}$ when $n = 10000$. This suggests that $\widehat{\sigma}^2_{0,\mathrm{F}}$ converges more slowly than $\widehat{\sigma}^2_{0,\mathrm{B}}$, but attains comparable performance once the sample size is sufficiently large.

Despite these differences in rRMSE, they do not translate into noticeable differences in statistical inference. Both in terms of empirical power (see Section~\ref{sec:powers}) and ECP, the two estimators yield nearly identical results, with their corresponding curves largely overlapping.

Finally, in terms of computational efficiency, the $k$NN-based estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$ is substantially faster to compute than the bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$. Even with a relatively small number of bootstrap replications ($B = 200$), the $m$-out-of-$n$ bootstrap incurs a much higher computational cost than the $k$NN-based estimator.

In summary, $\widehat{\sigma}^2_{0,\mathrm{B}}$ achieves lower rRMSE when the sample size is small, whereas $\widehat{\sigma}^2_{0,\mathrm{F}}$ becomes increasingly competitive, faster, and more stable as the sample size grows.


\begin{figure}[!htbp]
\centering
\raisebox{1 em}{
\includegraphics[width= \linewidth]
{figures/Model1_power.eps}}
\caption{\textbf{(Model 1: power comparison)}
Empirical powers of the four proposed tests in \eqref{eq:four_tests} under Model 1, at nominal significance level $\alpha = 0.05$. The $y$-axis shows rejection frequencies based on $1000$ replicates, while the $x$-axis represents the model parameter $\rho$ in \eqref{eq:model1}, which characterizes the strength of conditional correlation. The case $\rho = 0$ corresponds to the null hypothesis.
}
\label{Figure_1}
\end{figure}

\begin{figure}[!htbp]
\centering
\raisebox{1 em}{
\includegraphics[width= \linewidth]
{figures/Model2_power.eps}}
\caption{\textbf{(Model 2: power comparison)}
Empirical powers of the four proposed tests in \eqref{eq:four_tests} under Model 3, at  nominal  significance level $\alpha = 0.05$. The $y$-axis shows rejection frequencies based on $1000$ replicates, while the $x$-axis represents the model parameter $\rho$ in \eqref{eq:model1}, which characterizes the strength of conditional correlation. The case $\rho = 0$ corresponds to the null hypothesis.
}
\label{Figure_2}
\end{figure}


\begin{figure}[!htbp]
\centering
\raisebox{1 em}{
\includegraphics[width= \linewidth]
{figures/Model1_sigma.eps}}
\caption{
\textbf{(Model 1: comparison of variance estimators under $H_0$)}
For the fast $k$NN-based estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$ in \eqref{eq:sigma2_F} and the $m$-out-of-$n$ bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$ in \eqref{eq:sigma2_B}, the first row reports the rRMSE values, the second row reports the ECP values, and the third row presents the CPU time (in seconds). All values are averaged over $1000$ simulation replicates.
}
\label{Figure_4}
\end{figure}

\begin{figure}[!htbp]
\centering
\raisebox{1 em}{
\includegraphics[width= \linewidth]
{figures/Model2_sigma.eps}}
\caption{
\textbf{(Model 2: comparison of variance estimators under $H_0$)}
For the fast $k$NN-based estimator $\widehat{\sigma}^2_{0,\mathrm{F}}$ in \eqref{eq:sigma2_F} and the $m$-out-of-$n$ bootstrap estimator $\widehat{\sigma}^2_{0,\mathrm{B}}$ in \eqref{eq:sigma2_B}, the first row reports the rRMSE values, the second row reports the ECP values, and the third row presents the CPU time (in seconds). All values are averaged over $1000$ simulation replicates.
}
\label{Figure_5}
\end{figure}



\clearpage

{\centering\huge Supplement to ``Limit theorems of Azadkia-Chatterjee's conditional graph correlation'' \par}