# Distribution-Free Joint Independence Testing and Robust Independent
Component Analysis Using Optimal Transport

Ziang Niu  Applied Mathematics and Computational Science, University of
Pennsylvania, Philadelphia, USA<EMAIL_ADDRESS>and  Bhaswar B.
Bhattacharya  Department of Statistics and Data Science, University of
Pennsylvania, Philadelphia, USA<EMAIL_ADDRESS>
###### Abstract.

In this paper, we study the problem of measuring and testing joint
independence for a collection of multivariate random variables. Using the
emerging theory of optimal transport (OT)-based multivariate ranks, we propose
a distribution-free test for multivariate joint independence. Towards this we
introduce the notion of rank joint distance covariance $(\mathrm{RJdCov})$,
the higher-order rank analogue of the celebrated distance covariance measure,
which can capture the dependencies among all the subsets of the variables. The
$\mathrm{RJdCov}$ can be easily estimated from the data without any moment
assumptions and the associated test for joint independence is universally
consistent. We derive the asymptotic null distribution of the
$\mathrm{RJdCov}$ estimate, using which we can readily calibrate the test
without any knowledge of the (unknown) marginal distributions (due to the
distribution-free property). We also provide an efficient data-agnostic
resampling-based implementation of the test which controls Type-I error in
finite samples and is consistent with only a fixed number of resamples. In
addition to being distribution-free and universally consistent, the proposed
test is also statistically efficient, that is, it has non-trivial asymptotic
(Pitman) efficiency (power against $1/\sqrt{n}$ alternatives). We demonstrate
this by computing the limiting local power of the test for both mixture
alternatives and joint Konijn alternatives. We then use the $\mathrm{RJdCov}$
measure to develop a new method for independent component analysis (ICA) that
is easy to implement and robust to outliers and contamination. Extensive
simulations are performed to illustrate the efficacy of the proposed test in
comparison to other existing methods. Finally, we apply the proposed method to
learn the higher-order dependence structure among different US industries
based on stock prices. As a byproduct of our theoretical analysis, we develop
a version of Hoeffding’s classical combinatorial central theorem for mutiple
independent permutations and a multivariate Hájek representation result for
joint rank statistics, which might be of independent interest.

###### Key words and phrases:

Asymptotic efficiency, combinatorial central limit theorem, distance
correlation, joint dependence measures, multivariate ranks, optimal transport,
Stein’s method.

##  1\. Introduction

Given $r$ probability distributions on
$\mathbb{R}^{d_{1}},\mathbb{R}^{d_{2}},\ldots,\mathbb{R}^{d_{r}}$ with
marginal laws $\mu_{1},\mu_{2},\ldots,\mu_{r}$, respectively, and joint law
$\mu$ on $\mathbb{R}^{d}$, where $d=\sum_{s=1}^{r}d_{s}$, the mutual
independence testing problem is given by

$\displaystyle
H_{0}:\mu=\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}\qquad\text{versus}\qquad
H_{1}:\mu\neq\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r},$ (1.1)

where $\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}$ is the product of the
marginal distributions $\mu_{1},\mu_{2},\ldots,\mu_{r}$. This has been
extensively studied in the case $r=2$, which is the classical pairwise
independence testing problem for a collection of random variables. For
univariate distributions, that is, $d_{1}=d_{2}=1$, nonparametric tests for
pairwise independence begin with the celebrated results of Hoeffding [37] and
Blum et al. [7]. These tests are distribution-free, that is, the distributions
of the test statistics under $H_{0}$ do not depend on the (unknown) marginal
distributions, and consistent against a general class of alternatives. Since
then a slew of methods for measuring and testing univariate pairwise
independence have been proposed. These include, among several others, results
of Yanagimoto [75], Feuerverger [21], Bergsma and Dassios [5], Nandy et al.
[56], Heller et al. [36] and the recent breakthrough of Chatterjee [14, 15].

Pairwise independence testing is more challenging for dimensions greater than
1, that is, either $d_{1}>1$ or $d_{2}>1$, because of the lack of a natural
ordering in multivariate data. One of the most popular methods for pairwise
multivariate independence testing is the celebrated distance covariance
($\mathrm{dCov}$) of Székely et al. [70] and Székely and Rizzo [69], which
measures dependence between 2 random vectors based on the difference of their
joint and marginal characteristic functions and can be estimated using the
correlation among the pairwise distances. The distance covariance is zero if
and only if the null hypothesis of independence holds, provided $\mu_{1}$ and
$\mu_{2}$ have finite first moments. Another approach is based on the Hilbert-
Schmidt independence criterion (HSIC), which also provides a consistent
kernel-based test for pairwise independence (see Gretton et al. [28, 27, 29]).
Sejdinovic et al. [65] showed that distance covariance and the kernel-based
independence criterion are, in fact, equivalent if the kernel is chosen based
on the relevant distance function. Other popular methods include mutual
information-based tests [2, 6, 43, 74], graph-based methods [19, 34, 22], the
maximal information coefficient [61], ranking of interpoint distances [35,
55], ball covariance [58], and binning approaches based on partitions of the
sample space [35, 50, 77] (see [41, 49] for reviews of the various methods).
However, none of the aforementioned multivariate tests simultaneously inherit
the distribution-free and universal consistency properties of the rank-based
univariate tests. A breakthrough in this direction was made recently by Deb
and Sen [18] and Shi et al. [67] using optimal transport (OT) based
multivariate ranks [16, 31, 30]. The OT-based tests for pairwise independence
are distribution-free in finite samples, computationally feasible, universally
consistent (that is, the power of the tests converge to 1 as the sample size
increases), and enjoy attractive efficiency properties [68] (see also [20,
50]).

In this paper we consider the mutual independence testing problem (1.1) for
more than 2 variables. Measuring and testing higher-order dependence111We say
a collection of random variables $X_{1},\ldots,X_{r}$ has higher-order
dependency if they are pairwise independent but not jointly independent. and
understanding the dependence structure between the variables are much more
involved tasks than detecting pairwise dependence. Examples of such higher-
order dependencies abound in the literature. This includes, among others,
applications in diagnostic checking for directed acyclic graphs (DAGs), where
the noise variables are assumed to be jointly independent, hence inferring
pairwise independence is not enough, and independent component analysis (ICA),
which entails finding a suitable transformation of multivariate data with
mutually independent components (see Section 5 for details). In fact, Böttcher
[8] collects over 350 datasets featuring statistically significant higher-
order dependencies. While there are a few results for mutual independence
testing in the univariate case, where $d_{1}=d_{2}=\cdots=d_{r}=1$ (see [4,
17, 23, 42, 60] and the references therein), in higher dimensions the problem
is much more challenging. Towards this, Pfister et al. [59] generalized the
HSIC to multiple variables and obtained kernel-based tests joint independence
(see also Sejdinovic et al. [64] for a kernel-based 3-variable interaction
test). Recently, Böttcher et al. [9] and Chakraborty and Zhang [12] proposed
higher-order generalizations of the distance covariance, which provide
consistent tests for multivariate joint independence. This is referred to as
the total distance multivariance by Böttcher et al. [9] or the joint distance
covariance ($\mathrm{JdCov}$) by Chakraborty and Zhang [12]. One of the
attractive properties of the distance multivariance/$\mathrm{JdCov}$ is that
it has a hierarchical structure that can capture the dependence among any
subset of the variables. Recently, Roy et al. [63] proposed tests for joint
independence based on pairwise distances and linear projections (see also
Banerjee and Ghosh [3] for generalizations to functional data). However, none
of the aforementioned methods possesses the distribution-free property of the
univariate methods. In fact, these tests usually have intractable (non-
Gaussian) asymptotic distributions under $H_{0}$ that depend on the (unknown)
marginal distributions and are generally calibrated using permutation methods.

In this paper we propose a distribution-free test for joint independence using
the framework of OT-based multivariate ranks. Towards this, we introduce the
notion of rank joint distance covariance $(\mathrm{RJdCov})$, the rank
analogue of the $\mathrm{JdCov}$, which is obtained by aggregating the
(higher-order) rank distance covariances ($\mathrm{RdCov}$) over all the
subsets of variables (Section 2.4). The higher-order $\mathrm{RdCov}$
characterizes the mutual independence of a subset of variables given their
lower-order independence and the $\mathrm{RJdCov}$ characterizes the total
mutual independence of all the $r$ variables (Proposition 2.10). The
$\mathrm{RdCov}$ and, hence, the $\mathrm{RJdCov}$ can be consistently
estimated from data without any moment assumptions (Theorem 2.12).
Consequently, since OT-based multivariate ranks are themselves distribution-
free (see Section 2.3), we can construct a distribution-free test for
multivariate joint independence based on the estimated $\mathrm{RJdCov}$. The
proposed test has the following properties:

  * •

Distribution-free, universally consistent tests: Since $\mathrm{RJdCov}$ is
well-defined whenever the marginal distributions
$\mu_{1},\mu_{2},\ldots,\mu_{r}$ are absolutely continuous and is zero if and
only if the null hypothesis of joint independence holds, our proposed test is
consistent whenever the joint distribution has absolutely continuous marginals
(Proposition 3.2). Moreover, we can use the $\mathrm{RdCov}$ estimates to
obtain distribution-free, universally consistent tests for higher-order
independence of a subset of variables given their lower-independence. For
example, if a collection variables are pairwise independent, the
$\mathrm{RJdCov}$ measure can be used to construct distribution-free,
consistent tests for 3-way or higher-order dependence (see Remark 3.3 for a
discussion and Section 7 for an application in real data).

  * •

Asymptotic null distribution: In Section 3.1 we derive the limiting
distribution of the proposed test under the null hypothesis of joint
independence. In fact, we derive the limiting joint distribution of the
$\mathrm{RdCov}$ estimates for all the subsets of $r$ variables, which can be
expressed as a squared integral of a certain Gaussian process or,
equivalently, as an infinite weighted sum of independent chi-squared
distributions (Theorem 3.7). One remarkable property of the limiting
distribution is that the distributions of the $\mathrm{RdCov}$ estimates
across the various subsets are asymptotically independent, thereby
facilitating the hierarchical testing for higher-order dependence mentioned
above. Since the OT-based multivariate ranks are distributed as uniform random
permutations, the $\mathrm{RdCov}$ estimates can be expressed as a
combinatorial sum indexed by multiple (more than 1) independent permutations.
For this we use Stein’s method based on exchangeable pairs to develop a
version of Hoeffiding’s classical combinatorial central limit theorem for
multi-dimensional arrays (tensors) indexed by independent permutations (see
Theorem A.1 in Appendix A), a result which might be of independent interest.
The distribution-free property implies that the asymptotic null distribution
of the test statistic does not depend on the distribution of the data
generating mechanism, hence, can be readily used to calibrate the test
statistic.

  * •

Finite sample properties: The distribution-free property also allows us to
approximate the quantiles of the null distribution in finite samples using a
data-agnostic resampling technique. Consequently, we can calibrate the test
statistic to have finite sample level $\alpha$ (Section 3.2). Interestingly,
this finite sample implementation is consistent with only a finite number of
resamples (Theorem 3.11). Hence, we can implement the proposed test both
asymptotically and in finite samples efficiently, without having to estimate
any nuisance parameters of the marginal distributions or use computationally
intensive data-dependent permutation techniques to approximate the rejection
thresholds.

  * •

Asymptotic efficiency: In Section 4 we establish the asymptotic (Pitman)
efficiency of the proposed test by computing its limiting power for local
contiguous alternatives (alternatives shrinking towards $H_{0}$ at rate
$O(1/\sqrt{n})$) [71]. Specifically, we derive the asymptotic distribution of
the proposed test for two kinds of local alternatives: mixture alternatives
(Theorem 4.1) and joint Konijn alternatives (Theorem 4.4). To the best of our
knowledge, these are the first results on the efficiency properties of
nonparametric joint independence tests. This implies that the proposed test,
in addition to being distribution-free, universally consistent, and
computationally feasible, is also statistically efficient, making it
particularly attractive for modern data applications. The proof of the local
power analysis is based on a Hajék representation result for the
$\mathrm{RdCov}$ estimates (Proposition 4.5), which allows us to replace the
empirical rank maps with their population counterparts without altering the
limiting distribution of the test statistic under $H_{0}$.

Going beyond hypothesis testing, the $\mathrm{RJdCov}$ can be used more
generally to quantify deviations from joint independence. Specifically, in
Section 5 using the $\mathrm{RJdCov}$ measure as an objective function we
develop a new method for the independent component analysis (ICA) problem that
is robust to outliers and contamination. ICA is a method for extracting
independent components from multivariate data that emerged from research in
artificial neural networks and has found applications in blind source
separation, feature extraction, computational biology, and time series
analysis (see [39] for a review). Our estimator, which is obtained by
minimizing the $\mathrm{RJdCov}$ measure, can be computed efficiently and is
consistent for the independent components (Theorem 5.2). We illustrate the
effectiveness of the proposed ICA estimator with the approach of Matteson and
Tsay [52] based on the $\mathrm{dCov}$ in various simulation settings (Section
6.4).

In Section 6 we compare the performance of the proposed test with other
popular tests for joint independence, such the dHSIC [59], $\mathrm{JdCov}$
[12], and the $\mathrm{dCov}$ based measure in [52]. Our test performs well
across a variety of data distributions and is especially powerful compared to
the existing tests in contamination models and heavy-tailed distributions.

Finally, in Section 7 we apply our method to a dataset consisting of stock
prices of different companies in USA. Our approach sheds more insights into
the structure of the higher-order dependencies in this data, producing results
that are more interpretable than existing methods. All the implementation of
tests can be found in GitHub.

##  2\. Rank-Based Joint Independence Measures

In this section, we define a rank-based measure of joint independence by
combining higher-order distance covariances with optimal transportation of
measures. We review the relevant concepts regarding higher-order distance
covariance and joint dependence measures in Section 2.2 and OT-based
multivariate ranks in Section 2.3. The rank-based joint dependence measures
are introduced in Section 2.4 and their estimation is discussed in Section
2.5. We begin by introducing some notations.

###  2.1. Notations

Fix $r\geq 2$ and suppose $\bm{X}=(X_{1},\ldots,X_{r})$ is a random vector,
where each $X_{i}$ is a random variable taking values in $\mathbb{R}^{d_{i}}$,
for $1\leq i\leq r$ and $d_{0}=\sum_{i=1}^{r}d_{i}$. The characteristic
function of $X_{i}$ will be denoted by

$\phi_{X_{i}}(t)=\mathbb{E}\left[e^{\iota\langle t,X_{i}\rangle}\right],$

for $t\in\mathbb{R}^{d_{i}}$.

######  Definition 2.1.

The random vector $\bm{X}=(X_{1},\ldots,X_{r})$ is said to be $q$-independent
(for some $2\leq q\leq r$) if for any sub-family
$\\{i_{1},i_{2},\ldots,i_{q}\\}\subset\\{1,2,\ldots,r\\}$ the random variables
$X_{i_{1}},X_{i_{2}},\ldots,X_{i_{q}}$ are mutually independent.

Throughout we will denote

$w_{d_{i}}(t):=\frac{1}{c_{d_{i}}\|t\|^{1+d_{i}}}\quad\text{where
}c_{d_{i}}:=\frac{\pi^{\frac{1+d_{i}}{2}}}{\Gamma\left(\frac{1+d_{i}}{2}\right)},$

for $1\leq i\leq r$. Moreover, for $t_{i}\in\mathbb{R}^{d_{i}}$ define

$\displaystyle\mathrm{d}w_{i}:=w_{d_{i}}(t_{i})\mathrm{d}t_{i}=\frac{\mathrm{d}t_{i}}{c_{d_{i}}\|t_{i}\|^{1+d_{i}}}\quad\text{and
}\quad\mathrm{d}w=\prod_{i=1}^{r}\mathrm{d}w_{i}.$ (2.1)

Finally, for $p\geq 1$, let $\mathcal{P}(\mathbb{R}^{p})$ and
$\mathcal{P}_{ac}(\mathbb{R}^{p})$ denote the collection of all probability
distributions and Lebesgue absolutely continuous probability distributions on
$\mathbb{R}^{p}$, respectively.

###  2.2. Higher-Order Distance Covariance and Joint Dependence Measures

The celebrated distance covariance $(\mathrm{dCov})$ of Székely and Rizzo [70]
is a powerful measure for independence between two random vectors. It has been
widely used for testing of independence [70, 66], feature screening [48] and
general association analysis, including canonical component analysis [79] and
independent component analysis [52]. Recently, there has been efforts to
generalize the notion of $\mathrm{dCov}$ beyond pairwise independence to joint
independence of a collection of $r>2$ random vectors. Towards this,
independently and concurrently, Böttcher et al. [9] and Chakraborty and Zhang
[12] proposed the following higher-order generalization of $\mathrm{dCov}$.
This is referred to as the distance multivariance (by Böttcher et al. [9]) or
the $r$-th order $\mathrm{dCov}$ (by Chakraborty and Zhang [12]).

######  Definition 2.2 ([9, Definition 2.1] and [12, Definition 1]).

The $r$-th order $\mathrm{dCov}$ of $(X_{1},X_{2},\ldots,X_{r})$ is defined as
the positive square root of

$\displaystyle\mathrm{dCov}^{2}(X_{1},\ldots,X_{r})=\int_{\mathbb{R}^{d_{0}}}\left|\mathbb{E}\left[\prod_{i=1}^{r}(\phi_{X_{i}}(t_{i})-e^{\iota\langle
t_{i},X_{i}\rangle})\right]\right|^{2}\mathrm{d}w.$ (2.2)

where $\phi_{X_{i}}(t_{i})=\mathbb{E}\left[e^{\iota\langle
t_{i},X_{i}\rangle}\right]$ for $t_{i}\in\mathbb{R}^{d_{i}}$, and
$\mathrm{d}w$ is as defined in (2.1).

Clearly, $\mathrm{dCov}(X_{1},X_{2},\ldots,X_{r})=0$ whenever the collection
$(X_{1},X_{2},\ldots,X_{r})$ are mutually independent. However, the converse
is only true for $r=2$, in which case (2.2) reduces to the classical
$\mathrm{dCov}$ between two random vectors $X_{1}$ and $X_{2}$ [70]. In other
words, when $r\geq 3$ the joint independence of $(X_{1},X_{2},\ldots,X_{r})$
is not a necessary condition for $\mathrm{dCov}(X_{1},X_{2},\ldots,X_{r})$ to
be zero. (For example, $\mathrm{dCov}(X_{1},X_{2},X_{3})=0$ whenever
$X_{1}\perp\\!\\!\\!\perp(X_{2},X_{3})$.) One can conclude
$(X_{1},X_{2},\ldots,X_{r})$ are jointly independent when
$\mathrm{dCov}(X_{1},X_{2},\ldots,X_{r})=0$ under the additional assumption
that $(X_{1},X_{2},\ldots,X_{r})$ are $(r-1)$-independent (see Definition 2.1
and [9, Theorem 3.4]), but not in general. To characterize joint independence
one has to consider all possible interactions between the $r$ variables
$X_{1},X_{2},\ldots,X_{r}$. This leads to the notion of the total distance
multivariance [9] or the joint distance correlation ($\mathrm{JdCov}$) [12].
Towards this, we need the following notation: For
$S\subseteq\\{1,2,\ldots,r\\}$ with $|S|\geq 2$ denote by
$\bm{X}_{S}=(X_{i})_{i\in S}$ and

$\displaystyle\mathrm{dCov}^{2}(\bm{X}_{S})=\int_{\mathbb{R}^{d_{S}}}\left|\mathbb{E}\left[\prod_{i\in
S}(\phi_{X_{i}}(t_{i})-e^{\iota\langle
t_{i},X_{i}\rangle})\right]\right|^{2}\prod_{i\in S}\mathrm{d}w_{i},$ (2.3)

where $d_{S}:=\sum_{i\in S}d_{i}$. In other words,
$\mathrm{dCov}^{2}(\bm{X}_{S})$ is the $|S|$-th order $\mathrm{dCov}$ for the
variables in $\bm{X}_{S}$.

######  Definition 2.3 ([12, Definition 2]).

Given a vector non-negative weights $\bm{C}=(C_{2},C_{3},\ldots,C_{r})$, the
joint distance covariance $(\mathrm{JdCov})$ of the random vector
$\bm{X}=(X_{1},X_{2},\ldots,X_{r})$ is defined as:

$\displaystyle\mathrm{JdCov}^{2}(\bm{X};\bm{C})=\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\mathrm{dCov}^{2}(\bm{X}_{S}).$ (2.4)

$\mathrm{JdCov}$ completely characterizes the joint independence of
$(X_{1},X_{2},\ldots,X_{r})$, that is,

$\displaystyle\mathrm{JdCov}^{2}(\bm{X};\bm{C})=0\text{ if and only if
}(X_{1},X_{2},\ldots,X_{r})\text{ are mutually independent }$ (2.5)

(see [12, Proposition 3]). Moreover, by choosing all the weights to be equal,
that is, $C_{s}=1$, for $2\leq s\leq r$, one gets the total distance
multivariance as in [9, Definition 2.1]. Chakraborty and Zhang [12] suggests
choosing $C_{s}=c^{r-s}$, for some constant $c>0$, which allows the following
simplier expression that does not require evaluating all the $(2^{r}-r-1)$
dCov terms in (2.4).

######  Lemma 2.4 (Proposition 4 in [12]).

For any $c\geq 0$,

$\displaystyle\mathrm{JdCov}^{2}(\bm{X};c):=\mathrm{JdCov}^{2}(\bm{X};(c^{r-s})_{2\leq
s\leq
r})=\mathbb{E}\left[\prod_{i=1}^{r}(U_{i}(X_{i},X_{i}^{\prime})+c)\right]-c^{r},$

where $(X_{1}^{\prime},\ldots,X_{r}^{\prime})$ be an independent copy of
$\bm{X}=(X_{1},X_{2},\ldots,X_{r})$ and, for $1\leq i\leq r$,

$\displaystyle U_{i}(x,x^{\prime})$
$\displaystyle=\mathbb{E}\|x-X_{i}^{\prime}\|+\mathbb{E}\|X_{i}-x^{\prime}\|-\|x-x^{\prime}\|-\mathbb{E}\|X_{i}-X_{i}^{\prime}\|$
$\displaystyle=\int_{\mathbb{R}^{d_{i}}}\bigg{\\{}(\phi_{X_{i}}(t)-e^{\iota\langle
t,x\rangle})(\phi_{X_{i}}(-t)-e^{-\iota\langle
t,x^{\prime}\rangle})\bigg{\\}}w_{d_{i}}(t)\mathrm{d}t.$

Depending on whether $c>1$ or $c<1$, $\mathrm{JdCov}^{2}(\bm{X};c)$ puts more
or less weights on the lower-order dependence terms, respectively. For
example, if the joint distribution of $\\{X_{1},\ldots,X_{r}\\}$ is known to
be Gaussian, where mutual independence is equivalent to pairwise independence,
larger values of $c$ should be considered. Otherwise, a smaller $c<1$ makes
more sense.

###  2.3. Multivariate Ranks Based on Optimal Transport

To motivate the notion of multivariate ranks based on optimal transport, let
us recall some fundamental facts about univariate ranks. Suppose that
$X\sim\nu\in\mathcal{P}_{ac}(\mathbb{R})$ is a random variable with
distribution function $F$. Here, the distribution function itself serves as
the one-dimensional population rank function, which has the property that
$F(X)\sim\mathrm{Unif}([0,1])$, that is, $F$ transports $\mu$ to
$\mathrm{Unif}([0,1])$, the uniform distribution on $[0,1]$. In fact, when
$\mu$ has finite second moment, it can be shown that $F$ is the almost
everywhere unique map that transports $\mu$ to $\mathrm{Unif}([0,1])$ and
minimizes the expected squared-error cost, that is,

$F=\mathrm{arg}\min_{T:T(X)\sim\mathrm{Unif}([0,1])}\mathbb{E}_{X\sim\mu}[(X-T(X))^{2}],$

where the minimization is over all functions $T$ that transport the
distribution of $X$ to the uniform distribution on $[0,1]$. This shows that in
dimension 1, ranks can be thought of as the univariate analogue of the
celebrated Monge transportation problem [54, 73]. Hallin et al. [31, 30] used
this interpretation of ranks and proposed the following multivariate
generalization.

######  Definition 2.5 ([18, 31, 30]).

Given a random variable $X\sim\mu$ and a pre-specified reference distribution
$\nu$, with $\mu,\nu\in\mathcal{P}_{ac}(\mathbb{R}^{d})$, the multivariate
population rank map $R_{\mu}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is
defined as:

$\displaystyle
R_{\mu}:=\mathrm{arg}\min_{T:T(X)\sim\nu}\mathbb{E}_{X\sim\mu}[||X-T(X)||^{2}],$
(2.6)

where the minimum is over all functions
$T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that transports the distribution
of $X$ to $\mathbb{R}^{d}$ and $||\cdot||$ denotes the usual Euclidean norm in
$\mathbb{R}^{d}$.222Note that, even though the definition in (2.6) is not
meaningful when $X\sim\mu$ does not have finite second moment, OT-based
multivariate ranks can still be defined using the Brenier-McCann theorem [53,
10], which says that there exists an almost everywhere unique map
$R_{\mu}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, which is the gradient of a
convex function and satisfies $R_{\mu,\nu}(X)\sim\nu$. This notion coincides
with that in Definition 2.5, whenever $\mu$ and $\nu$ have finite second
moment.

To define the empirical analogue of the population multivariate rank we begin
with a ‘natural’ discretization
$\mathcal{H}_{N}^{d}:=\\{h_{1}^{d},\ldots,h_{N}^{d}\\}$ of the pre-specified
reference distribution $\nu\in\mathcal{P}_{ac}(\mathbb{R}^{d})$, that is, the
empirical measure $\nu_{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{h_{i}^{d}}$
converges converges weakly to $\nu$. The natural choice for $d=1$ is
$\mathcal{H}_{N}^{1}=\\{1/N,2/N,\ldots,1\\}$, whose empirical distribution
converges to $\mathrm{Unif}([0,1])$. For higher dimensions one can choose
$\mathcal{H}_{N}^{d}$ as $N$ i.i.d. points from $\nu$ or, more commonly, a
deterministic quasi-Monte Carlo sequence such as the Halton sequence [18, 38].
Then, given i.i.d. samples $X_{1},X_{2},\ldots,X_{N}$ from a distribution
$\mu\in\mathcal{P}_{ac}(\mathbb{R}^{d})$, the multivariate empirical rank map
$\hat{R}(X_{i})$ is the optimal transport map from the empirical distribution
of the data $\hat{\mu}_{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{i}}$ to
$\hat{\mu}_{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{h_{i}^{d}}$. In other words,

$\displaystyle\hat{R}(X_{i})=h_{\hat{\sigma}_{N}(i)}^{d},$ (2.7)

where

$\displaystyle\mathbf{\hat{\sigma}}_{N}:=\mathrm{arg}\min_{\sigma\in
S_{N}}\sum_{i=1}^{N}||X_{i}-h_{{\sigma(i)}}^{d}||^{2},$

and $S_{N}$ denotes the set of all $N!$ permutations of $\\{1,2,\ldots,N\\}$.
This is a minimum bipartite matching problem (also known as the assignment
problem), which is a fundamental problem in combinatorial optimization that
can be solved exactly $O(N^{3})$ time [40] and approximately in near-linear
time (see [1] and the references therein).

A remarkable feature of the optimal transport based multivariate empirical
ranks defined in (2.7) is that they mimic the distribution-free property of
univariate ranks [18, 31].

######  Proposition 2.6 ([31, Proposition 1.6.1] and [18, Proposition 2.2]).

Suppose that $X_{1},\ldots,X_{N}$ are i.i.d. samples from a distribution
$\mu\in\mathcal{P}_{ac}(\mathbb{R}^{d})$. Then the vector of empirical ranks

$(\hat{R}(X_{1}),\hat{R}(X_{2}),\ldots,\hat{R}(X_{N}))$

is uniformly distributed over the $N!$ permutations of the fixed grid
$\mathcal{H}_{N}^{d}$.

The result above suggests a general strategy for constructing non-parametric
distribution-free tests, by replacing the data points with their empirical
ranks (appropriately defined depending on the testing problem). This strategy
was recently adopted in [18, 67] to construct distribution-free,
computationally efficient, and universally consistent multivariate two-sample
and independence tests. For other interesting properties and applications of
optimal transport based ranks see [16, 24, 33, 32] and the references therein.

######  Remark 2.7.

Common choices of the reference distribution $\nu$ are
$\mathrm{Unif}([0,1]^{d})$, the uniform distribution on the $d$-dimensional
cube $[0,1]^{d}$ [18], the spherical uniform distribution $U_{d}$, which is
the product of the uniform distribution on $[0,1)$ and the uniform
distribution on the unit sphere $S^{d-1}$ [30, 67], and
$N(\bm{0},\bm{I}_{d})$, the $d$-dimensional standard Gaussian [20]. Hereafter,
for concreteness we will choose $\nu=\mathrm{Unif}([0,1]^{d})$ and
$\mathcal{H}_{N}^{d}:=\\{h_{1}^{d},\ldots,h_{N}^{d}\\}$ will denote the Halton
sequence corresponding to this distribution. However, our results continue to
hold for any reference distribution $\nu$ with finite moments, which, in
particular, includes the spherical uniform and the standard multivariate
Gaussian.

###  2.4. Rank-Based Joint Dependence Measures

We are now ready to define the rank-based counterparts of the $d$-th order
$\mathrm{dCov}$ and $\mathrm{JdCov}$. Towards this, suppose, as before,
$(X_{1},\ldots,X_{r})$ is a $r$-dimensional random vector, where each $X_{i}$
is a random variable taking in $\mathbb{R}^{d_{i}}$ with distribution
$\mu_{i}$, for $1\leq i\leq r$ and $d_{0}=\sum_{i=1}^{r}d_{i}$. Throughout we
will assume $\mu_{i}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{i}})$ and denote the
population rank map of $\mu_{i}$ by $R_{\mu_{i}}$ (recall Definition 2.5), for
$1\leq i\leq r$.

######  Definition 2.8.

The $r$-th order rank $\mathrm{dCov}$ $(\mathrm{RdCov})$ of
$(X_{1},X_{2},\ldots,X_{r})$ is defined as the positive square root of

$\displaystyle\mathrm{RdCov}^{2}(X_{1},\ldots,X_{r})=\int_{\mathbb{R}^{d_{0}}}\left|\mathbb{E}\left[\prod_{i=1}^{r}\left\\{\mathbb{E}\left(e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right)-e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right\\}\right]\right|\mathrm{d}w$ (2.8)

where $t_{i}\in\mathbb{R}^{d_{i}}$, for $1\leq i\leq r$, and $\mathrm{d}w$ is
as defined in (2.1).

Note that the $r$-th order $\mathrm{RdCov}$ is obtained from the $r$-th order
$\mathrm{dCov}$ by replacing $(X_{1},X_{2},\ldots,X_{r})$ with their
population rank maps
$(R_{\mu_{1}}(X_{1}),R_{\mu_{2}}(X_{2}),\ldots,R_{\mu_{r}}(X_{r}))$. For
$r=2$, (2.8) reduces to the rank distance covariance of $(X_{1},X_{2})$
defined in [18, 67], which, like the classical $\mathrm{dCov}$, is zero if and
only if $X_{1}\perp\\!\\!\\!\perp X_{2}$. For the $r\geq 3$, as in
$\mathrm{JdCov}$, we need to consider all the interactions between the $r$
variables $(X_{1},X_{2},\ldots,X_{r})$ to characterize the joint indepedence.
This leads to our the notion of rank joint distance covariance
($\mathrm{RJdCov}$):

######  Definition 2.9.

Given a vector non-negative weights $\bm{C}=(C_{2},C_{2},\ldots,C_{r})$, the
rank joint distance covariance $(\mathrm{RJdCov})$ of the random vector
$\bm{X}=(X_{1},X_{2},\ldots,X_{r})$ is defined as:

$\displaystyle\mathrm{RJdCov}^{2}(\bm{X};\bm{C})=\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\mathrm{RdCov}^{2}(\bm{X}_{S}),$ (2.9)

where

$\displaystyle\mathrm{RdCov}^{2}(\bm{X}_{S})=\int_{\mathbb{R}^{d_{S}}}\left|\mathbb{E}\left[\prod_{i\in
S}\left\\{\mathbb{E}\left(e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right)-e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right\\}\right]\right|\prod_{i\in
S}\mathrm{d}w_{i}.$ (2.10)

Moreover, for any $c\geq 0$,

$\displaystyle\mathrm{RJdCov}^{2}(\bm{X};c):=\mathrm{RJdCov}^{2}(\bm{X};(c^{r-s})_{2\leq
s\leq
r})=\sum_{s=2}^{r}c^{r-s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\mathrm{RdCov}^{2}(\bm{X}_{S}).$ (2.11)

$($Note that
$\mathrm{RdCov}^{2}(X_{1},\ldots,X_{r})=\mathrm{RJdCov}^{2}(\bm{X};0)$.$)$

As in Lemma 2.4, $\mathrm{RJdCov}^{2}(\bm{X};c)$ has the following compact
representation:

$\displaystyle\mathrm{RJdCov}^{2}(\bm{X};c):=\mathbb{E}\left[\prod_{i=1}^{r}(W_{i}(R_{\mu_{i}}(X_{i}),R_{\mu_{i}}(X_{i}^{\prime}))+c)\right]-c^{r},$
(2.12)

where $(X_{1}^{\prime},\ldots,X_{r}^{\prime})$ be an independent copy of
$\bm{X}=(X_{1},X_{2},\ldots,X_{r})$ and, for $1\leq i\leq r$,

$\displaystyle W_{i}(x,x^{\prime})$
$\displaystyle=\mathbb{E}\|x-R_{\mu_{i}}(X_{i}^{\prime})\|+\mathbb{E}\|R_{\mu_{i}}(X_{i})-x^{\prime}\|-\|x-x^{\prime}\|-\mathbb{E}\|R_{\mu_{i}}(X_{i})-R_{\mu_{i}}(X_{i}^{\prime})\|.$
(2.13)

Note that $R_{\mu_{1}}(X_{1}),R_{\mu_{2}}(X_{2}),\ldots,R_{\mu_{r}}(X_{r})$
are i.i.d. $\mathrm{Unif}([0,1]^{d})$ when the variables
$X_{1},X_{2},\ldots,X_{r}$ are mutually independent (by the definition of the
optimal transport maps). Therefore, the distribution of
$\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$ does not depend on the marginal
distributions of $X_{1},X_{2},\ldots,X_{r}$ under mutual independence.
Moreover, $\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$ characterizes the joint
independence of $X_{1},X_{2},\ldots,X_{r}$.

######  Proposition 2.10.

Let $\mathrm{RdCov}$ and $\mathrm{RJdCov}$ be as defined in (2.10) and (2.9),
respectively. Then the following hold:

  * $(1)$

For any $S\subseteq\\{1,2,\ldots,r\\}$ with $|S|\geq 2$,
$\mathrm{RdCov}(\bm{X}_{S})=0$ if the variables $(X_{i})_{i\in S}$ are
independent. Moreover, the converse holds whenever the variables
$(X_{i})_{i\in S}$ are $|S|-1$ independent.

  * $(2)$

For any vector positive weights $\bm{C}$,
$\mathrm{RJdCov}^{2}(\bm{X};\bm{C})=0$ if and only if
$(X_{1},X_{2},\ldots,X_{r})$ are mutually independent.

###### Proof.

Note that $\\{R_{\mu_{i}}(X_{i})\\}_{i\in S}$ are independent when
$(X_{i})_{i\in S}$ are mutually independent, which implies
$\mathrm{RdCov}^{2}(\bm{X}_{S})=0$, by definition. For the converse, note that
$\mathrm{RJdCov}^{2}(\bm{X}_{S})=0$ and $\\{R_{\mu_{i}}(X_{i})\\}_{i\in S}$
are $|S|-1$ independent implies $\\{R_{\mu_{i}}(X_{i})\\}_{i\in S}$ are
mutually independent by [9, Theorem 3.4]. Then, since
$\mu_{i}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{i}})$, by McCann’s theorem [53]
(see also [18, Proposition 2.1]) there exists a measurable function
$Q_{\mu_{i}}$ such that $Q_{\mu_{i}}(R_{\mu_{i}}(X_{i}))=X_{i}$ almost
everywhere $\mu_{i}$, for $1\leq i\leq r$. This implies,
$X_{1},X_{2},\ldots,X_{r}$ are mutually independent. This completes the proof
of (1).

For (2) note that $\mathrm{RJdCov}^{2}(\bm{X};\bm{C})=0$ if and only if
$\mathrm{RdCov}^{2}(\bm{X}_{S})=0$ for all $S\subseteq\\{1,2,\ldots,r\\}$ with
$|S|\geq 2$. The result then follows from (1). ∎

###  2.5. Consistent Estimation of $\mathrm{RJdCov}$

In this section we discuss how $\mathrm{RdCov}$ and $\mathrm{RJdCov}$ can be
consistently estimated given $n$ samples i.i.d. $\\{\bm{X}_{a}\\}_{a=1}^{n}$,
with $\bm{X}_{a}=(X_{1}^{(a)},\ldots,X_{r}^{(a)})$, from a distribution
$\mu\in\mathcal{P}_{ac}(\mathbb{R}^{d_{0}})$. The natural plug-in estimator of
$\mathrm{RdCov}^{2}(\bm{X}_{S})$, for $S\subseteq\\{1,2,\ldots,r\\}$ with
$|S|\geq 2$, is:

$\displaystyle\mathrm{RdCov}^{2}_{n}(\bm{X}_{S}):=\int_{\mathbb{R}^{d_{S}}}\left|\frac{1}{n}\sum_{b=1}^{n}\left\\{\prod_{i\in
S}\left(\frac{1}{n}\sum_{a=1}^{n}e^{\iota\langle
t_{i},\hat{R}_{i}(X_{i}^{(a)})\rangle}-e^{\iota\langle
t_{i},\hat{R}_{i}(X_{i}^{(b)})\rangle}\right)\right\\}\right|^{2}\prod_{i\in
S}\mathrm{d}w_{i},$ (2.14)

where $\hat{R}_{i}$ is the empirical rank map for the $i$-th marginal
distribution, that is, the optimal transport map from
$\frac{1}{n}\sum_{a=1}^{n}\delta_{X_{i}^{(a)}}$ to
$\frac{1}{n}\sum_{a=1}^{n}\delta_{h_{a}^{d_{i}}}$, with
$\mathcal{H}_{n}^{d_{i}}=\\{h_{1}^{d_{i}},h_{2}^{d_{i}},\ldots,h_{n}^{d_{i}}\\}$
the Halton sequence in $[0,1]^{d_{i}}$, for $1\leq i\leq r$. Using the
representation in (2.12) and (2.13) the estimate in (2.14) can be written as:

$\displaystyle\mathrm{RdCov}^{2}_{n}(\bm{X}_{S})$
$\displaystyle=\frac{1}{n^{2}}\sum_{1\leq a,b\leq n}\prod_{i\in
S}\hat{\mathcal{E}}_{i}(a,b).$ (2.15)

where

$\displaystyle\hat{\mathcal{E}}_{i}(a,b)$
$\displaystyle:=\frac{1}{n}\sum_{v=1}^{n}\|\hat{R}_{i}(X_{i}^{(a)})-\hat{R}_{i}(X_{i}^{(v)})\|+\frac{1}{n}\sum_{u=1}^{n}\|\hat{R}_{i}(X_{i}^{(u)})-\hat{R}_{i}(X_{i}^{(b)})\|$
$\displaystyle\hskip
90.3375pt-\|\hat{R}_{i}(X_{i}^{(a)})-\hat{R}_{i}(X_{i}^{(b)})\|-\frac{1}{n^{2}}\sum_{1\leq
u,v\leq n}\|\hat{R}_{i}(X_{i}^{(u)})-\hat{R}_{i}(X_{i}^{(v)})\|,$ (2.16)

for $1\leq a,b\leq n$. Consequently, the natural plug-in estimate of (2.9) is

$\displaystyle\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})=\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\mathrm{RdCov}^{2}_{n}(\bm{X}_{S}),$ (2.17)

and that of (2.11) (using the representation in (2.12) is

$\displaystyle\mathrm{RJdCov}^{2}_{n}(\bm{X};c)$
$\displaystyle:=\frac{1}{n^{2}}\sum_{1\leq a,b\leq
n}\prod_{i=1}^{r}\left(\hat{\mathcal{E}}_{i}(a,b)+c\right)-c^{r}.$

######  Remark 2.11.

For $S\subseteq\\{1,2,\ldots,r\\}$ with $|S|\geq 2$, consider the function
$\theta_{S}:(\mathbb{R}^{d_{S}})^{n}\rightarrow\mathbb{R}_{\geq 0}$

$\displaystyle\theta_{S}((x_{i}^{(1)})_{i\in S},\ldots,(x_{i}^{(n)})_{i\in
S}):=\int_{\mathbb{R}^{d_{S}}}\left|\frac{1}{n}\sum_{b=1}^{n}\left\\{\prod_{i\in
S}\left(\frac{1}{n}\sum_{a=1}^{n}e^{\iota\langle
t_{i},x_{i}^{(a)})\rangle}-e^{\iota\langle
t_{i},x_{i}^{(b)}\rangle}\right)\right\\}\right|^{2}\prod_{i\in
S}\mathrm{d}w_{i},$ (2.18)

where $x_{i}^{(a)}\in\mathbb{R}^{d_{i}}$, for $1\leq i\leq r$ and, hence,
$(x_{i}^{(a)})_{i\in S}\in\mathbb{R}^{d_{S}}$, for $1\leq a\leq n$. Evaluating
$\theta_{S}$ at the data points gives us the natural plug-in estimate of
$\mathrm{dCov}(\bm{X}_{S})$ (recall (2.3)) which is denoted by:

$\mathrm{dCov}_{n}(\bm{X}_{S}):=\theta_{S}(\bm{X}_{S}^{(1)},\ldots,\bm{X}_{S}^{(n)}),$

where $\bm{X}_{S}^{(a)}:=(X_{i}^{(a)})_{i\in S}\in\mathbb{R}^{d_{S}}$, for
$1\leq a\leq n$. On the other hand, evaluating $\theta_{S}$ on the the rank
transformed data gives us the estimate of $\mathrm{RdCov}$ in (2.14).
Specifically, recalling (2.14) note that

$\displaystyle\mathrm{RdCov}_{n}(\bm{X}_{S})=\theta_{S}(\hat{\bm{R}}_{S}^{(1)},\ldots,\hat{\bm{R}}_{S}^{(n)}),$
(2.19)

where $\hat{\bm{R}}_{S}^{(a)}:=(\hat{R}_{i}(X_{i}^{(a)}))_{i\in
S}\in\mathbb{R}^{d_{S}}$, for $1\leq a\leq n$. The function $\bm{\theta}_{S}$
will play an important role in the calibration of the independence tests on
$\mathrm{RJdCov}_{n}$ discussed in Section 3.

To establish the consistency of $\mathrm{RJdCov}_{n}$ we will assume the
following:

######  Assumption 1.

For every $1\leq i\leq p$, the empirical distribution of
$\mathcal{H}_{n}^{d_{i}}=\\{h_{1}^{d_{i}},h_{2}^{d_{i}},\ldots,h_{n}^{d_{i}}\\}$
converges weakly to $\mathrm{Unif}([0,1]^{d_{i}})$.

As mentioned before, Assumption 1 holds whenever $\mathcal{H}_{n}^{d_{i}}$ is
the Halton sequence in $[0,1]^{d_{i}}$ (which is our default choice), for a
random sample $n$ i.i.d. $\mathrm{Unif}([0,1]^{d_{i}})$ points, or other
quasi-Monte Carlo sequences as well (see [18, Appendix D] for a discussion on
the various of choices $\mathcal{H}_{n}^{d_{i}}$). Under this assumption we
have the following result:

######  Theorem 2.12.

Suppose Assumption 1 holds. Then for any $S\subseteq\\{1,2,\ldots,d\\}$ with
$|S|\geq 2$,

$\mathrm{RdCov}^{2}_{n}(\bm{X}_{S})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\mathrm{RdCov}^{2}(\bm{X}_{S}).$

Consequently, $\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$.

The proof of Theorem 2.12 is given in Appendix B. One of the highlights of the
above result is that it holds without any moment assumptions on the data
generating distribution. This is in contrast to consistency results for
$\mathrm{JdCov}$ which require the distribution to satisfy certain moment
conditions. Consequently, the rank-based joint dependence measures are more
robust to heavy-tail distributions, as will been seen from the simulations in
Section 6.

##  3\. Distribution-Free Joint Independence Testing

In section we discuss new distribution-free tests for the multivariate joint
independence testing. Given i.i.d.samples $\\{\bm{X}_{a}\\}_{a=1}^{n}$, with
$\bm{X}_{a}=(X_{1}^{(a)},\ldots,X_{r}^{(a)})$, from a distribution
$\mu\in\mathcal{P}_{ac}(\mathbb{R}^{d_{0}})$ with marginal distributions
$\mu_{1},\mu_{2},\ldots,\mu_{r}$ such that
$\mu_{i}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{i}})$, for $1\leq i\leq r$,
consider the joint independence testing problem:

$\displaystyle
H_{0}:\mu=\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}\quad\text{ versus
}\quad H_{1}:\mu\neq\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}.$ (3.1)

Our test for (3.1) will be based on rejecting $H_{0}$ for ‘large’ values of
$\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$, for a given choice non-negative weights
$\bm{C}=(C_{2},C_{3},\ldots,C_{r})$. Note that the vector of ranks
$\\{\hat{R}_{i}(X_{i}^{(a)})\\}_{1\leq a\leq n}$ are distributed uniformly
over the $n!$ permutations of $\mathcal{H}_{n}^{d_{i}}$, for $1\leq i\leq r$
(by Proposition 2.6) and under $H_{0}$ independent over $1\leq i\leq r$.
Hence, the distribution of $\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$ under $H_{0}$
does not depend on the marginal distributions
$\mu_{1},\mu_{2},\ldots,\mu_{r}$.

######  Proposition 3.1 (Distribution-freeness under $H_{0}$).

Under $H_{0}$, the distribution of $\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})$ is
universal, that is, it is free of $\mu_{1},\ldots,\mu_{r}$.

The above result automatically gives a finite sample distribution-free
independence test that uniformly controls the Type-I error for (3.1). To this
end, fix a level $\alpha\in(0,1)$ and let $c_{\alpha,n}$ denote the upper
$\alpha$ quantile of the universal distribution in Proposition 3.1. Consider
the test function:

$\phi_{n}(\bm{C}):=\bm{1}\left\\{\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\geq
c_{\alpha,n}\right\\}.$ (3.2)

This test is exactly distribution-free for all $n\geq 1$ and uniformly of
level $\alpha$ under $H_{0}$,333Strictly speaking, to guarantee exact level
$\alpha$, we have to randomize $\phi_{n}(\bm{C})$, as the exact distribution
of $\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})$ is discrete. However, this makes
no practical difference unless $n$ is very small. that is,

$\sup_{\mu=\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}}\mathbb{E}_{\mu}\left[\phi_{n}(\bm{C})\right]=\alpha.$
(3.3)

Moreover, since $\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})$ is a consistent
estimate of $\mathrm{RJdCov}^{2}(\bm{X};\bm{C})$, which characterizes joint
independence, the test $\phi_{n}(\bm{C})$ is consistent against all fixed
alternatives. This is summarized in the following proposition (see Appendix B
for the proof).

######  Proposition 3.2 (Universal consistency).

Suppose Assumption 1 holds. Then the test $\phi_{n}(\bm{C})$ is universally
consistent, that is, for any
$\mu\neq\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}$,

$\lim_{n\rightarrow\infty}\mathbb{E}_{\mu}\left[\phi_{n}(\bm{C})\right]=1.$
(3.4)

######  Remark 3.3.

Note that, while $\mathrm{RJdCov}_{n}(\bm{X};\bm{C})$ provides a distribution-
free, universally consistent test for joint independence, the individual
$\mathrm{RdCov}(\bm{X}_{S})$ can be used for testing the joint independence of
the subset of variables $\bm{X}_{S}$ given that they are $|S|-1$ independent.
This is because $\mathrm{RdCov}_{n}(\bm{X}_{S})=0$ if and if only if the
variables $\bm{X}_{S}$ are mutually independent, provided they are $|S|-1$
independent. For example, if we are interested in testing pairwise
independence:

$H_{0}:X_{1},X_{2},\ldots,X_{r}\text{ are pairwise independent }\quad\text{
versus }\quad H_{1}:\text{ not }H_{0},$

then the test which rejects $H_{0}$ for ‘large’ values of

$\sum_{1\leq i<j\leq r}\mathrm{RdCov}^{2}_{n}(X_{i},X_{j}),$

will be distribution-free and consistent. This is, in fact, the rank-analogue
of the test based on pairwise distance covariances studied in [76]. Similarly,
if we know that 3 variables $(X_{1},X_{2},X_{3})$ are pairwise independent,
then $\mathrm{RdCov}^{2}_{n}(X_{1},X_{2},X_{3})$ can be used to obtain a
distribution-free, consistent test for the mutual independence of
$(X_{1},X_{2},X_{3})$. Another attractive property of $\mathrm{RdCov}$ is that
the collection $\mathrm{RdCov}_{n}(\bm{X}_{S})$ is asymptotically independent
over $S\subseteq\\{1,2,\ldots,d\\}$ under the null hypothesis of joint
independence (see Theorem 3.7), which can be leveraged to learn the higher-
order dependencies among the variables $X_{1},X_{2},\ldots,X_{r}$.

######  Remark 3.4.

As mentioned in the Introduction, another related measure for joint
independence based on the second-order $\mathrm{dCov}$ was proposed by
Matteson and Tsay [52], and its rank version was discussed in [18]. Although
these measures characteristic joint independence, unlike
$\mathrm{JdCov}/\mathrm{RJdCov}$ it is unable to capture the dependence
structure among the variables. In Section 6 we compare our method with test
proposed in [52] in simulations.

###  3.1. Asymptotic Null Distribution

In this section we will derive the asymptotic null distribution of the
collection of $2^{r}-r-1$ random variables:

$\displaystyle\\{\mathrm{RdCov}_{n}(\bm{X}_{S}):S\in\mathcal{T}\\},$ (3.5)

where $\mathcal{T}:=\\{S\subseteq\\{1,2,\ldots,r\\}\text{ such that }|S|\geq
2\\}$, and consequently that of $\mathrm{RJdCov}_{n}(\bm{X};\bm{C})$. To
describe the limit we need the following definition:

######  Definition 3.5.

Denote by $\\{Z_{S}\\}_{S\in\mathcal{T}}$ a collection of mutually independent
complex-valued Gaussian processes such that, for each $S\in\mathcal{T}$,

$\displaystyle\\{Z_{S}(\bm{t}):\bm{t}\in\mathbb{R}^{d_{S}}\\}$ (3.6)

a complex-valued Gaussian process indexed by $\mathbb{R}^{d_{S}}$ (recall
$d_{S}=\sum_{i\in S}d_{i}$) with zero mean and covariance function:

$\displaystyle
C_{S}(\bm{t},\bm{t}^{\prime}):=\operatorname{Cov}[Z_{S}(\bm{t}),Z_{S}(\bm{t}^{\prime})]:=\prod_{i\in
S}\left(\mathbb{E}\left[e^{\iota\langle
t_{i}-t_{i}^{\prime},U_{i}\rangle}\right]-\mathbb{E}\left[e^{\iota\langle
t_{i},U_{i}\rangle}\right]\mathbb{E}\left[e^{\iota\langle-
t_{i}^{\prime},U_{i}\rangle}\right]\right),$ (3.7)

where

  * •

$\bm{t}=(t_{i})_{i\in S}\in\mathbb{R}^{d_{S}}$ and
$\bm{t}^{\prime}=(t_{i}^{\prime})_{i\in S}\in\mathbb{R}^{d_{S}}$, with
$t_{i},t_{i}^{\prime}\in\mathbb{R}^{d_{i}}$, for $i\in S$, and

  * •

$U_{1},U_{2},\ldots,U_{r}$ are independent with
$U_{i}\sim\mathrm{Unif}([0,1]^{d_{i}})$, for $1\leq i\leq r$.

######  Remark 3.6.

The covariance function (3.7) can be expressed in closed form using the fact,

$\mathbb{E}\left[e^{\iota\langle
t_{i},U_{i}\rangle}\right]=\prod_{j=1}^{d_{i}}\frac{e^{\iota t_{ij}}-1}{\iota
t_{ij}},$

for $U_{i}\sim\mathrm{Unif}([0,1]^{d_{i}})$ and $t_{i}=(t_{ij})_{1\leq j\leq
d_{i}}\in\mathbb{R}^{d_{i}}$. Hence,

$C_{S}(\bm{t},\bm{t}^{\prime})=\prod_{i\in
S}\left(\prod_{j=1}^{d_{i}}\frac{e^{\iota(t_{ij}-t_{ij}^{\prime})}-1}{\iota(t_{ij}-t_{ij}^{\prime})}-\prod_{j=1}^{d_{i}}\frac{e^{\iota
t_{ij}}-1}{t_{ij}}\prod_{j=1}^{d_{i}}\frac{e^{-\iota
t_{ij}^{\prime}}-1}{t_{ij}^{\prime}}\right),$

where $t_{i}=(t_{ij})_{1\leq j\leq d_{i}}\in\mathbb{R}^{d_{i}}$ and
$t_{i}^{\prime}=(t_{ij}^{\prime})_{1\leq j\leq d_{i}}\in\mathbb{R}^{d_{i}}$.

With the above definition we can now describe the asymptotic distribution of
the collection (3.5) under the null hypothesis of joint independence.

######  Theorem 3.7.

Suppose Assumption 1 holds. Then under $H_{0}$ as in (3.1), as
$n\rightarrow\infty$,

$\displaystyle\\{n\mathrm{RdCov}_{n}(\bm{X}_{S}):S\in\mathcal{T}\\}\stackrel{{\scriptstyle
D}}{{\to}}\left\\{\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}\right\\}_{S\in\mathcal{T}},$ (3.8)

for $Z_{S}(\bm{t})$ as defined in (3.6). Consequently, for any vector of non-
negative weights $\bm{C}=(C_{1},C_{2},\ldots,C_{r})$, under $H_{0}$ as in
(3.1),

$\displaystyle n\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\stackrel{{\scriptstyle
D}}{{\to}}\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}.$ (3.9)

The proof of Theorem 3.7 is given in Appendix C. One of the challenges in
dealing with $\mathrm{RdCov}_{n}$ and $\mathrm{RJdCov}_{n}$ is that they are
functions of dependent multivariate ranks. To solve the problem, we develop a
version of the classical Hoeffiding’s combinatorial central limit theorem for
$r$-tensors indexed by multiple independent permutations (see Theorem A.1 in
Appendix A), a result which might be more broadly useful in the analysis of
other nonparametric tests. Our proof of the combinatorial CLT uses Stein’s
method based on exchangeable pairs, combined with techniques from empirical
process theory gives us the result in Theorem 3.7.

######  Remark 3.8.

(Equivalent representation of the limiting distribution) By considering the
process $Z_{S}(\bm{t})$ as a random element in $L^{2}(\mathbb{R}^{d_{S}})$ and
an eigen-expansion of the integral operator on associated with the covariance
kernel (3.7) the limit in (3.8) can be expressed as:

$\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}\stackrel{{\scriptstyle
D}}{{=}}\sum_{j=1}^{\infty}\lambda_{j}Z_{j}^{2},$

where $\\{\lambda_{j}\\}_{j\geq 1}$ are positive constants and
$Z_{1},Z_{2},\ldots,$ are i.i.d. $N(0,1)$ (see [46, Chapter 1, Section 2]).
Hence, by the independence of the processes $\\{Z_{S}\\}_{S\in\mathcal{T}}$,
the limit in (3.9) can be represented as an infinite sum weighted chi-squares:

$n\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\stackrel{{\scriptstyle
D}}{{\to}}\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}\stackrel{{\scriptstyle
D}}{{=}}\sum_{j=1}^{\infty}\lambda_{j}^{\prime}Z_{j}^{2},$

for a sequence of positive constants $\\{\lambda_{j}^{\prime}\\}_{j\geq 1}$
and $Z_{1},Z_{2},\ldots,$ are i.i.d. $N(0,1)$.

######  Remark 3.9.

Another interesting consequence of Theorem 3.7 is that collection
$\\{\mathrm{RdCov}_{n}(\bm{X}_{S}):S\in\mathcal{T}\\}$ is asymptotically
independent. This can be attributed to the representation of
$\mathrm{RdCov}_{n}(\bm{X}_{S})$ as a product of zero mean (under $H_{0}$)
random variables (recall (2.14)). This property is shared by the
$\mathrm{JdCov}$ (equivalently, the distance multivariance) which has a
similar product structure (see [9, Theorem 4.10]).

###  3.2. Approximating the Cut-off and Finite Sample Properties

Note that the limit in (3.8), as expected, does not depend on the distribution
of the data or the specific choices of
$\mathcal{H}_{n}^{d_{1}},\mathcal{H}_{n}^{d_{2}},\ldots,\mathcal{H}_{n}^{d_{r}}$.
This is a manifestation of the distribution-free property of the multivariate
ranks. Hence, if $c_{\alpha}$ denotes the $(1-\alpha)$-th quantile of the
distribution in the RHS of (3.9), the test which rejects $H_{0}$ when

$\displaystyle\phi^{\mathrm{asymp}}(\bm{C})=\bm{1}\\{n\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})>c_{\alpha}\\},$
(3.10)

will have asymptotic level $\alpha$ and universally consistent for the
hypothesis (3.1). Although this test is asymptotically valid, there is, in
general, no tractable form of the quantile $c_{\alpha}$ for the distribution
(3.9). Nevertheless, using the distribution-free property of the multivariate
ranks we can still approximate the quantiles of (3.8) (and consequently (3.9))
in a data-agnostic manner. This is because $\mathrm{RdCov}_{n}(\bm{X}_{S})$ is
a function of the multivariate ranks (recall (2.14) and (2.19)), and the ranks
are independent and the uniformly distributed over the $n!$ permutations of
$\mathcal{H}_{n}^{d_{i}}$, for $1\leq i\leq r$ (by Proposition 2.6).
Therefore, we can compute an approximate quantiles of
$\mathrm{RdCov}_{n}(\bm{X};\bm{C})$ as follows: For $1\leq b\leq B$ repeat the
following two steps:

  * $(1)$

Generate $r$ i.i.d. uniform random permutations
$\bm{\pi}:=\\{\pi_{1},\pi_{2},\ldots,\pi_{r}\\}$ from $S_{n}$.

  * $(2)$

Compute the value of the statistic $\mathrm{RdCov}_{n}(\bm{X};\bm{C})$ on the
permuted Halton sequences
$\pi_{1}\mathcal{H}_{n}^{d_{1}},\pi_{2}\mathcal{H}_{n}^{d_{2}},\ldots,\pi_{r}\mathcal{H}_{n}^{d_{r}}$,
where

$\pi_{i}\mathcal{H}_{n}^{d_{i}}:=\\{h_{\pi_{i}(1)}^{d_{i}},h_{\pi_{i}(2)}^{d_{i}},\ldots,h_{\pi_{i}(n)}^{d_{i}}\\},$

for $1\leq i\leq r$. More precisely, we evaluate

$\displaystyle\mathrm{RdCov}_{n}^{(b)}(\bm{X};\bm{C})=\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\theta_{S}(\bm{h}_{S,\bm{\pi}}^{(1)},\ldots,\bm{h}_{S,\bm{\pi}}^{(n)}),$
(3.11)

where $\theta_{S}$ is the function in (2.18) and
$\bm{h}_{S,\bm{\pi}}^{(a)}:=(h_{\pi_{i}(a)}^{d_{i}})_{i\in
S}\in\mathbb{R}^{d_{S}}$, for $1\leq a\leq n$.

Then the permutation $(1-\alpha)$-th quantile of
$\mathrm{RdCov}_{n}(\bm{X};\bm{C})$ is obtained as:

$\displaystyle
c_{n,M,\alpha}:=\min\left\\{\mathrm{RdCov}_{n}^{(s)}(\bm{X};\bm{C}):\frac{1}{B}\sum_{b=1}^{B}\bm{1}\left\\{\mathrm{RdCov}_{n}^{(b)}(\bm{X};\bm{C})>\mathrm{RdCov}_{n}^{(s)}(\bm{X};\bm{C})\right\\}\leq\alpha\right\\}.$
(3.12)

To compute the permutation $p$-value, let $R$ be the rank of
$\mathrm{RdCov}_{n}(\bm{X};\bm{C})$ (the original value computed from the
data) in the sequence

$\displaystyle(\mathrm{RdCov}_{n}^{(1)}(\bm{X};\bm{C}),\mathrm{RdCov}_{n}^{(2)}(\bm{X};\bm{C}),\ldots,\mathrm{RdCov}_{n}^{(B)}(\bm{X};\bm{C})),$

by breaking ties at random and where $R=1$ denotes the rank of the largest
element. Then the permutation $p$-value as $p_{B}:=R/(B+1)$. Note that if
$H_{0}$ is true,
$\mathbb{P}(p_{B}\leq\alpha)=\lfloor\alpha(B+1)\rfloor/(B+1)\leq\alpha$. In
other words,

$\displaystyle\tilde{\phi}_{n}(\bm{C}):=\bm{1}\\{p_{B}\leq\alpha\\}.$ (3.13)

is a finite sample level $\alpha$ test.

######  Remark 3.10.

Note that the resampling method described above is different from the (data-
dependent) bootstrap/permutation method used to calibrate
$\mathrm{JdCov}(\bm{X},\bm{C})$, or other independence tests which are not
distribution-free), where the asymptotic distribution depends on the (unknown)
marginal distributions of $X_{1},X_{2},\ldots,X_{r}$ (see [12, Proposition
9]). On the other hand, our method can be used to compute the cut-offs by only
permuting the Halton sequences, without any knowledge of the data, and hence
is much more efficient when one has to do multiple tests.

Although (3.13) produces a test which has level $\alpha$ in finite samples,
its power properties with a finite number of permutations is a-priori unclear.
This is due to fact that for any collection of $r$ random permutations
$\bm{\pi}=\\{\pi_{1},\pi_{2},\ldots,\pi_{r}\\}$, the RHS of (3.11) tends to
zero with high probability (see Proposition D.1 in Appendix D). This combined
with Theorem 2.12 leads to the following result:

######  Theorem 3.11 (Consistency with finite number of resamples).

Suppose Assumption 1 holds. Fix $\alpha\in(0,1)$ and suppose $B\geq
1/\alpha-1$. Then, for any fixed
$\mu\neq\mu_{1}\otimes\mu_{2}\otimes\cdots\otimes\mu_{r}$,

$\lim_{n\rightarrow\infty}\mathbb{E}_{\mu}[\tilde{\phi}(\bm{C})]=1.$

###### Proof.

Define $T:=\mathrm{RdCov}^{2}(\bm{X},\bm{C})$,
$T_{n}:=\mathrm{RdCov}_{n}(\bm{X};\bm{C})$, and
$T_{n}^{(b)}:=\mathrm{RdCov}_{n}^{(b)}(\bm{X};\bm{C})$, for $1\leq b\leq B$
(as defined in (3.11)). Then, since $B\geq 1/\alpha-1$,

$\displaystyle\mathbb{E}_{\mu}[\tilde{\phi}_{n}(\bm{C})]$
$\displaystyle=\mathbb{P}(p_{B}\leq\alpha)$
$\displaystyle\geq\mathbb{P}\left(T_{n}>T/2\text{ and }T_{n}^{(b)}<T/2\text{
for }1\leq b\leq B\right)$ $\displaystyle\geq
1-\left\\{\mathbb{P}(T_{n}<T/2)+B\cdot\mathbb{P}(T_{n}^{(1)}>T/2)\right\\}$
(by the union bound) $\displaystyle\rightarrow 1,$

where the last convergence is due to the consistency result in Theorem 2.12
and Proposition D.1 in Appendix D. ∎

##  4\. Local Power Analysis

In this section, we derive the asymptotic local power of the statistic
$\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})$ under the contiguous alternatives. In
the context of testing pairwise independence two types of local alternatives
have been considered in the literature: (1) mixture alternatives [20, 68] and
(2) Konijn alternatives which was originally proposed in [44] and has been
later investigated in [26] and used more recently in [20, 68]. To the best of
our knowledge, local power analysis has not been carried out for joint
independence testing of more than 2 variables. The first step towards this is
to define contiguous versions of the mixture and Konijn alternatives for
multiple random vectors. For this suppose the joint distribution $\mu$ has
density $f$ with respect to the Lebesgue measure in $\mathbb{R}^{d_{0}}$ and
the marginal distributions $\mu_{1},\mu_{2},\ldots,\mu_{r}$ have densities
$f_{X_{1}},f_{X_{2}},\ldots,f_{X_{r}}$ with respect to the Lebesgue measures
in $\mathbb{R}^{d_{1}},\mathbb{R}^{d_{2}},\ldots,\mathbb{R}^{d_{r}}$,
respectively.

###  4.1. Mixture Alternatives

The mixture alternatives constructed as follows: Given a $\delta>0$ consider
the following mixture density in $\mathbb{R}^{d_{0}}$:

$\displaystyle f(\bm{x})=(1-\delta)\prod_{i=1}^{r}f_{X_{i}}(x_{i})+\delta
g(\bm{x}),$

where $x_{i}\in\mathbb{R}^{d_{i}}$ for $1\leq i\leq r$,
$\bm{x}=(x_{1},x_{2},\ldots,x_{r})\in\mathbb{R}^{d_{0}}$, and
$g\neq\prod_{i=1}^{r}f_{i}$ is a probability density function with respect to
the Lebsegue measure in $\mathbb{R}^{d_{0}}$ such that the following holds:

######  Assumption 2.

The support of $g$ is contained in that of $\prod_{i=1}^{r}f_{i}$ and

$0<\mathbb{E}\left[\left(\frac{g(\bm{X})}{\prod_{i=1}^{r}f_{X_{i}}(X_{i})}-1\right)^{2}\right]<\infty,$

where the expectation is taken under $H_{0}$, that is,
$\bm{X}=(X_{1},X_{2},\ldots X_{r})\sim\prod_{i=1}^{r}f_{X_{i}}$.

Under this assumption contiguous local alternatives are obtained by
considering local perturbations of the mixing proportion $\delta$ as follows:

$H_{0}:\delta=0\qquad\mathrm{versus}\qquad H_{1}:\delta=h/\sqrt{n},$ (4.1)

for some $h>0$. This type of alternative captures all local additive
perturbations from $H_{0}$ by compressing potential lower-order dependence
among the variables in the functions $g$. The following theorem derives the
distribution of $\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})$ under $H_{1}$ as in
(4.1).

######  Theorem 4.1.

Suppose Assumptions 1 and 2 hold. Then under $H_{1}$ as in (4.1), as
$n\rightarrow\infty$,

$\displaystyle\\{n\mathrm{RdCov}_{n}(\bm{X}_{S}):S\in\mathcal{T}\\}\stackrel{{\scriptstyle
D}}{{\to}}\left\\{\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})+\mu_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}\right\\}_{S\in\mathcal{T}},$ (4.2)

with $Z_{S}(\bm{t})$ as defined in (3.6) and

$\displaystyle\mu_{S}(\bm{t})=h\cdot\mathbb{E}_{\bm{X}\sim
H_{0}}\left[\left(\frac{g(\bm{X})}{\prod_{i=1}^{r}f_{X_{i}}(X_{i})}-1\right)\prod_{i\in
S}\left\\{\mathbb{E}\left[e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right]-e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right\\}\right],$ (4.3)

where $\bm{t}=(t_{i})_{i\in S}$ and $t_{i}\in\mathbb{R}^{d_{i}}$ for $i\in S$.
Consequently, for any vector of non-negative weights
$\bm{C}=(C_{1},C_{2},\ldots,C_{r})$, under $H_{1}$ as in (4.1),

$\displaystyle n\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\stackrel{{\scriptstyle
D}}{{\to}}\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})+\mu_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}.$ (4.4)

The proof of this theorem is given in Appendix E. Using this result we can
derive the limiting local power of the test based on
$\mathrm{RdCov}^{2}(\bm{X};\bm{C})$. Specifically, if $\mathcal{H}_{\bm{C}}$
denotes the CDF of the limiting distribution in (4.4) then local asymptotic
power of the test $\phi_{n}^{\mathrm{asymp}}(\bm{C})$ (recall (3.10)) is given
by

$\lim_{n\rightarrow\infty}\mathbb{E}_{H_{1}}[\phi_{n}^{\mathrm{asymp}}(\bm{C})]=1-\mathcal{H}_{\bm{C}}^{(h)}(c_{\alpha}).$

This implies, $\phi_{n}^{\mathrm{asymp}}(\bm{C})$ (and similarly
$\phi_{n}(\bm{C})$ and $\tilde{\phi}_{n}(\bm{C})$ in (3.2) and (3.13),
respectively) have non-trivial Pitman efficiency and is rate-optimal, in the
sense that,

$\lim_{|h|\rightarrow\infty}\lim_{n\rightarrow\infty}\mathbb{E}_{H_{1}}[\phi_{n}^{\mathrm{asymp}}(\bm{C})]=1.$

###  4.2. Konijn Alternative

The Konijn family of alternatives [44] for 2 random vectors can be defined as
follows: Given two independent random vectors
$X_{1}^{\prime}\sim\mu_{1}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{1}})$,
$X_{2}^{\prime}\sim\mu_{1}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{2}})$, and
$\delta>0$, consider the law of

$\displaystyle\begin{pmatrix}X_{1}\\\
X_{2}\end{pmatrix}=\begin{pmatrix}(1-\delta)\bm{I}_{d_{1}}&\delta\bm{M}_{1}\\\
\delta\bm{M}_{2}&(1-\delta)\bm{I}_{d_{2}}\end{pmatrix}\begin{pmatrix}X_{1}^{\prime}\\\
X_{2}^{\prime}\end{pmatrix}$ (4.5)

where $\bm{M}_{1}$ and $\bm{M}_{2}$ are $d_{1}\times d_{2}$ and $d_{2}\times
d_{1}$ dimensional deterministic matrices, respectively. Note that when
$\delta=0$, the random vectors $X_{1},X_{2}$ are independent and as $\delta$
increases introduces dependence/correlation through the matrix $\bm{M}$.
Gieser [25] showed that if one choses $\delta=h/\sqrt{n}$, for some $h>0$,
then the distributions of $(X_{1},X_{2})$ and
$(X_{1}^{\prime},X_{2}^{\prime})$ are contiguous when $X_{1}^{\prime}$ and
$X_{2}^{\prime}$ are elliptically symmetric distributions centered at
$\theta_{1}\in\mathbb{R}^{d_{1}}$ and $\theta_{2}\in\mathbb{R}^{d_{2}}$ with
covariance matrices $\Sigma_{1}$ and $\Sigma_{2}$, respectively. Here, we
extend the Konijn family to multiple random vectors as follows:

######  Definition 4.2.

Given $r$ independent random vectors
$X_{1}^{\prime},X_{2}^{\prime},\ldots,X_{r}^{\prime}$, where
$X_{i}^{\prime}\sim\mu_{i}\in\mathcal{P}_{ac}(\mathbb{R}^{d_{i}})$, for $1\leq
i\leq r$, and $\delta>0$, consider the law of

$\displaystyle\bm{X}=\begin{pmatrix}X_{1}\\\ X_{2}\\\ \vdots\\\
X_{r}\end{pmatrix}=\begin{pmatrix}(1-\delta)\bm{I}_{d_{1}}&\delta\bm{M}_{1,2}&\cdots&\delta\bm{M}_{1,r}\\\
\delta\bm{M}_{2,1}&(1-\delta)\bm{I}_{d_{2}}&\cdots&\delta\bm{M}_{2,r}\\\
\vdots&\vdots&\ddots&\vdots\\\
\delta\bm{M}_{r,1}^{\top}&\delta\bm{M}_{r,2}&\cdots&(1-\delta)\bm{I}_{d_{r}}\\\
\end{pmatrix}\begin{pmatrix}X_{1}^{\prime}\\\ X_{2}^{\prime}\\\ \vdots\\\
X_{r}^{\prime}\end{pmatrix}:=\bm{A}_{\delta}\begin{pmatrix}X_{1}^{\prime}\\\
X_{2}^{\prime}\\\ \vdots\\\ X_{r}^{\prime}\end{pmatrix},$ (4.6)

where $\bm{M}_{i,j}$, for $1\leq i\neq j\leq r$, is a $d_{i}\times d_{j}$ a
dimensional deterministic matrix.

Note that the matrix $\bm{A}_{\delta}$ in (4.6) is a $d_{0}\times d_{0}$ times
block matrix with $\bm{I}_{d_{i}}$ in the $i$-th diagonal block for $1\leq
i\leq r$ and $\bm{M}_{i,j}$ in the $(i,j)$-th block for $1\leq i\neq j\leq r$.
When $\delta=0$, the matrix $\bm{A}_{\delta}$ is identity and, hence, by
continuity there exists a neighborhood $\Theta$ of zero such that
$\bm{A}_{\delta}$ is invertible for $\delta\in\Theta$. Therefore, by a change
of variable formula, for $\delta\in\Theta$, the density of $\bm{X}$ can be
written as:

$\displaystyle
f_{\bm{X}}(\bm{x}|\delta)=\frac{1}{|\mathrm{det}(\bm{A}_{\delta})|}f_{\bm{X}^{\prime}}(\bm{A}_{\delta}\bm{x}),$
(4.7)

where $f_{\bm{X}^{\prime}}(\cdot)=\prod_{i=1}^{r}f_{X_{i}^{\prime}}(\cdot)$ is
the density of $(X_{1}^{\prime},X_{2}^{\prime},\ldots,X_{r}^{\prime})$.
Hereafter, we will denote
$f^{\prime}_{\bm{X}}(\bm{x}|\delta)=\frac{\mathrm{d}}{\mathrm{d}\theta}f_{\bm{X}}(\bm{x}|\bm{\theta})|_{\theta=\delta}$
and assume the following:

######  Assumption 3.

The family of distributions $\\{f_{\bm{X}}(\cdot|\delta)\\}_{\delta\in\Theta}$
satisfies the following:

  * •

The support of $f_{\bm{X}}(\cdot|\delta)$ does not depend on
$\delta\in\Theta$. In other words, the set
$\mathcal{S}:=\\{\bm{x}\in\mathbb{R}^{d_{0}}:f_{\bm{X}}(\bm{x}|\delta)>0\\}$
does not depend on $\delta$.

  * •

The map $\bm{x}\mapsto\sqrt{f_{\bm{X}}(\bm{x}|0)}$ is continuously
differentiable.

  * •

The Fisher information

$I(\delta)=\int_{\mathbb{R}^{d_{0}}}\left(\frac{f^{\prime}_{\bm{X}}(\bm{x}|\delta)}{f_{\bm{X}}(\bm{x}|\delta)}\right)^{2}f_{\bm{X}^{\prime}}(x)\mathrm{d}x$

is well-define, finite, strictly positive, and continuous at $\delta=0$.

Under the above assumption, we consider the following sequence of local
hypotheses:

$H_{0}:\delta=0\qquad\mathrm{versus}\qquad H_{1}:\delta=h/\sqrt{n},$ (4.8)

for some $h>0$. We show in Lemma E.1 that $H_{0}$ and $H_{1}$ as in (4.8) are
contiguous under Assumption 3. This assumption also ensures that the family
$\\{f_{\bm{X}}(\cdot|\delta)\\}_{\delta\in\Theta}$ is quadratic mean-
differentiable (QMD) (see [47, Definition 12.2.1]) and, hence, Le Cam’s theory
of contiguity applies (see [47, Chapter 12] and [71, Chapter 6] ) .

######  Remark 4.3.

Assumption 3 is satisfied, for example, when $X_{1},X_{2},\ldots,X_{r}$ are
centered elliptical distributions, under certain non-degeneracy conditions on
the respective scale matrices (as discussed [68, Example 5.1] for the 2
variables case).

To describe the limiting distribution of $\mathrm{RJdCov}(\bm{X};\bm{C})$
under $H_{1}$, denote the likelihood ratio and its derivative as:

$L(\bm{x}|\delta):=\frac{f_{\bm{X}}(\bm{x}|\delta)}{f_{\bm{X}}(\bm{x}|0)}=\frac{f_{\bm{X}}(\bm{x}|\delta)}{f_{\bm{X}^{\prime}}(\bm{x})}\quad\text{and}\quad
L^{\prime}(\bm{x}|\delta):=\frac{\nabla
f_{\bm{X}}(\bm{x}|\delta)}{f_{\bm{X}^{\prime}}(\bm{x})},$

where $\nabla
f_{\bm{X}}(\bm{x}|\delta)=\frac{\partial}{\partial\delta}f_{\bm{X}}(\bm{x}|\delta)$.

######  Theorem 4.4 (Local power for Konijn alternatives).

Suppose Assumptions 1 and 3 hold. Then under $H_{1}$ as in (4.8), as
$n\rightarrow\infty$,

$\displaystyle\\{n\mathrm{RdCov}_{n}(\bm{X}_{S}):S\in\mathcal{T}\\}\stackrel{{\scriptstyle
D}}{{\to}}\left\\{\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})+\nu_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}\right\\}_{S\in\mathcal{T}},$ (4.9)

with $Z_{S}(\bm{t})$ as defined in (3.6) and

$\displaystyle\nu_{S}(\bm{t})=h\cdot\mathbb{E}_{\bm{X}\sim
H_{0}}\left[L^{\prime}(\bm{X}|0)\prod_{i\in
S}\left\\{\mathbb{E}\left[e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right]-e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i})\rangle}\right\\}\right],$ (4.10)

where $\bm{t}=(t_{i})_{i\in S}$ and $t_{i}\in\mathbb{R}^{d_{i}}$ for $i\in S$.
Consequently, for any vector of non-negative weights
$\bm{C}=(C_{1},C_{2},\ldots,C_{r})$, under $H_{1}$ as in (4.1),

$\displaystyle n\mathrm{RJdCov}^{2}_{n}(\bm{X};\bm{C})\stackrel{{\scriptstyle
D}}{{\to}}\sum_{s=2}^{r}C_{s}\sum_{\begin{subarray}{c}S\subseteq\\{1,2,\ldots,r\\}\\\
|S|=s\end{subarray}}\int_{\mathbb{R}^{d_{S}}}|Z_{S}(\bm{t})+\nu_{S}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}.$ (4.11)

###  4.3. Hájek Representation and Proof Outline

The proofs of Theorems 4.1 and 4.4 are given in Appendix E. One of the main
ingredients of the proof is a Hájek representation for the empirical process
associated with $\mathrm{RdCov}_{n}(\bm{X}_{S})$, which shows that we can
replace the empirical rank maps in the definition of
$\mathrm{RdCov}_{n}(\bm{X}_{S})$ with their population counterparts without
altering its limiting distribution under $H_{0}$. Since this result could be
of independent interest we summarized this as a proposition below. To describe
the result, consider the following empirical process:

$\displaystyle Z_{S,n}(\bm{t})$
$\displaystyle=\frac{1}{n}\sum_{b=1}^{n}\left\\{\prod_{i\in
S}\left(\frac{1}{n}\sum_{a=1}^{n}e^{\iota\langle
t_{i},\hat{R}_{i}(X_{i}^{(a)})\rangle}-e^{\iota\langle
t_{i},\hat{R}_{i}(X_{i}^{(b)})\rangle}\right)\right\\},$ (4.12)

where $S\in\mathcal{T}$, $\bm{t}=(t_{i})_{i\in S}\in\mathbb{R}^{d_{S}}$ and
$t_{i}\in\mathbb{R}^{d_{i}}$, for $i\in S$. Recalling (2.14), note that

$\displaystyle\mathrm{RdCov}^{2}_{n}(\bm{X}_{S})=\int|Z_{S,n}(\bm{t})|^{2}\prod_{i\in
S}\mathrm{d}w_{i}.$ (4.13)

One of the challenges of dealing with the process $Z_{S,n}(\bm{t})$ is the
dependence across the indices $i\in S$ arising from the empirical rank maps
$\\{\hat{R}_{i}(\cdot)\\}_{i\in S}$. The Hájek representation result shows
that this dependence is asymptotically negligible under $H_{0}$. Towards this,
define the oracle version of the empirical process $Z_{S,n}$ where the
empirical rank maps $\\{\hat{R}_{i}(\cdot)\\}_{i\in S}$ in replaced with their
population analogues empirical rank maps $\\{R_{\mu_{i}}(\cdot)\\}_{i\in S}$

$\displaystyle Z_{S,n}^{\textrm{oracle}}(\bm{t})$
$\displaystyle=\frac{1}{n}\sum_{b=1}^{n}\left\\{\prod_{i\in
S}\left(\frac{1}{n}\sum_{a=1}^{n}e^{\iota\langle
t_{i},R_{\mu_{i}}(X_{i}^{(a)})\rangle}-e^{\iota\langle
t_{i}R_{\mu_{i}}(X_{i}^{(b)})\rangle}\right)\right\\},$ (4.14)

for $\bm{t}$ as in (3.6). We are now ready the state the Hájek representation
result which shows that the difference between $Z_{S,n}$ and its oracle
version is $o_{P}(1/\sqrt{n})$ under $H_{0}$.

######  Proposition 4.5.

Suppose Assumption 1 holds. Then under $H_{0}$ as $n\rightarrow\infty$,

$\sqrt{n}\left(Z_{S,n}(\bm{t})-Z_{S,n}^{\mathrm{oracle}}(\bm{t})\right)=o_{P}(1),$

for all $S\in\mathcal{T}$ and $\bm{t}\in\mathbb{R}^{d_{S}}$.

Once the Hájek representation is established the limiting distribution of
$\\{\mathrm{RdCov}(\bm{X}_{S})\\}_{S\in\mathcal{T}}$ under contiguous local
alternatives can be obtained as follows:

  * •

The first step is to show the alternatives considered in (4.1) and (4.8) are
mutually contiguous with $H_{0}$. While this is well-known for the case of 2
variables and relatively straightforward to extend to multiple variables for
the mixture alternative (4.1), for the Konjin alternative (4.8) establishing
contiguity is much more involved. Invoking Le Cam’s second lemma [71], we show
in Lemma E.1 that Konjin alternatives as in (4.6) are contiguous under
Assumption 3.

  * •

Next, we derive the joint distribution of the collection of empirical
processes $\\{\sqrt{n}Z_{S,n}^{\mathrm{oracle}}(\bm{t})\\}_{S\in\mathcal{T}}$
and the log-likelihood under $H_{0}$. Proposition 4.5 together with contiguity
and Le Cam’s third lemma [71] gives the joint distribution of
$\\{\sqrt{n}Z_{S,n}^{\mathrm{oracle}}(\bm{t})\\}_{S\in\mathcal{T}}$ under
$H_{1}$ as in (4.1) and (4.8). The results in (4.2) and (4.9) then follow by
the representation in (4.13) and the continuous mapping theorem.

######  Remark 4.6.

Note that, as expected, one obtains the limiting null distribution of
$n\mathrm{RJdCov}_{n}(\bm{X};\bm{C})$ by substituting $h=0$ in Theorems 4.1
and 4.4. In fact the Hájek representation result provides an alternative
strategy to proving of the asymptotic null distribution, which does require
the combinatorial CLT for multiple permutations. Nevertheless, we present both
approaches because the technical by-products emerging from them, namely the
Hájek representation and the combinatorial CLT, are results of independent
interest. In particular, the scope of the combinatorial CLT with multiple
permutations goes beyond rank-based methods to other nonparametric
independence tests, such as methods based on geometric graphs [22].

To validate the asymptotic results discussed above we present a finite sample
simulation comparing the empirical power of the $\mathrm{RJdCov}$ based test
with test based on $\mathrm{JdCov}$ studied in [12]. Towards this we consider
the following 3 distributions: (1) the multivariate normal (MVN), (2) the
third-order multivariate normal copula (MVN copula),444A random variable $X$
is said to follow a $p$-dimensional $r$-th order, for $r>0$, multivariate
normal copula if $X\stackrel{{\scriptstyle D}}{{=}}Z^{r}$, where $Z\sim
N_{p}(\bm{\mu},\bm{\Sigma})$ as in [12, Example 1], and (3) the multivariate
student distribution (MVT) with 5 degrees of freedom, with dimensions
$d_{1}=d_{2}=d_{3}=2$. In all the 3 cases the mean vectors are set to zero and
the covariance matrix $\Sigma=((\sigma_{ij}))_{1\leq i,j\leq 2}$ is
$\sigma_{11}=\sigma_{22}=1$ and $\sigma_{12}=\sigma_{21}=0.5$. We set the
number of variables $r=3$, the matrices $\bm{M}_{i,j}=\bm{I}_{2}$, for $1\leq
i\neq j\leq 3$ and compute the matrix $\bm{A}_{\delta}$ as in (4.6) by varying
$\delta$ as follows: $\delta\in\\{0,0.4,0.8,1.2,1.6\\}$ for the Gaussian case,
$\delta\in\\{0,0.1,0.2,0.3,0.4\\}$ for the Gaussian copula case, and
$\delta\in\\{0,0.25,0.5,0.75,1\\}$ for the multivariate $t$-distribution.
Table 1 shows the empirical power (over 200 repetitions) with sample size
$n=300$ as a function of the signal strength ($\delta$ varying as above), for
tests based on the $\mathrm{JdCov}$ and the $\mathrm{RJdCov}$. The plots show
that our method has comparable power with $\mathrm{JdCov}$ for the Gaussian
and $t$-distributions, while it is significantly better than the
$\mathrm{JdCov}$ in the more heavy-tailed Gaussian copula model, illustrating
the attractive efficiency properties of rank-based distribution-free methods.

(a)

(b)

(c)

Figure 1. Empirical power of against Konijn alternatives with 3 variables: (a)
bivariate normal, (b) bivariate normal copula model of order $3$, and (c)
bivariate $t$-distribution with $5$ degrees of freedom.

##  5\. Robust Independent Component Analysis

In this section, we apply the $\mathrm{RJdCov}$ measure to develop a robust
method for independent component analysis (ICA). Towards this, suppose
$X\in\mathbb{R}^{r}$ be a random vector satisfying the following the moment
conditions:

######  Assumption 4.

The random vector $X$ has a continuous distribution with $\mathbb{E}[X]=0$ and
$\mathbb{E}[\|X\|^{2}]<\infty$.

The ICA model assume that there is a random vector $S\in\mathbb{R}^{r}$ with
jointly independent coordinates such that $X=\bm{M}S$, where $\bm{M}$ is a
nonsingular mixing matrix. Without loss of generality, we assume $S$ is
standardized such that $\mathbb{E}[S]=0$ and
$\operatorname{Var}[S]=\bm{I}_{r}$. Then by the spectral decomposition:

$\Sigma_{X}=\operatorname{Var}[X]=\bm{P}\Lambda\bm{P}^{\top},$

where $\bm{P}$ and $\Lambda$ are the $r\times r$ matrix of eigenvectors and
the diagonal matrix of the corresponding eigenvalues of $\Sigma_{\bm{X}}$,
respectively. Define $\bm{O}:=\Lambda^{-\frac{1}{2}}\bm{P}^{\top}$ and
$Z=\bm{O}X$. Note that $\mathbb{E}[Z]=0$ and
$\operatorname{Var}[Z]=\bm{I}_{r}$. Therefore, when the ICA model holds we can
write

$\displaystyle S=\bm{M}^{-1}X=\bm{M}^{-1}\bm{O}^{-1}Z=\bm{W}Z,$ (5.1)

where $\bm{W}=\bm{M}^{-1}\bm{O}^{-1}$. This implies,
$\bm{I}=\operatorname{Var}[S]=\bm{W}\operatorname{Var}[Z]\bm{W}^{\top}=\bm{W}\bm{W}^{\top}$,
that is, $\bm{W}$ is an orthogonal matrix.

Following Matteson and Tsay [52] we restrict $\bm{W}\in\mathcal{SO}(r)$, where
$\mathcal{SO}(r)$ denotes the set $r\times r$ orthogonal matrices with
determinant 1 (the rotation group) and parameterize $\bm{W}$ by a vector
$\theta$, where $\theta$ is a vectorized triangular array of rotation angles
of length $r(r-1)/2$ indexed by $\\{i,j:1\leq i<j\leq r\\}$, as follows:

$\displaystyle\bm{W}=\bm{W}(\theta)=\bm{Q}^{(r-1)}\cdots\bm{Q}^{(1)},\quad\text{with
}\bm{Q}^{(i)}=\bm{Q}_{i,r}(\theta_{i,r})\cdots\bm{Q}_{i,i+1}(\theta_{i,i+1}),$
(5.2)

where $\bm{Q}_{i,j}(\theta)$ is the identity matrix $\bm{I}_{r}$ with the
$(i,i)$ and $(j,j)$ elements replaced by $\cos(\theta)$, the $(i,j)$ element
replaced by $-\sin(\theta)$ and the $(j,i)$ element replaced by
$\sin(\theta)$. To ensure the existence and uniqueness of the representation
(5.2) define

$\displaystyle\Theta=\Bigg{\\{}\theta_{i,j}:\begin{cases}0\leq\theta_{1,j}<2\pi,\\\
0\leq\theta_{i,j}<\pi,&i\neq 1.\end{cases}\Bigg{\\}}.$ (5.3)

Then by [51, Theorem 2.3.2] there exists a unique inverse map of
$\bm{W}\in\mathcal{SO}(r)$ into $\theta\in\Theta$ such that the mapping is
assured to be continuous if either all elements on the main-diagonal of $W$
are positive, or all elements of $W$ are nonzero (see [52, 51] for more
details on such decompositions).

Another identification issue with the ICA model is the sign and order of the
independent components. Suppose $\bm{P}_{\pm}$ is a signed permutation matrix
and then note that $X=\bm{M}S=\bm{M}\bm{P}^{\top}_{\pm}\bm{P}_{\pm}S$ in which
case $\bm{P}_{\pm}S$ are the new independent components and
$\bm{M}\bm{P}^{\top}_{\pm}$ is the new mixing matrix. This necessitates
considering metrics that are invariant to the scale, sign, and order of the
independent components, when comparing different estimates (see Section 6.4
for details).

In general the ICA model (5.1) is misspecified, in which case $\theta$ is
chosen such that the components of $\bm{W}(\theta)Z$ are as close to mutually
independent as possible. To this end, the $\mathrm{RJdCov}$ provides a robust
measure for quantifying deviations from joint independence. Hence, one choose
$\theta$ by minimizing the following objective function:

$\displaystyle\theta_{0}={\arg\min}_{\theta}\mathrm{RJdCov}^{2}(\bm{W}(\theta)Z;c).$
(5.4)

Recall that the (population) rank map in the 1-dimensional case is the
distribution function, hence, $\mathrm{RJdCov}^{2}(\bm{W}(\theta)Z;c)$ can be
computed by replacing $R_{\mu_{i}}$ in (2.8) with $F_{W_{i}(\theta)^{\top}Z}$,
the distribution function of $F_{W_{i}(\theta)^{\top}Z}$, for $1\leq i\leq r$,
where $W_{i}(\theta)$ denotes the transpose of $i$-th row of $\bm{W}(\theta)$.

###  5.1. $\mathrm{RJdCov}$-Based ICA Estimator

We now discuss how one can estimate $\theta$ based on samples by minimizing
the empirical version of (5.4). Towards this, suppose
$\bm{X}_{n}=\\{X_{1},X_{2},\ldots X_{n}\\}$ be i.i.d. samples distributed as
$X$ satisfying Assumption 4. We define approximately uncorrelated observations
$\hat{\bm{Z}}_{n}=\bm{X}_{n}\hat{\bm{O}}^{\top}$, where $\hat{\bm{O}}_{n}$ is
the sample covariance matrix of $\bm{X}_{n}$. Assumption 4 implies that
$\operatorname{Var}[\hat{\bm{Z}}_{n}]\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\bm{I}_{r}$. To simplify notation, we omit the above
steps described and assume $\bm{Z}_{n}=\\{Z_{1},Z_{2},\ldots Z_{n}\\}$ are
i.i.d. samples distributed as $Z$, which has zero mean vector and identity
covariance matrix. Also, we arrange $\\{Z_{1},Z_{2},\ldots Z_{n}\\}$ such that
the $i$-th row of $\bm{Z}_{n}$ is $Z_{i}^{\top}$, that is, $\bm{Z}_{n}$ has
dimension $n\times r$.

To estimate $\theta$ from $\bm{Z}_{n}$ we replace $\mathrm{RJdCov}^{2}$ in
(5.4) with its empirical version

$\displaystyle\hat{\theta}_{n}={\arg\min}_{\theta}\mathrm{RJdCov}^{2}_{n}(\bm{Z}_{n}\bm{W}(\theta)^{\top};c)$
(5.5)

Note that for any value of $\theta$ the objective function
$\mathrm{RJdCov}^{2}_{n}(\bm{Z}_{n}\bm{W}(\theta)^{\top};c)$ can be computed
easily, since in dimension 1 the empirical rank map is the empirical
distribution function. Specifically, the empirical rank map for the $i$-th
component is
$\hat{F}_{i,n}(t|\theta)=\frac{1}{n}\sum_{a=1}^{n}\bm{1}\\{W_{i}^{\top}(\theta)Z_{a}\leq
t\\}$. To establish the consistency of $\hat{\theta}_{n}$ we need the
following definition:

######  Definition 5.1.

We say $\hat{\theta}_{n}$ is consistent for $\theta$ over $\mathcal{SO}(r)$ if
the function $\bm{W}(\theta)$ is continuous in $\theta$ and
$\mathcal{D}(\bm{W}(\hat{\theta}_{n}),\bm{W}(\theta))\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$, as $n\rightarrow\infty$, for any metric $\mathcal{D}$
on $\mathcal{SO}(r)$ satisfying $\mathcal{D}(\bm{W},\bm{A})=0$ if and only if
there exists a signed permutation matrix $\bm{P}_{\pm}$ such that
$\bm{W}=\bm{P}_{\pm}\bm{A}$.

With this definition we can now state the consistency result of
$\hat{\theta}_{n}$. Towards this, recall the definition of $\Theta$ from (5.3)
and let $\bar{\Theta}$ denote a sufficiently large compact subset of the space
$\Theta$.

######  Theorem 5.2.

Suppose Assumption 4 holds and the density function of $Z$ is uniformly
bounded. Moreover, assume there exists a unique minimizer
$\theta_{0}\in\bar{\Theta}$ for (5.4) and $\bm{W}(\theta)$ satisfies the
conditions for a unique continuous inverse to exist, then $\hat{\theta}_{n}$
is consistent for $\theta_{0}$ over $\mathcal{SO}(d)$.

The proof of Theorem 5.2 is given in Appendix F.2. This shows that the
proposed estimator $\hat{\theta}_{n}$ is consistent for $\theta_{0}$ for
metrics which are invariant to the non-identifiability issues of the ICA
problem. The assumption of uniform boundedness is required to guarantee the
distribution function is Lipschitz, which ensures the equicontinuity of the
objective function.

###  5.2. Practical Implementation

One of the challenges in directly implementing gradient based methods to
minimize the objective in (5.5) is that, the empirical distribution functions
are not differentiable with respect to $\theta$. Matteson and Tsay [52]
suggested a practical way to circumvent this issue by estimating the
univariate distribution functions $F_{W_{i}(\theta)^{\top}Z}$, for $1\leq
i\leq r$, using kernel-based methods. Specifically, we consider the following
problem:

$\displaystyle\tilde{\theta}_{n}=\arg\min_{\theta}\mathrm{JdCov}^{2}_{n}((\tilde{F}_{1}(\bm{Z}_{n}W_{1}(\theta)),\ldots,\tilde{F}_{r}(\bm{Z}_{n}W_{r}(\theta));c),$
(5.6)

with
$\tilde{F}_{1}(\bm{Z}_{n}W_{1}(\theta))=(\check{F}_{1}(W_{1}(\theta)^{\top}Z_{1}),\ldots,\tilde{F}_{1}(W_{1}(\theta)^{\top}Z_{n}))^{\top}$
and

$\displaystyle\tilde{F}_{i}(s)=\frac{1}{n}\sum_{a=1}^{n}G\bigg{(}\frac{s-W_{i}(\theta)^{\top}Z_{a}}{h_{n}(i)}\bigg{)},$
(5.7)

where $G$ is the integral of a density kernel and $h_{n}(i)$ is a data
dependent bandwidth for the $i$-th component, for $1\leq i\leq r$ (see Section
2 in [52] for practical choices of $h_{n}(i)$). The optimization problem (5.6)
can now be solved using a gradient-based method (see Section F.1 for the
computation of the gradient).

The advantages of using kernel-based estimates is that they have nice
theoretical properties, which can be used to establish the consistency of the
estimate $\tilde{\theta}_{n}$. In particular, if we assume that $h_{n}(i)$ is
a positive measurable function of the $i$-th coordinates of
$Z_{1},Z_{2},\ldots,Z_{n}$ such that $h_{n}(i)\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$ as $n\rightarrow\infty$, then

$\displaystyle\sup_{s\in\mathbb{R}}\big{|}\tilde{F}_{i}(s)-F_{i}(s)\big{|}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0,$

as $n\rightarrow\infty$, for $1\leq i\leq r$ (see Corollary 1 in [11]). Then
whenever $G$ is Lipschitz and Assumption 4 and the conditions in Theorem 5.2
holds, similar to the proof Theorem 2.1 in [52], by replacing their objective
function with $\mathrm{RJdCov}_{n}^{2}$, it can be shown that
$\tilde{\theta}_{n}$ is consistent for $\theta_{0}$ over $\mathcal{SO}(r)$.

##  6\. Simulation Study

In this section, we showcase the efficacy of the proposed test, both in terms
of Type-I error control and power, across a range of simulation settings. For
our test we implement the finite sample version discussed in Section 3.2 with
cut-off chosen as in (3.12) and all the weights set to 1. We compare our
method (hereafter referred to as $\mathrm{RJdCov}$) with the JdCov based test
in [12], the dCov based test studied by Matteson and Tsay [52] (hereafter
denoted as MT), and the dHSIC based test [59]. The section is organized as
follows: In Section 6.1 we present the results on Type-I error control. The
empirical power for testing joint independence is given in Section 6.2 and
results for detecting higher-order dependence are in Section 6.3. The
performance of the ICA estimator is studied in Section 6.4.

###  6.1. Type-I Error Control

To illustrate the finite sample Type-I error control of the proposed test we
consider $\bm{X}=(X_{1}^{\top},X_{2}^{\top},X_{3}^{\top})^{\top}$ in the
following three settings:

  * $(1)$

Multivariate Gaussian: $X_{1},X_{2},X_{3}$ are i.i.d.
$N_{3}(\bm{0},\bm{I}_{3})$.

  * $(2)$

Third-order Gaussian copula: $X_{1},X_{2},X_{3}$ are independent random
vectors in $\mathbb{R}^{3}$ with i.i.d. coordinates distributed as
$N(0,1)^{3}$.

  * $(3)$

Cauchy: $X_{1},X_{2},X_{3}$ are independent random vectors in $\mathbb{R}^{3}$
with i.i.d. coordinates distributed as $\mathcal{C}(0,1)$, the Cauchy
distribution with location parameter 0 and scale to parameter 1.

Table 1 shows the empirical Type-I error probability (over 500 repetitions)
with sample size $n=500$ for the different tests in the above 3 situations.
Our test satisfactorily controls Type-I error in finite samples, validating
the theoretical results. The other tests also have good Type-I error control.

Table 1. Empirical probability of Type-I error of the different tests.  Distribution | Multivariate Gaussian | Gaussian Copula | Cauchy
---|---|---|---
Proposed | 0.038 | 0.042 | 0.042
JdCov | 0.036 | 0.030 | 0.046
Matteson-Tsay (MT) | 0.056 | 0.046 | 0.056
dHSIC | 0.032 | 0.074 | 0.050

###  6.2. Power Results

In this section we compare the empirical power of $\mathrm{RJdCov}$ with the
existing methods in the 4 settings described below. In the first 2 cases we
have (light-tailed) multivariate Gaussian distributions with covariance
structure as in [12, Example 2], the third setting is a (heavy-tailed) Cauchy
regression model, and the fourth is a Cauchy sine dependence model (similar to
[78]). Throughout we set the sample size $n=500$ and compute the empirical
power over 500 repititions.

(a)

(b)

Figure 2. Power of the different tests in (a) the multivariate Gaussian model
with autoregressive correlation structure and (b) the multivariate Gaussian
model with banded covariance matrix.

  * $(1)$

Multivariate Gaussian with autoregressive covariance structure: In this case
the data is generated as follows: $\bm{X}\sim N_{9}(\bm{0},\Sigma)$ where
$\Sigma=((\sigma_{ij}))_{1\leq i,j\leq 9}$ has an autocorrelation structure,
that is, $\sigma_{ij}=\rho^{|i-j|}$. Then the first 3 coordinates of $\bm{X}$
are set as $X_{1}$, second 3 coordinates as $X_{2}$, and the last 3
coordinates as $X_{3}$. Figure 2 (a) shows the empirical power of the
different tests as $\rho$ varies between $0.05$ to $0.25$ over a grid of 5
equally spaced values.

  * $(2)$

Multivariate Gaussian with banded covariance matrix: Here, $\bm{X}\sim
N_{9}(\bm{0},\Sigma)$ where $\Sigma=((\sigma_{ij}))_{1\leq i,j\leq 9}$ has a
banded structure, that is, $\sigma_{ii}=1$ for $1\leq i\leq 9$,
$\sigma_{ij}=\rho$ for $1\leq|i-j|\leq 2$, and $\sigma_{ij}=0$ otherwise. As
before, the first 3 coordinates of $\bm{X}$ are set as $X_{1}$, second 3
coordinates as $X_{2}$, and the last 3 coordinates as $X_{3}$. Figure 2 (b)
shows the empirical power of the different tests as $\rho$ varies between
$0.05$ to $0.25$ over a grid of 5 equally spaced values.

  * $(3)$

Cauchy regression: In this case
$\bm{X}=(X_{1}^{\top},X_{2}^{\top},X_{3}^{\top})^{\top}$ is generated as
follows:

$\displaystyle X_{i}=Z_{i}+a\cdot V,$ (6.1)

for $1\leq i\leq 3$, where $Z_{i}\in\mathbb{R}^{3}$ have i.i.d. coordinates
distributed according to $\mathcal{C}(0,1)$, and
$V=(V_{1},V_{2},V_{3})^{\top}$ with $V_{1}=V_{2}=V_{3}\sim\mathcal{C}(0,1)$.
Note that the variables $X_{1},X_{2},X_{3}$ are mutually dependent because of
the shared noise vector. Figure 3 (a) shows the empirical power of the
different tests as $a$ varies between $0.2$ to 1 over a grid of 5 equally
spaced values.

  * $(4)$

Sine dependence: In this case
$\bm{X}=(X_{1}^{\top},X_{2}^{\top},X_{3}^{\top})^{\top}$ is generated as in
(6.1), but with a different noise distribution. Specifically, we fix $a=1$,
generate $Z_{i}\in\mathbb{R}^{3}$ with i.i.d. coordinates distributed
according to $\mathcal{C}(0,1)$ as before, but change
$V=(V_{1},V_{2},V_{3})^{\top}$ as $V_{1}=V_{2}=V_{3}=\sin(b\cdot W)$ where
$W\sim\mathcal{C}(0,1)$. Figure 3 (b) shows the empirical power of the
different tests as $b$ varies between $0.1$ to $0.5$ over a grid of 5 equally
spaced values.

We observe from the plots in Figure 2 that the $\mathrm{RJdCov}$ along with
the other methods have very similar performance in the Gaussian settings. On
the other hand, Figure 3 shows in the heavy-tailed Cauchy and sine dependence
settings the $\mathrm{RJdCov}$ significantly outperforms the other tests.
Overall the $\mathrm{RJdCov}$ emerges as the preferred method by exhibiting
good power over a range of data distributions.

(a)

(b)

Figure 3. Power of the different tests in (a) the Cauchy regression model and
(b) the sine dependence model.

###  6.3. Higher-Order Dependence

In this section we consider situations where the variables are pairwise
independent but jointly dependent. Specifically, generate $(X,Y,Z^{\prime})$
as follows:

  * •

Suppose $X$, $Y$ and $Z$ are independent random vectors with i.i.d.
coordinates distributed according to $F$, where $F$ is a distribution in
$\mathbb{R}$ that is symmetric around $0$.

  * •

For $1\leq s\leq d-1$, set $Z_{s}^{\prime}=Z_{s}$ and set
$Z^{\prime}_{d}=-Z_{d}$ if $X_{d}Y_{d}Z_{d}$ is positive, otherwise set
$Z^{\prime}_{d}=Z_{d}$.

It can easily be checked that $X,Y,Z^{\prime}$ are pairwise independent.
However, they are mutually dependent, since, for example,
$Z_{d}^{\prime}X_{d}Y_{d}$ is always non-negative. In the simulations we
choose $F$ to the following 4 distributions: (1) $N(0,1)$, (2) the
$t$-distribution with 2 degrees of freedom $(t_{(2)})$, (3) the
$t$-distribution with 3 degrees of freedom $(t_{(3)})$, and (4) the Cauchy
distribution $\mathcal{C}(0,1)$, equivalently, the the $t$-distribution with 1
degree of freedom.

In Table 2 we show the empirical power (over 500 iterations) with sample size
$n=500$ of the following 3 tests (and their corresponding versions based on
$\mathrm{dCov}$/$\mathrm{JdCov}$):

  * •

The test for pairwise independence that rejects for large values of

$\mathrm{RdCov}^{2}_{n}(X,Y)+\mathrm{RdCov}^{2}_{n}(Y,Z^{\prime})+\mathrm{RdCov}^{2}_{n}(X,Z^{\prime}).$

  * •

The test for higher-order independence that rejects for large values of
$\mathrm{RdCov}^{2}_{n}(X,Y,Z^{\prime})$.

  * •

The test for joint independence based on $\mathrm{RJdCov}_{n}^{2}$.

The results in Table 2 show that for all the 4 distributions considered both
the $\mathrm{RJdCov}$ and $\mathrm{JdCov}$ based tests successfully detect the
null hypothesis of pairwise independence (probability of Type-I error).
However, the $\mathrm{JdCov}$ based test fails to detect higher-order
dependence for the Cauchy distribution and joint dependence for the $t_{(2)}$
and the Cauchy distribution. This is not unexpected because the consistency of
the $\mathrm{dCov}$ based rely on certain moment assumptions which are not
satisfied by heavy-tailed distributions like the $t_{(2)}$ and the Cauchy
(recall that $t_{(2)}$ does not have finite variance and $\mathcal{C}(0,1)$
does not have finite mean). On the other hand, the $\mathrm{RJdCov}$ based
tests easily detects both the higher-order and joint dependences in the all 4
cases, which once again highlights the robustness and broad application of the
proposed method.

Table 2.  Empirical power for detecting higher-order dependence.  Dependence | Pariwise | Higher-Order | Joint
---|---|---|---
Tests | JdCov | Proposed | JdCov | Proposed | JdCov | Proposed
Gaussian | 0.058 | 0.054 | 1 | 1 | 1.000 | 0.942
$t_{(3)}$ | 0.052 | 0.064 | 1 | 1 | 0.868 | 0.904
$t_{(2)}$ | 0.056 | 0.040 | 1 | 1 | 0.098 | 0.888
Cauchy | 0.064 | 0.044 | 0.04 | 1 | 0.050 | 0.788

###  6.4. Performance of the ICA Estimator

To evaluate the performance of the distribution-free ICA estimator proposed in
Section 5, we consider the following model:

$\displaystyle X=\bm{M}Z,$

where $Z=(Z_{1},Z_{2},\ldots,Z_{d})$, $Z_{1},Z_{2},\ldots,Z_{d}$ i.i.d. for
some distribution $F$, and $M$ is a $d\times d$ random mixing matrix with
condition number between 1 and 2, using the R package
ProDenICA.555https://cran.r-project.org/web/packages/ProDenICA/index.html We
consider 12 candidate distributions for $F$ as described Table 3. For each
candidate distribution we apply the gradient descent algorithm to solve the
optimization program (5.6) (where the gradient is derived in Section 5.2). We
compare our results with the method proposed in Matteson and Tsay [52]. As in
[52] the estimation error is measured in the following metric, which is
designed to resolve the non-identifiability issues:

$\displaystyle
D(\hat{M},M)=\frac{1}{\sqrt{d-1}}\inf_{C\in\mathcal{C}}\|C\hat{M}^{-1}M-\bm{I}_{d}\|_{F}$

where

  * •

$\hat{M}$ is the estimated mixing matrix,

  * •

$\|\cdot\|_{F}$ denotes the Frobenius norm,

  * •

$\mathcal{C}=\\{C\in\mathcal{M}:C=P_{\pm}B\text{ for some }P_{\pm}\text{ and
}B\\}$, where $\mathcal{M}$ is the the set of $d\times d$ nonsingular
matrices, $P_{\pm}$ is a $d\times d$ signed permutation matrix, and $B$ is a
$d\times d$ diagonal matrix with positive diagonal elements.

The computation for metric $D$ is available in the R package JADE [57]. We set
dimension $d=3$ and the sample size $n=500$ and present the results in 4,
where the indices in $x$-axis denote different distributions as described in
Table 3. The results in Table 3 show that our estimator have smaller
estimation error and variance except for distribution (f).

Table 3. Description of the distributions for the ICA simulations and the means and standard deviations of the estimation errors.  |  | Estimation Error (Standard Deviation)
---|---|---
Index | Distribution | Proposed | MT
$\mathsf{(a)}$ | $\mathsf{N(0,1)^{3}}$ | 0.129 (0.077) | 0.192 (0.094)
$\mathsf{(b)}$ | $\mathsf{N(0,1)^{5}}$ | 0.170 (0.109) | 0.209 (0.120)
$\mathsf{(c)}$ | $\mathsf{Gamma(5,1)}$ | 0.067 (0.028) | 0.073 (0.049)
$\mathsf{(d)}$ | $\mathsf{Gamma(10,1)}$ | 0.127 (0.069) | 0.148 (0.088)
$\mathsf{(e)}$ | $\mathsf{\frac{3}{10}Exp(1)+\frac{7}{19}Exp(5)}$ | 0.034 (0.046) | 0.038 (0.055)
$\mathsf{(f)}$ | $\mathsf{\frac{3}{10}N(-2,1)+\frac{7}{10}N(2,1)}$ | 0.157 (0.121) | 0.181 (0.120)
$\mathsf{(g)}$ | $\mathsf{Uniform(0,1)}$ | 0.189 (0.110) | 0.310 (0.132)
$\mathsf{(h)}$ | $\mathsf{\frac{7}{10}N(-2,3)+\frac{3}{10}N(2,1)}$ | 0.115 (0.117) | 0.174 (0.134)
$\mathsf{(i)}$ | $\mathsf{\frac{5}{10}N(-2,2)+\frac{5}{10}N(2,2)}$ | 0.390 (0.105) | 0.435 (0.104)
$\mathsf{(j)}$ | $\mathsf{(\frac{5}{10}N(-2,2)+\frac{5}{10}N(2,2))^{3}}$ | 0.121 (0.085) | 0.205 (0.093)
$\mathsf{(k)}$ | $\mathsf{(\frac{5}{10}N(-2,2)+\frac{5}{10}N(2,2))^{5}}$ | 0.191 (0.097) | 0.236 (0.107)
$\mathsf{(l)}$ | $\mathsf{(\frac{5}{10}N(-2,2)+\frac{5}{10}N(2,2))^{7}}$ | 0.131 (0.084) | 0.157 (0.101)

Figure 4.  ICA estimation error across 12 distributions: The blue boxplots
correspond to the Matteson-Tsay estimator and the green boxplots corresponds
to the proposed estimator based on $\mathrm{RJdCov}$.

##  7\. Real data application

In this section, we apply our method to the stock price data available in R
package gmm. The data can be accessed using the command data(finance). The
dataset has prices for 24 stocks corresponding to 20 companies in US starting
from January 5th, 1993 to January 30th, 2009. Since the dataset has several
missing entries before 2003, we use data from the 2003-2009 period for 9 of
these companies (described below) for our analysis. To account for the
temporal dependence, we average the stock prices over the week, after which we
obtain $n=300$ data points. We categorize the different companies into 6
different industries according to North American Industry Classification
System (NAICS) as follows: (1) Real Estate (RE) category which contains the
company NHP, (2) the Food and Retail (FR) category which includes the
companies WMK and GCO; (3) Manufacture (M) category including ROG, MAT, and
VOXX; (4) the Finance (F) category with the company FNM; (5) Communications
Systems (JCS); and (6) Zoom Technologies (ZOOM). We created 2 separate
categories for the companies JCS and ZOOM because their stock prices seem to
vary independently of the other companies, possibly due to their relatively
small scales, even though they both belong to the information industry.

###  7.1. Results and Analysis

To understand the dependency structure between the 6 categories (industries)
described above, we first test the pairwise dependence of the stock prices
between different industries. Then we test for higher-order dependencies
between triples of categories for which the null hypothesis of pairwise
independence is accepted. We report the results for the tests based on
$\mathrm{RdCov}/\mathrm{RJdCov}$ and $\mathrm{dCov}/\mathrm{JdCov}$ for
comparison. Throughout we set the significance level to $0.05$.

The $p$-values for the pairwise independence tests are summarized in Table 4.
Note that this entails performing ${6\choose 2}=15$, so the $p$-values are
adjusted using the BH-procedure. To better visualize the results, we plot the
dependency graph between 6 categories using function dependence.structure in R
package multivariance [8] and the results are shown in Figure 5. From these we
can see that hypothesis of pairwise independence is accepted by
$\mathrm{RdCov}$ for the following 3 triples: (I) Real Estate-JCS-ZOOM, (II)
Finance-JCS-ZOOM, and (III) Manufacture-JCS-ZOOM. On the other hand,
$\mathrm{dCov}$ only accepts the pairwise independence hypothesis for the
triple: Real Estate-JCS-ZOOM when using $\mathrm{dCov}$. This might be due to
the small scale of companies JCS and ZOOM, which tends to make them
independent of each other and among others. Moreover, observe that the
Manufacture and Food and Retail categories are prone to be dependent with all
the others, which is also reasonable because their services are essential to
the other industries.

Table 4. $p$-values for pairwise independence testing among the 6 categories.  |  $p$-values
---|---
Industry-Industry | JdCov | Proposed
RE-JCS | 0.096 | 0.136
RE-ZOOM | 0.188 | 0.107
RE-FR | 0.002 | 0.000
RE-M | 0.002 | 0.000
RE-F | 0.002 | 0.000
JCS-ZOOM | 0.152 | 0.136
JCS-FR | 0.002 | 0.004
JCS-M | 0.002 | 0.258
JCS-F | 0.102 | 0.408
ZOOM-FR | 0.008 | 0.010
ZOOM-M | 0.012 | 0.058
ZOOM-F | 0.032 | 0.968
FR-M | 0.002 | 0.000
FR-F | 0.002 | 0.000
M-F | 0.002 | 0.000

The $p$-values for testing third-order independence for the 3 triples for
which the hypothesis of pairwise independence are accepted are given in Table
5. We observe from Table 5 that the third-order $\mathrm{dCov}$ based
$p$-value, which tests for the null hypothesis of joint independence of RE,
JCS, and ZOOM given their pairwise independence, accepts the null hypothesis
at level $0.05$. However, the $\mathrm{JdCov}$ based $p$-value, which tests
for the null hypothesis of joint independence of RE, JCS, and ZOOM (without
any assumptions on pairwise independence), rejects the null hypothesis of
joint independence at $0.05$. Given that the pairwise independence hypothesis
has been accepted, for consistent hierarchical interpretability one would have
ideally expected the third-order $\mathrm{dCov}$ and the $\mathrm{JdCov}$ to
have the same conclusions. On the other hand, the conclusions from the
$\mathrm{RdCov}$ and $\mathrm{RJdCov}$ $p$-values are hierarchically
consistent. Moreover, Figure 5(b) shows that there is higher-order dependence
between JCS, ZOOM, and Manufacture industries when using $\mathrm{RdCov}$ and
$\mathrm{RJdCov}$. This may be because JCS and ZOOM, and Manufacturing will
typically serve as downstream industries but direct links between them are not
obvious due to the scales of the companies.

(a)

(b)

Figure 5.  Stock prices data: (a) Dependency graph obtained using higher-order
$\mathrm{dCov}$, (b) Dependency graph obtained using higher-order
$\mathrm{RdCov}$.

.

Table 5. $p$-values for third order independence given pairwise independence.  | $\mathrm{dCov}$ | $\mathrm{JdCov}$ | $\mathrm{RdCov}$ | $\mathrm{RJdCov}$
---|---|---|---|---
Real estate-JCS- ZOOM | 0.101 | 0.020 | 0.187 | 0.063
Finance-JCS-ZOOM |  |  | 0.274 | 0.425
Manufacture-JCS-ZOOM |  |  | 0.003 | 0.006

### Acknowledgment

The authors are grateful to Björn Böttcher and David Matteson for sharing
their codes and datasets. BBB was partly supported by NSF CAREER Grant
DMS-2046393 and a Sloan research fellowship.

## References

  * Agarwal and Sharathkumar [2014] P. K. Agarwal and R. Sharathkumar.  Approximation algorithms for bipartite matching with metric and geometric costs.  In _STOC’14—Proceedings of the 2014 ACM Symposium on Theory of Computing_ , pages 555–564. ACM, New York, 2014.
  * Ai et al. [2022] C. Ai, L.-H. Sun, Z. Zhang, and L. Zhu.  Testing unconditional and conditional independence via mutual information.  _Journal of Econometrics_ , 2022.
  * Banerjee and Ghosh [2022] B. Banerjee and A. K. Ghosh.  Test of independence for hilbertian random variables.  _Stat_ , 11(1):e474, 2022.
  * Beran et al. [2007] R. Beran, M. Bilodeau, and P. L. de Micheaux.  Nonparametric tests of independence between random vectors.  _Journal of Multivariate Analysis_ , 98(9):1805–1824, 2007.
  * Bergsma and Dassios [2014] W. Bergsma and A. Dassios.  A consistent test of independence based on a sign covariance related to kendall’s tau.  _Bernoulli_ , 20(2):1006–1028, 2014.
  * Berrett and Samworth [2019] T. B. Berrett and R. J. Samworth.  Nonparametric independence testing via mutual information.  _Biometrika_ , 106(3):547–566, 2019.
  * Blum et al. [1961] J. R. Blum, J. Kiefer, and M. Rosenblatt.  Distribution free tests of independence based on the sample distribution function.  _Ann. Math. Statist._ , 32:485–498, 1961.
  * Böttcher [2020] B. Böttcher.  Dependence and dependence structures: estimation and visualization using the unifying concept of distance multivariance.  _Open Statistics_ , 1(1):1–48, 2020.
  * Böttcher et al. [2019] B. Böttcher, M. Keller-Ressel, and R. L. Schilling.  Distance multivariance: New dependence measures for random vectors.  _The Annals of Statistics_ , 47(5):2757–2789, 2019.
  * Brenier [1991] Y. Brenier.  Polar factorization and monotone rearrangement of vector-valued functions.  _Communications on pure and applied mathematics_ , 44(4):375–417, 1991.
  * Chacon and Rodriguez-Casal [2010] J. E. Chacon and A. Rodriguez-Casal.  A note on the universal consistency of the kernel distribution function estimator.  _Statistics & Probability Letters_, 80(17-18):1414–1419, 2010.
  * Chakraborty and Zhang [2019] S. Chakraborty and X. Zhang.  Distance metrics for measuring joint dependence with application to causal inference.  _Journal of the American Statistical Association_ , 2019.
  * Chatterjee [2007] S. Chatterjee.  Stein’s method, 2007.  URL https://souravchatterjee.su.domains/AllLectures.pdf.
  * Chatterjee [2021] S. Chatterjee.  A new coefficient of correlation.  _Journal of the American Statistical Association_ , 116(536):2009–2022, 2021.
  * Chatterjee [2022] S. Chatterjee.  A survey of some recent developments in measures of association.  _arXiv preprint arXiv:2211.04702_ , 2022.
  * Chernozhukov et al. [2017] V. Chernozhukov, A. Galichon, M. Hallin, and M. Henry.  Monge–kantorovich depth, quantiles, ranks and signs.  _The Annals of Statistics_ , 45(1):223–256, 2017\.
  * Csörgő [1985] S. Csörgő.  Testing for independence by the empirical characteristic function.  _Journal of Multivariate Analysis_ , 16(3):290–299, 1985.
  * Deb and Sen [2021] N. Deb and B. Sen.  Multivariate rank-based distribution-free nonparametric testing using measure transportation.  _Journal of the American Statistical Association_ , pages 1–16, 2021\.
  * Deb et al. [2020] N. Deb, P. Ghosal, and B. Sen.  Measuring association on topological spaces using kernels and geometric graphs.  _arXiv preprint arXiv:2010.01768_ , 2020.
  * Deb et al. [2021] N. Deb, B. B. Bhattacharya, and B. Sen.  Efficiency lower bounds for distribution-free hotelling-type two-sample tests based on optimal transport.  _arXiv preprint arXiv:2104.01986_ , 2021.
  * Feuerverger [1993] A. Feuerverger.  A consistent test for bivariate dependence.  _International Statistical Review/Revue Internationale de Statistique_ , pages 419–433, 1993.
  * Friedman and Rafsky [1983] J. H. Friedman and L. C. Rafsky.  Graph-theoretic measures of multivariate association and prediction.  _The Annals of Statistics_ , pages 377–391, 1983.
  * Gaißer et al. [2010] S. Gaißer, M. Ruppert, and F. Schmid.  A multivariate version of Hoeffding’s phi-square.  _Journal of Multivariate Analysis_ , 101(10):2571–2586, 2010.
  * Ghosal and Sen [2022] P. Ghosal and B. Sen.  Multivariate ranks and quantiles using optimal transport: Consistency, rates and nonparametric testing.  _The Annals of Statistics_ , 50(2):1012–1037, 2022.
  * Gieser [1993] P. W. Gieser.  _A new nonparametric test for independence between two sets of variates_.  PhD thesis, University of Florida, 1993.
  * Gieser and Randles [1997] P. W. Gieser and R. H. Randles.  A nonparametric test of independence between two vectors.  _Journal of the American Statistical Association_ , 92(438):561–567, 1997.
  * Gretton et al. [2005a] A. Gretton, O. Bousquet, A. Smola, and B. Schölkopf.  Measuring statistical dependence with hilbert-schmidt norms.  In _International conference on algorithmic learning theory_ , pages 63–77. Springer, 2005a.
  * Gretton et al. [2005b] A. Gretton, R. Herbrich, A. Smola, O. Bousquet, B. Schölkopf, et al.  Kernel methods for measuring independence.  2005b.
  * Gretton et al. [2007] A. Gretton, K. Fukumizu, C. Teo, L. Song, B. Schölkopf, and A. Smola.  A kernel statistical test of independence.  _Advances in neural information processing systems_ , 20, 2007.
  * Hallin [2017] M. Hallin.  On distribution and quantile functions, ranks and signs in $\mathbb{R}^{d}$.  _ECARES Working Papers_ , 2017.
  * Hallin et al. [2020a] M. Hallin, E. del Barrio, J. A. Cuesta-Albertos, and C. Matrán.  On distribution and quantile functions, ranks, and signs in $\mathbb{R}^{d}$: a measure transportation approach.  _Ann.Stat. (to appear)_ , 2020a.
  * Hallin et al. [2020b] M. Hallin, D. La Vecchia, and H. Liu.  Center-outward R-estimation for semiparametric VARMA models.  _Journal of the American Statistical Association_ , pages 1–14, 2020b.
  * Hallin et al. [2022] M. Hallin, D. Hlubinka, and Š. Hudecová.  Efficient fully distribution-free center-outward rank tests for multiple-output regression and manova.  _Journal of the American Statistical Association_ , pages 1–17, 2022\.
  * Heller et al. [2012] R. Heller, M. Gorfine, and Y. Heller.  A class of multivariate distribution-free tests of independence based on graphs.  _Journal of Statistical Planning and Inference_ , 142(12):3097–3106, 2012.
  * Heller et al. [2013] R. Heller, Y. Heller, and M. Gorfine.  A consistent multivariate test of association based on ranks of distances.  _Biometrika_ , 100(2):503–510, 2013.
  * Heller et al. [2016] R. Heller, Y. Heller, S. Kaufman, B. Brill, and M. Gorfine.  Consistent distribution-free $k$-sample and independence tests for univariate random variables.  _The Journal of Machine Learning Research_ , 17(1):978–1031, 2016.
  * Hoeffding [1948] W. Hoeffding.  A Non-Parametric Test of Independence.  _The Annals of Mathematical Statistics_ , 19(4):546 – 557, 1948.
  * Hofer [2009] R. Hofer.  On the distribution properties of Niederreiter-Halton sequences.  _J. Number Theory_ , 129(2):451–463, 2009.  ISSN 0022-314X.
  * Hyvärinen and Oja [2000] A. Hyvärinen and E. Oja.  Independent component analysis: algorithms and applications.  _Neural networks_ , 13(4-5):411–430, 2000.
  * Jonker and Volgenant [1987] R. Jonker and A. Volgenant.  A shortest augmenting path algorithm for dense and sparse linear assignment problems.  _Computing_ , 38(4):325–340, 1987.
  * Josse and Holmes [2016] J. Josse and S. Holmes.  Measuring multivariate association and beyond.  _Stat. Surv._ , 10:132–167, 2016.
  * Kankainen [1995] A. Kankainen.  _Consistent testing of total independence based on the empirical characteristic function_.  University of Jyväskylä, 1995.
  * Kinney and Atwal [2014] J. B. Kinney and G. S. Atwal.  Equitability, mutual information, and the maximal information coefficient.  _Proceedings of the National Academy of Sciences_ , 111(9):3354–3359, 2014.
  * Konijn [1956] H. Konijn.  On the power of certain tests for independence in bivariate populations.  _The Annals of Mathematical Statistics_ , 27(2):300–323, 1956.
  * Kosorok [2008] M. R. Kosorok.  _Introduction to empirical processes and semiparametric inference._ Springer, 2008.
  * Kuo [1975] H.-H. Kuo.  Gaussian measures in banach spaces.  _Gaussian Measures in Banach Spaces_ , pages 1–109, 1975.
  * Lehmann and Romano [2005] E. L. Lehmann and J. P. Romano.  _Testing Statistical Hypotheses_.  Springer Texts in Statistics. Springer, New York, third edition, 2005\.  ISBN 0-387-98864-5.
  * Li et al. [2012] R. Li, W. Zhong, and L. Zhu.  Feature screening via distance correlation learning.  _Journal of the American Statistical Association_ , 107(499):1129–1139, 2012.
  * Ma [2022] J. Ma.  Evaluating independence and conditional independence measures.  _arXiv preprint arXiv:2205.07253_ , 2022.
  * Ma and Mao [2019] L. Ma and J. Mao.  Fisher exact scanning for dependency.  _Journal of the American Statistical Association_ , 114(525):245–258, 2019.
  * Matteson [2008] D. S. Matteson.  _Statistical inference for multivariate nonlinear time series_.  The University of Chicago, 2008.
  * Matteson and Tsay [2017] D. S. Matteson and R. S. Tsay.  Independent component analysis via distance covariance.  _Journal of the American Statistical Association_ , 112(518):623–637, 2017.
  * McCann [1995] R. J. McCann.  Existence and uniqueness of monotone measure-preserving maps.  _Duke Mathematical Journal_ , 80(2):309–323, 1995\.
  * Monge [1781] G. Monge.  Mémoire sur la théorie des déblais et des remblais.  _Histoire de l’Académie Royale des Sciences de Paris_ , 1781.
  * Moon and Chen [2022] H. Moon and K. Chen.  Interpoint-ranking sign covariance for the test of independence.  _Biometrika_ , 109(1):165–179, 2022.
  * Nandy et al. [2016] P. Nandy, L. Weihs, and M. Drton.  Large-sample theory for the bergsma-dassios sign covariance.  _Electronic Journal of Statistics_ , 10(2):2287–2311, 2016.
  * Nordhausen et al. [2014] K. Nordhausen, J. Cardoso, J. Miettinen, H. Oja, E. Ollila, and S. Taskinen.  Jade: Jade and other bss methods as well as some bss performance criteria.  _R package version_ , pages 1–9, 2014.
  * Pan et al. [2019] W. Pan, X. Wang, H. Zhang, H. Zhu, and J. Zhu.  Ball covariance: A generic measure of dependence in banach space.  _Journal of the American Statistical Association_ , 2019.
  * Pfister et al. [2018] N. Pfister, P. Bühlmann, B. Schölkopf, and J. Peters.  Kernel-based tests for joint independence.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(1):5–31, 2018.
  * Póczos et al. [2012] B. Póczos, Z. Ghahramani, and J. Schneider.  Copula-based kernel dependency measures.  _arXiv preprint arXiv:1206.4682_ , 2012.
  * Reshef et al. [2011] D. N. Reshef, Y. A. Reshef, H. K. Finucane, S. R. Grossman, G. McVean, P. J. Turnbaugh, E. S. Lander, M. Mitzenmacher, and P. C. Sabeti.  Detecting novel associations in large data sets.  _science_ , 334(6062):1518–1524, 2011.
  * Resnick [2013] S. I. Resnick.  _A Probability Path_.  Springer Science & Business Media, 2013.
  * Roy et al. [2020] A. Roy, S. Sarkar, A. K. Ghosh, and A. Goswami.  On some consistent tests of mutual independence among several random vectors of arbitrary dimensions.  _Statistics and Computing_ , 30(6):1707–1723, 2020.
  * Sejdinovic et al. [2013a] D. Sejdinovic, A. Gretton, and W. Bergsma.  A kernel test for three-variable interactions.  _Advances in Neural Information Processing Systems_ , 26, 2013a.
  * Sejdinovic et al. [2013b] D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu.  Equivalence of distance-based and RKHS-based statistics in hypothesis testing.  _Ann. Statist._ , 41(5):2263–2291, 2013b.
  * Sejdinovic et al. [2013c] D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu.  Equivalence of distance-based and rkhs-based statistics in hypothesis testing.  _The annals of statistics_ , pages 2263–2291, 2013c.
  * Shi et al. [2020a] H. Shi, M. Drton, and F. Han.  Distribution-free consistent independence tests via center-outward ranks and signs.  _Journal of the American Statistical Association_ , pages 1–16, 2020a.
  * Shi et al. [2020b] H. Shi, M. Hallin, M. Drton, and F. Han.  On universally consistent and fully distribution-free rank tests of vector independence.  _arXiv preprint arXiv:2007.02186_ , 2020b.
  * Székely and Rizzo [2009] G. J. Székely and M. L. Rizzo.  Brownian distance covariance.  _The annals of applied statistics_ , 3(4):1236–1265, 2009.
  * Székely et al. [2007] G. J. Székely, M. L. Rizzo, and N. K. Bakirov.  Measuring and testing dependence by correlation of distances.  _The annals of statistics_ , 35(6):2769–2794, 2007.
  * Van der Vaart [2000] A. W. Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Van Der Vaart and Wellner [1996] A. W. Van Der Vaart and J. Wellner.  _Weak convergence and empirical processes: with applications to statistics_.  Springer Science & Business Media, 1996.
  * Villani [2003] C. Villani.  _Topics in optimal transportation_ , volume 58 of _Graduate Studies in Mathematics_.  American Mathematical Society, Providence, RI, 2003.  ISBN 0-8218-3312-X.
  * Wu et al. [2009] E. H. Wu, L. Philip, and W. K. Li.  A smoothed bootstrap test for independence based on mutual information.  _Computational statistics & data analysis_, 53(7):2524–2536, 2009.
  * Yanagimoto [1970] T. Yanagimoto.  On measures of association and a related problem.  _Annals of the Institute of Statistical Mathematics_ , 22(1):57–63, 1970.
  * Yao et al. [2018] S. Yao, X. Zhang, and X. Shao.  Testing mutual independence in high dimension via distance covariance.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(3):455–480, 2018.
  * Zhang [2019] K. Zhang.  Bet on independence.  _Journal of the American Statistical Association_ , 114(528):1620–1637, 2019.
  * Zhang et al. [2018] Q. Zhang, S. Filippi, A. Gretton, and D. Sejdinovic.  Large-scale kernel methods for independence testing.  _Statistics and Computing_ , 28(1):113–130, 2018\.
  * Zhu et al. [2020] H. Zhu, R. Li, R. Zhang, and H. Lian.  Nonlinear functional canonical correlation analysis via distance covariance.  _Journal of Multivariate Analysis_ , 180:104662, 2020.

##  Appendix A Combinatorial CLT with multiple permutations

In this section we derive an analogue of Hoeffding’s combinatorial CLT [37]
for tensors with multiple random permutations. We begin with the following
assumption:

######  Assumption 5.

Fix $r\geq 2$ and suppose $\bm{A}_{n}=((a_{i_{1},i_{2},\ldots,i_{r}}))_{1\leq
i_{1},i_{2},\ldots,i_{r}\leq n}$ is a sequence of $r$-tensors satisfying the
following conditions:

  * $(1)$

For all $1\leq s\leq r$ and $1\leq
i_{1},i_{2},\ldots,i_{s-1},i_{s+1},\ldots,i_{r}\leq n$,

$\displaystyle\sum_{i_{s}=1}^{n}a_{i_{1},i_{2},\cdots,i_{r}}=0.$ (A.1)

  * $(2)$

There is a universal constants $K_{1},K_{2}>0$ such that
$|a_{i_{1},i_{2},\ldots,i_{r}}|\leq K_{1}/\sqrt{n}$, for all $1\leq
i_{1},i_{2},\ldots,i_{r}\leq n$, and $\sum_{1\leq i_{1},i_{2},\ldots,i_{r}\leq
n}a_{i_{1},i_{2},\ldots,i_{r}}^{2}\geq K_{2}n^{r-1}$.

Under this assumption we have the following theorem:

######  Theorem A.1.

Fix $r\geq 2$ and suppose $\bm{A}_{n}=((a_{i_{1},i_{2},\ldots,i_{r}}))_{1\leq
i_{1},i_{2},\cdots,i_{r}\leq n}$ is a sequence of $r$-tensors satisfying
Assumption 5. Consider $r-1$ independent random permutations
$\pi_{1},\pi_{2},\ldots,\pi_{r-1}$ of the set $\\{1,\ldots,n\\}$ and define

$C_{n}:=\sum_{i=1}^{n}a_{i,\pi_{1}(i),\ldots,\pi_{r-1}(i)}.$

Then, as $n\rightarrow\infty$,

$\displaystyle\frac{C_{n}}{\sqrt{\operatorname{Var}[C_{n}]}}\stackrel{{\scriptstyle
D}}{{\to}}N(0,1).$

The proof of Theorem A.1 is given in Section A.1. We first show that by
centering $\bm{A}_{n}$, condition (1) in Assumption 5 holds without loss of
generality.

######  Remark A.2.

Suppose $\bm{A}_{n}$ is an $r$-tensor not satisfying (A.1), then we can
construct a tensor
$\tilde{\bm{A}_{n}}=((\tilde{a}_{i_{1},i_{2},\ldots,i_{r}}))_{1\leq
i_{1},i_{2},\ldots,i_{r}\leq n}$ which satisfies Assumption 5. We illustrate
this for $r=3$. For $1\leq i_{1},i_{2},i_{3}\leq r$, define

$\displaystyle\tilde{a}_{i_{1},i_{2},i_{3}}=a_{i_{1},i_{2},i_{3}}-a_{\bullet,i_{2},i_{3}}-a_{i_{1},\bullet,i_{3}}-a_{i_{1},i_{2},\bullet}+a_{i_{1},\bullet,\bullet}+a_{i_{2},\bullet,\bullet}+a_{i_{3},\bullet,\bullet}-a_{\bullet,\bullet,\bullet}$

where

  * •

$a_{\bullet,i_{2},i_{3}}=\frac{1}{n}\sum_{i_{1}=1}^{n}a_{i_{1},i_{2},i_{3}}$,
$a_{i_{1},\bullet,i_{3}}=\frac{1}{n}\sum_{i_{2}=1}^{n}a_{i_{1},i_{2},i_{3}}$,
$a_{i_{1},i_{2},\bullet}=\frac{1}{n}\sum_{i_{3}=1}^{n}a_{i_{1},i_{2},i_{3}}$;

  * •

$a_{i_{1},\bullet,\bullet}=\frac{1}{n^{2}}\sum_{1\leq i_{2},i_{3}\leq
n}a_{i_{1},i_{2},i_{3}}$,
$a_{\bullet,i_{2},\bullet}=\frac{1}{n^{2}}\sum_{1\leq i_{1},i_{3}\leq
n}a_{i_{1},i_{2},i_{3}}$,
$a_{\bullet,\bullet,i_{3}}=\frac{1}{n^{2}}\sum_{1\leq i_{1},i_{2}\leq
n}a_{i_{1},i_{2},i_{3}}$;

  * •

$a_{\bullet,\bullet,\bullet}=\frac{1}{n^{3}}\sum_{1\leq i_{1},i_{2},i_{3}\leq
n}a_{i_{1},i_{2},i_{3}}$.

Clearly, $\tilde{\bm{A}_{n}}=((\tilde{a}_{i_{1},i_{2},i_{3}}))_{1\leq
i_{1},i_{2},i_{3}\leq n}$ satisfies (A.1).

###  A.1. Proof of Theorem A.1

For notational convenience we will present the proof of Theorem A.1 for $r=3$.
The proof can be easily extended to general $r$ by following the arguments.
Hereafter, we will assume $\bm{A}_{n}=((a_{i_{1},i_{2},i_{3}}))_{1\leq
i_{1},i_{2},\cdots,i_{r}\leq n}$ is a sequence of $r$-tensors satisfying
Assumption 5. We begin by computing $\operatorname{Var}[C_{n}]$.

######  Lemma A.3.

Suppose Assumption 5 holds. Then

$\displaystyle\operatorname{Var}[C_{n}]=\frac{n-2}{n(n-1)^{2}}\sum_{1\leq
i_{1},i_{2},i_{3}\leq n}a_{i_{1},i_{2},i_{3}}^{2}=\Theta(1).$

###### Proof.

Note that

$\displaystyle\operatorname{Var}[C_{n}]=\sum_{i=1}^{n}\operatorname{Var}[a_{i,\pi_{1}(i),\pi_{2}(i)}]+\sum_{1\leq
i_{1}\neq i_{1}^{\prime}\leq
n}\operatorname{Cov}[a_{i_{1},\pi_{1}(i_{1}),\pi_{2}(i_{1})},a_{i_{1}^{\prime},\pi_{1}(i_{1}^{\prime}),\pi_{2}(i_{1}^{\prime})}].$
(A.2)

Assumption 5 implies that $\mathbb{E}[a_{i,\pi_{1}(i),\pi_{2}(i)}]=0$, hence,

$\displaystyle\sum_{i=1}^{n}\operatorname{Var}[a_{i,\pi_{1}(i),\pi_{2}(i)}]=\sum_{i=1}^{n}\mathbb{E}[a_{i,\pi_{1}(i),\pi_{2}(i)}^{2}]=\frac{1}{n^{2}}\sum_{1\leq
i_{1},i_{2},i_{3}\leq n}a_{i_{1},i_{2},i_{3}}.$ (A.3)

Moreover,

$\displaystyle\operatorname{Cov}[a_{i_{1},\pi_{1}(i_{1}),\pi_{2}(i_{1})},a_{i_{1}^{\prime},\pi_{1}(i_{1}^{\prime}),\pi_{2}(i_{1}^{\prime})}]$
$\displaystyle=\mathbb{E}[a_{i_{1},\pi_{1}(i_{1}),\pi_{2}(i_{1})}a_{i_{1}^{\prime}\pi_{1}(i_{1}^{\prime})\pi_{2}(i^{\prime}_{1})}]$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{2}\neq
i_{2}^{\prime}\leq n}\sum_{1\leq i_{3}\neq i_{3}^{\prime}\leq
n}a_{i_{1},i_{2},i_{3}}a_{i_{1}^{\prime},i_{2}^{\prime},i_{3}^{\prime}}$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{2},i_{3}\leq
n}\left\\{a_{i_{1},i_{2},i_{3}}\sum_{\begin{subarray}{c}1\leq
i_{2}^{\prime}\leq n\\\ i_{2}^{\prime}\neq
i_{2}\end{subarray}}\sum_{\begin{subarray}{c}1\leq i_{3}^{\prime}\leq n\\\
i_{3}^{\prime}\neq
i_{3}\end{subarray}}a_{i_{1}^{\prime},i_{2}^{\prime},i_{3}^{\prime}}\right\\}$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{2},i_{3}\leq
n}\left\\{a_{i_{1},i_{2},i_{3}}\sum_{\begin{subarray}{c}1\leq
i_{2}^{\prime}\leq n\\\ i_{2}^{\prime}\neq
i_{2}\end{subarray}}(-a_{i_{1}^{\prime},i_{2}^{\prime},i_{3}})\right\\}$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{2},i_{3}\leq
n}\left\\{a_{i_{1},i_{2},i_{3}}a_{i_{1}^{\prime},i_{2},i_{3}}\right\\}.$

Hence,

$\displaystyle\sum_{1\leq i_{1}\neq i_{1}^{\prime}\leq
n}\operatorname{Cov}[a_{i_{1},\pi_{1}(i_{1}),\pi_{2}(i_{1})},a_{i_{1}^{\prime},\pi_{1}(i_{1}^{\prime}),\pi_{2}(i_{1}^{\prime})}]$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{1}\neq
i_{1}^{\prime}\leq n}\sum_{1\leq i_{2},i_{3}\leq
n}\left\\{a_{i_{1},i_{2},i_{3}}a_{i_{1}^{\prime},i_{2},i_{3}}\right\\}$
$\displaystyle=\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{1},i_{2},i_{3}\leq
n}a_{i_{1},i_{2},i_{3}}\sum_{\begin{subarray}{c}1\leq i_{1}^{\prime}\leq n\\\
i_{1}^{\prime}\neq i_{1}\end{subarray}}a_{i_{1}^{\prime},i_{2},i_{3}}$
$\displaystyle=-\frac{1}{n^{2}(n-1)^{2}}\sum_{1\leq i_{1},i_{2},i_{3}\leq
n}a_{i_{1},i_{2},i_{3}}^{2}.$ (A.4)

Therefore, combining (A.3) and (A.1) with (A.2) the lemma follows. ∎

Since $\operatorname{Var}[C_{n}]=\Theta(1)$ by Assumption 5, we can always
normalize the elements of $\bm{A}_{n}$ such that
$\operatorname{Var}[C_{n}]=1$. Hereafter, we will assume $\bm{A}_{n}$ is
normalized such that $\operatorname{Var}[C_{n}]=1$. To prove the CLT in
Theorem A.1 will use Stein’s method based on exchangeable pairs [13]. The
first step towards this is to construct an exchangeable pair for the random
permutations $(\pi_{1},\pi_{2})$. This is done as follows:

  * •

We define $\pi_{1}^{\prime}=\pi_{1}\circ(I,J)$. Choose a transposition $(I,J)$
uniformly at random from the set of transpositions on $\\{1,\ldots,n\\}$ and
define $\pi_{1}^{\prime}:=\pi_{1}\circ(I,J)$, that is,
$\pi_{1}^{\prime}(I)=\pi_{1}(J)$, $\pi_{1}^{\prime}(J)=\pi_{1}(I)$, and
$\pi_{1}^{\prime}(\ell)=\pi_{1}(\ell)$ for $\ell\neq\\{I,J\\}$.

  * •

Similarly, choose an independent transposition $(K,L)$ uniformly at random
from the set of transpositions on $\\{1,\ldots,n\\}$ and define
$\pi_{2}^{\prime}:=\pi_{2}\circ(K,L)$.

######  Lemma A.4.

$\\{(\pi_{1},\pi_{2}),(\pi_{1}^{\prime},\pi_{2}^{\prime})\\}$ is an
exchangeable pair.

###### Proof.

Denote by $S_{n}$ the collection of all permutations of $\\{1,2,\ldots,n\\}$
and suppose $A_{1},A_{2},A_{1}^{\prime}A_{2}^{\prime}\subseteq S_{n}$. Then by
the independence between $(\pi_{1},\pi_{1}^{\prime})$ and
$(\pi_{2},\pi_{2}^{\prime})$,

$\displaystyle\mathbb{P}((\pi_{1},\pi_{2})\in A_{1}\times
A_{2},(\pi_{1}^{\prime},\pi_{2}^{\prime})\in A_{1}^{\prime}\times
A_{2}^{\prime})$ $\displaystyle=\mathbb{P}(\pi_{1}\in
A_{1},\pi_{1}^{\prime}\in A_{1}^{\prime})\mathbb{P}(\pi_{2}\in
A_{2},\pi_{2}^{\prime}\in A_{2}^{\prime})$
$\displaystyle=\mathbb{P}(\pi_{1}^{\prime}\in A_{1},\pi_{1}\in
A_{1}^{\prime})\mathbb{P}(\pi_{2}^{\prime}\in A_{2},\pi_{2}\in
A_{2}^{\prime})$
$\displaystyle=\mathbb{P}((\pi_{1}^{\prime},\pi_{2}^{\prime})\in A_{1}\times
A_{2},(\pi_{1},\pi_{2})\in A_{1}^{\prime}\times A_{2}^{\prime}),$

where the second inequality uses the fact that $(\pi_{1},\pi_{1}^{\prime})$
and $(\pi_{2},\pi_{2}^{\prime})$ are marginally two exchangeable pairs. ∎

Recall that $C_{n}=\sum_{i=1}^{n}a_{i\pi_{1}(i)\pi_{2}(i)}$ and define

$\displaystyle
C_{n}^{\prime}=\sum_{i=1}^{n}a_{i\pi_{1}^{\prime}(i)\pi_{2}^{\prime}(i)}.$
(A.5)

By Lemma A.4, $(C_{n},C_{n}^{\prime})$ is an exchangeable pair. Hence, by
Stein’s method based on exchangeable pairs we can obtain bounds on the
Wasserstein distance between $C_{n}$ and $N(0,1)$ in terms of the moments of
the difference $C_{n}^{\prime}-C_{n}$. To this end, recall that the
Wasserstein distance between random variables $X\sim\mu$ and $Y\sim\nu$ on
$\mathbb{R}$ is defined as

$\mathrm{Wass}(X,Y)=\sup\left\\{\left|\int f\mathrm{d}\mu-\int
f\mathrm{d}\nu\right|:f\text{ is }1-\text{Lipschitz}\right\\},$

where a function $f:\mathbb{R}\rightarrow\mathbb{R}$ is 1-Lipschitz if
$|f(x)-f(y)|\leq|x-y|$, for all $x,y\in\mathbb{R}$. We now invoke the
following theorem:

######  Theorem A.5.

[13] Let $(C_{n},C_{n}^{\prime})$ be as defined above. If
$\mathbb{E}[C_{n}^{\prime}-C_{n}|C_{n}]=-\lambda C_{n}$, for some
$0<\lambda<1$, and $\mathbb{E}[C_{n}^{2}]=1$, then

$\displaystyle\mathrm{Wass}(C_{n},N(0,1))\leq\sqrt{\frac{2}{\pi}\operatorname{Var}\left[\mathbb{E}\left[\frac{(C_{n}^{\prime}-C_{n})^{2}}{2\lambda}\Big{|}C_{n}\right]\right]}+\frac{1}{3\lambda}\mathbb{E}\left[|C_{n}^{\prime}-C_{n}|^{3}\right].$

The rest of the proof is organized as follows:

  * (1)

In Lemma A.6 we evaluate $\mathbb{E}[C_{n}^{\prime}-C_{n}|C_{n}]$ and show
that $\mathbb{E}[C_{n}^{\prime}-C_{n}|C_{n}]=-\lambda C_{n}$ holds with
$\lambda=\Theta(1/n)$.

  * (2)

Next, in Lemma A.7 we show
$\frac{1}{\lambda}\mathbb{E}[|C_{n}^{\prime}-C_{n}|^{3}]=O(1/\sqrt{n})$.

  * (3)

Finally, in Lemma A.8 we show that
$\operatorname{Var}[\mathbb{E}[\frac{(C_{n}^{\prime}-C_{n})^{2}}{2\lambda}|C_{n}]]=O(1/n)$.

These three steps combined with Theorem A.5 implies,
$\mathrm{Wass}(C_{n},N(0,1))=O(1/\sqrt{n})$, thus, completing the proof of
Theorem A.1.

######  Lemma A.6.

$\mathbb{E}[C_{n}^{\prime}-C_{n}|C_{n}]=-\lambda C_{n}$, where

$\displaystyle\lambda=\frac{4(n^{2}-7n+11)}{n^{2}(n-1)}=\Theta\left(\frac{1}{n}\right).$

###### Proof.

Note that

$\displaystyle
C_{n}^{\prime}-C_{n}=\begin{cases}\Delta_{1}(I,J)-\Delta_{0}(I,J)&\text{if
}K=I,L=J\text{ or }K=J,L=I,\\\ \Delta_{2}(I,J,K)-\Delta_{0}(I,J,K)&\text{if
}K\neq J,L=I,\\\ \Delta_{2}(I,J,L)-\Delta_{0}(I,J,L)&\text{if }K=I,L\neq J,\\\
\Delta_{3}(I,J,K)-\Delta_{0}(I,J,K)&\text{if }K\neq I,L=J,\\\
\Delta_{3}(I,J,L)-\Delta_{0}(I,J,L)&\text{if }K=J,L\neq I,\\\
\Delta_{4}(I,J,K,L)-\Delta_{0}(I,J,K,L)&\text{otherwise;}\end{cases}$ (A.6)

where

$\displaystyle\Delta_{0}(I,J)$
$\displaystyle:=a_{I,\pi_{1}(I),\pi_{2}(I)}+a_{J,\pi_{1}(J),\pi_{2}(J)},$
$\displaystyle\Delta_{0}(I,J,K)$
$\displaystyle:=a_{I,\pi_{1}(I),\pi_{2}(I)}+a_{J,\pi_{1}(J),\pi_{2}(J)}+a_{K,\pi_{1}(K),\pi_{2}(K)},$
$\displaystyle\Delta_{0}(I,J,K,L)$
$\displaystyle:=a_{I,\pi_{1}(I),\pi_{2}(I)}+a_{J,\pi_{1}(J),\pi_{2}(J)}+a_{K,\pi_{1}(K),\pi_{2}(K)}+a_{L,\pi_{1}(L),\pi_{2}(L)};$

and

$\displaystyle\Delta_{1}(I,J)$
$\displaystyle:=a_{I,\pi_{1}(J),\pi_{2}(J)}+a_{J,\pi_{1}(I),\pi_{2}(I)}$
$\displaystyle\Delta_{2}(I,J,\cdot)$
$\displaystyle:=a_{I,\pi_{1}(J),\pi_{2}(\cdot)}+a_{J,\pi_{1}(I),\pi_{2}(J)}+a_{\cdot,\pi_{1}(\cdot),\pi_{2}(I)}$
$\displaystyle\Delta_{3}(I,J,\cdot)$
$\displaystyle:=a_{I,\pi_{1}(J),\pi_{2}(I)}+a_{J,\pi_{1}(I),\pi_{2}(\cdot)}+a_{\cdot,\pi_{1}(\cdot),\pi_{2}(J)}$
$\displaystyle\Delta_{4}(I,J,K,L)$
$\displaystyle:=a_{I,\pi_{1}(J),\pi_{2}(I)}+a_{J,\pi_{1}(I),\pi_{2}(J)}+a_{K,\pi_{1}(K),\pi_{2}(L)}+a_{L,\pi_{1}(L),\pi_{2}(K)}.$

Denote $\eta=\frac{1}{n(n-1)}$. Then (A.6) implies,

$\displaystyle\mathbb{E}[C_{n}^{\prime}-C_{n}|\pi_{1},\pi_{2}]=\eta^{2}\left(S_{1}+S_{2}+S_{3}+S_{4}\right)$

where

$\displaystyle S_{1}$ $\displaystyle:=2\sum_{1\leq i\neq j\leq
n}\left\\{\Delta_{1}(i,j)-\Delta_{0}(i,j)\right\\}$ $\displaystyle S_{2}$
$\displaystyle:=\sum_{1\leq i\neq j\neq k\leq
n}\left\\{\Delta_{2}(i,j,k)+\Delta_{3}(i,j,k)-2\Delta_{0}(i,j,k)\right\\}$
$\displaystyle S_{3}$ $\displaystyle:=\sum_{1\leq i\neq j\neq\ell\leq
n}\left\\{\Delta_{2}(i,j,\ell)+\Delta_{3}(i,j,\ell)-2\Delta_{0}(i,j,\ell)\right\\}$
$\displaystyle S_{4}$ $\displaystyle:=\sum_{1\leq i\neq j\neq k\neq\ell\leq
n}\left\\{\Delta_{4}(i,j,k,\ell)-\Delta_{0}(i,j,k,\ell)\right\\}.$

Note that, by Assumption 5, $\sum_{1\leq i\neq j\leq
n}a_{i,\pi_{1}(j),\pi_{2}(j)}=-\sum_{j=1}^{n}a_{j,\pi_{1}(j),\pi_{2}(j)}=-C_{n}$.
Hence,

$\displaystyle S_{1}=-4C_{n}-4(n-1)C_{n}.$

Computing the other terms similarly and simplifying gives

$\displaystyle\mathbb{E}[C_{n}^{\prime}-C_{n}|C_{n}]=\mathbb{E}[C_{n}^{\prime}-C_{n}|\pi_{1},\pi_{2}]=-\lambda
C_{n}$

where

$\displaystyle\lambda=\frac{4(n-2)(n-3)}{n(n-1)^{2}}+\frac{4(n-2)(5-3n)}{n^{2}(n-1)^{2}}+\frac{4}{n^{2}(n-1)}=\frac{4(n^{2}-7n+11)}{n^{2}(n-1)}=\Theta\left(\frac{1}{n}\right).$

This completes the proof of Lemma A.6. ∎

######  Lemma A.7.

$\frac{1}{\lambda}\mathbb{E}[|C_{n}^{\prime}-C_{n}|^{3}]=O(1/\sqrt{n})$.

###### Proof.

Denote $\eta=\frac{1}{n(n-1)}$. Then (A.6) implies,

$\displaystyle\mathbb{E}[|C_{n}^{\prime}-C_{n}|^{3}|\pi_{1},\pi_{2}]\lesssim\eta^{2}(W_{1}+W_{2}+W_{3}),$

where

$\displaystyle W_{1}$ $\displaystyle:=\sum_{1\leq i\neq j\leq
n}\left(|\Delta_{1}(i,j)|^{3}+|\Delta_{0}(i,j)|^{3}\right),$ $\displaystyle
W_{2}$ $\displaystyle:=\sum_{1\leq i\neq j\neq k\leq
n}\left(|\Delta_{2}(i,j,k)|^{3}+|\Delta_{3}(i,j,k)|^{3}+|\Delta_{0}(i,j,k)|^{3}\right),$
$\displaystyle W_{3}$ $\displaystyle:=\sum_{1\leq i\neq j\neq k\neq\ell\leq