---
abstract: 'Understanding and developing a correlation measure that can detect general dependencies is not only imperative to statistics and machine learning, but also crucial to general scientific discovery in the big data age. In this paper, we establish a new framework that generalizes distance correlation — a correlation measure that was recently proposed and shown to be universally consistent for dependence testing against all joint distributions of finite moments — to the Multiscale Graph Correlation ([[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}). By utilizing the characteristic functions and incorporating the nearest neighbor machinery, we formalize the population version of local distance correlations, define the optimal scale in a given dependency, and name the optimal local correlation as [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}. The new theoretical framework motivates a theoretically sound Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and allows a number of desirable properties to be proved, including the universal consistency, convergence and almost unbiasedness of the sample version. The advantages of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} are illustrated via a comprehensive set of simulations with linear, nonlinear, univariate, multivariate, and noisy dependencies, where it loses almost no power in monotone dependencies while achieving better performance in general dependencies, compared to distance correlation and other popular methods.'
author:
- 'Cencheng Shen[^1]'
- 'Carey E. Priebe[^2]'
- 'Joshua T. Vogelstein[^3]'
bibliography:
- 'MGCbib.bib'
title: '**From Distance Correlation to Multiscale Graph Correlation**'
---

\#1

1

[1]{}

0

[1]{}

[**From Distance Correlation to Multiscale Graph Correlation**]{}

[*Keywords:*]{} testing independence, generalized distance correlation, nearest neighbor graph

Introduction
============

Given pairs of observations $(x_{i},y_{i}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{q}$ for $i=1,\ldots,n$, assume they are generated by independently identically distributed (*iid*) $F_{{\ensuremath{X}}{\ensuremath{Y}}}$. A fundamental statistical question prior to the pursuit of any meaningful joint inference is the independence testing problem: the two random variables are independent if and only if $F_{{\ensuremath{X}}{\ensuremath{Y}}} = F_{{\ensuremath{X}}} F_{{\ensuremath{Y}}}$, i.e., the joint distribution equals the product of the marginals. The statistical hypothesis is formulated as: $$\begin{aligned}
& H_{0}: F_{{\ensuremath{X}}{\ensuremath{Y}}}=F_{{\ensuremath{X}}}F_{{\ensuremath{Y}}},\\
& H_{A}: F_{{\ensuremath{X}}{\ensuremath{Y}}} \neq F_{{\ensuremath{X}}}F_{{\ensuremath{Y}}}.\end{aligned}$$ For any test statistic, the testing power at a given type $1$ error level equals the probability of correctly rejecting the null hypothesis when the random variables are dependent. A test is consistent if and only if the testing power converges to $1$ as the sample size increases to infinity, and a valid test must properly control the type $1$ error level. Modern datasets are often nonlinear, high-dimensional, and noisy, where density estimation and traditional statistical methods fail to be applicable. As multi-modal data are prevalent in much data-intensive research, a powerful, intuitive, and easy-to-use method for detecting general relationships is pivotal.

The classical Pearson’s correlation [@Pearson1895] is still extensively employed in statistics, machine learning, and real-world applications. It is an intuitive statistic that quantifies the linear association, a special but extremely important relationship. A recent surge of interests has been placed on using distance metrics and kernel transformations to achieve consistent independence testing against all dependencies. A notable example is the distance correlation ([[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}) [@SzekelyRizzoBakirov2007; @SzekelyRizzo2009; @SzekelyRizzo2013a; @SzekelyRizzo2014]: the population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is defined via the characteristic functions of the underlying random variables, while the sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} can be conveniently computed via the pairwise Euclidean distances of given observations. [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} enjoys universal consistency against any joint distribution of finite second moments, and is applicable to any metric space of strong negative type [@Lyons2013]. Notably, the idea of distance-based correlation measure can be traced back to the Mantel coefficient [@Mantel1967; @JosseHolmes2013]: the sample version differs from sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} only in centering, garnered popularity in ecology and biology applications, but does not have the consistency property of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}.

Developed almost in parallel from the machine learning community, the kernel-based method ([[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}) [@GrettonEtAl2005; @GrettonGyorfi2010] has a striking similarity with [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}: it is formulated by kernels instead of distances, can be estimated on sample data via the sample kernel matrix, and is universally consistent when using any characteristic kernel. Indeed, it is shown in [@SejdinovicEtAl2013] that there exists a mapping from kernel to metric (and vice versa) such that [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{} equals [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}. Another competitive method is the Heller-Heller-Gorfine method ([[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{}) [@HellerGorfine2013; @heller2016consistent]: it is also universally consistent by utilizing the rank information and the Pearson’s chi-square test, but has better finite-sample testing powers over [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in a collection of common nonlinear dependencies. There are other consistent methods available, such as the [[[<span style="font-variant:small-caps;">Copula</span>]{}]{}]{} method that tests independence based on the empirical copula process [@Genest06; @Genest07; @Holmes09], entropy-based methods [@Mendes06], and methods tailored for univariate data [@Reshef2011].

As the number of observations in many real world problems (e.g., genetics and biology) are often limited and very costly to increase, finite-sample testing power is crucial for certain data exploration tasks: [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} has been shown to perform well in monotone relationships, but not so well in nonlinear dependencies such as circles and parabolas; the performance of [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{} are often the opposite of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, which perform slightly inferior to [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in monotone relationships but excel in various nonlinear dependencies.

From another point of view, unraveling the nonlinear structure has been intensively studied in the manifold learning literature [@TenenbaumSilvaLangford2000; @SaulRoweis2000; @BelkinNiyogi2003]: by approximating a linear manifold locally via the k-nearest neighbors at each point, these nonlinear techniques can produce better embedding results than linear methods (like PCA) in nonlinear data. The main downside of manifold learning often lies in the parameter choice, i.e., the number of neighbor or the correct embedding dimension is often hard to estimate and requires cross-validation. Therefore, assuming a satisfactory neighborhood size can be efficiently determined in a given nonlinear relationship, the local correlation measure shall work better than the global correlation measure; and if the parameter selection is sufficiently adaptive, the optimal local correlation shall equal the global correlation in linear relationships.

In this manuscript we formalize the notion of population local distance correlations and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, explore their theoretical properties both asymptotically and in finite-sample, and propose an improved Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} algorithm. By combing distance correlation with the locality principle, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} inherits the universal consistency in testing, is able to efficiently search over all local scales and determine the optimal correlation, and enjoys the best testing powers throughout the simulations. A number of real data applications via [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} are pursued in [@ShenEtAl2016], e.g., testing brain images versus personality and disease, identify potential protein biomarkers for cancer, etc. And [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} are employed for vertex dependence testing and screening in [@mgc3; @mgc4].

The paper is organized as follows: In Section \[sec:main1\], we define the population local distance correlation and population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} via the characteristic functions of the underlying random variables and the nearest neighbor graphs, and show how the local variants are related to the distance correlation. In Section \[sec:main2\], we consider the sample local correlation on finite-samples, prove its convergence to the population version, and discuss the centering and ranking scheme. In Section \[sec:main3\], we present a thresholding-based algorithm for Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, prove its convergence property, propose a theoretically sound threshold choice, manifest that [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is valid and consistent under the permutation test, and finish the section with a number of fundamental properties for the local correlations and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}. The comprehensive simulations in Section \[sec:exp\] exhibits the empirical advantage of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, and the paper is concluded in Section \[sec:dis\]. All proofs are in Appendix \[sec:proofs\], the simulation functions are presented in Appendix \[appen:function\], and the code are available on Github 1[1]{}[[^4] ]{}0[1]{}and CRAN 1[1]{}[[^5]]{}0[1]{}.

Multiscale Graph Correlation for Random Variables {#sec:main1}
=================================================

Distance Correlation Review
---------------------------

We first review the original distance correlation in [@SzekelyRizzoBakirov2007]. A non-negative weight function $w(t,s)$ on $(t,s) \in \mathbb{R}^{p} \times \mathbb{R}^{q}$ is defined as: $$\begin{aligned}
        w(t,s) &=  (c_{p}c_{q} |t|^{1+p}|s|^{1+q})^{-1},\end{aligned}$$ where $c_{p}=\frac{\pi^{(1+p)/2}}{\Gamma((1+p)/2)}$ is a non-negative constant tied to the dimensionality $p$, and $\Gamma(\cdot)$ is the complete Gamma function. Then the population distance covariance, variance and correlation are defined by $$\begin{aligned}
dCov({\ensuremath{X}},{\ensuremath{Y}}) &=  \int_{\mathbb{R}^{p} \times \mathbb{R}^{q}} |E(g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s))-E(g_{{\ensuremath{X}}}(t))E(g_{{\ensuremath{Y}}}(s))|^{2} w(t, s)dtds, \\
dVar({\ensuremath{X}}) &= dCov({\ensuremath{X}},{\ensuremath{X}}), \\
dVar({\ensuremath{Y}}) &= dCov({\ensuremath{Y}},{\ensuremath{Y}}), \\
dCorr({\ensuremath{X}},{\ensuremath{Y}}) &= \frac{dCov({\ensuremath{X}},{\ensuremath{Y}})}{\sqrt{dVar({\ensuremath{X}}) \cdot dVar({\ensuremath{Y}})}},\end{aligned}$$ where $|\cdot|$ is the complex modulus, $g_{\cdot}(\cdot)$ denotes the exponential transformation within the expectation of the characteristic function, i.e., $g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s) = e^{\textbf{i} \left\langle t,{\ensuremath{X}}\right\rangle  +\textbf{i} \left\langle  s,{\ensuremath{Y}}\right\rangle }$ ($\textbf{i}$ represents the imaginary unit) and $E(g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s))$ is the characteristic function. Note that distance variance equals $0$ if and only if the random variable is a constant, in which case distance correlation shall be set to $0$. The main property of population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is the following.

For any two random variables $({\ensuremath{X}},{\ensuremath{Y}})$ with finite first moments, $dCorr({\ensuremath{X}},{\ensuremath{Y}})=0$ if and only if ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent.

To estimate the population version on sample data, the sample distance covariance is computed by double centering the pairwise Euclidean distance matrix of each data, followed by summing over the entry-wise product of the two centered distance matrices. When the underlying random variables have finite second moments, the sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is shown to converge to the population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} , and is thus universally consistent for testing independence against all joint distributions of finite second moments.

Population Local Correlations {#main1a}
-----------------------------

Next we formally define the population local distance covariance, variance, correlation by combining the k-nearest neighbor graphs with the distance covariance. For simplicity, they are named the local covariance, local variance, and local correlation from now on, and we always assume the following regularity conditions: $$\begin{aligned}
& 1) \mbox{ $({\ensuremath{X}},{\ensuremath{Y}})$ have finite second moments}, \\
& 2) \mbox{ Neither random variable is a constant}, \\
& 3) \mbox{ $({\ensuremath{X}},{\ensuremath{Y}})$ are continuous random variables}. \end{aligned}$$ The finite second moments assumption is required by [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, and also required by the local version to establish convergence and consistency. The non-constant condition is to avoid the trivial case and make sure population local correlations behave well. The continuous assumption is for ease of presentation, so the definition and related properties can be presented in a more elegant manner. Indeed, for any discrete random variable one can always apply jittering (i.e., add trivial white noise) to make it continuous without altering the independence testing.

Suppose $({\ensuremath{X}},{\ensuremath{Y}}), ({\ensuremath{X}}',{\ensuremath{Y}}'), ({\ensuremath{X}}'',{\ensuremath{Y}}''), ({\ensuremath{X}}''',{\ensuremath{Y}}''')$ are *iid* as $F_{{\ensuremath{X}}{\ensuremath{Y}}}$. Let ${\boldsymbol{I}}(\cdot)$ be the indicator function, define two random variables $$\begin{aligned}
{\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} &={\boldsymbol{I}}(\int_{B({\ensuremath{X}},\|{\ensuremath{X}}'-{\ensuremath{X}}\|)} dF_{\ensuremath{X}}(u) \leq \rho_k)  \\
{\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} &={\boldsymbol{I}}(\int_{B({\ensuremath{Y}}',\|{\ensuremath{Y}}'-{\ensuremath{Y}}\|)} dF_{\ensuremath{Y}}(v) \leq \rho_l)\end{aligned}$$ with respect to the closed balls $B({\ensuremath{X}},\|{\ensuremath{X}}'-{\ensuremath{X}}\|)$ and $B({\ensuremath{Y}}',\|{\ensuremath{Y}}-{\ensuremath{Y}}'\|)$ centered at ${\ensuremath{X}}$ and ${\ensuremath{Y}}'$ respectively. Then let $\overline{\cdot}$ denote the complex conjugate, define $$\begin{aligned}
h^{\rho_{k}}_{{\ensuremath{X}}}(t) &=(g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}'}(t)}-g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}''}(t)}) {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} \\
h^{\rho_{l}}_{{\ensuremath{Y}}'}(s) &=(g_{{\ensuremath{Y}}'}(s)\overline{g_{{\ensuremath{Y}}}(s)}-g_{{\ensuremath{Y}}'}(s)\overline{g_{{\ensuremath{Y}}'''}(s)}) {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}\end{aligned}$$ as functions of $t \in \mathbb{R}^{p}$ and $s \in \mathbb{R}^{q}$ respectively,

The population local covariance, variance, correlation at any $(\rho_k,\rho_l) \in [0,1] \times [0,1]$ are defined as $$\begin{aligned}
\label{eq:dcov1}
dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) &= \int_{\mathbb{R}^{p} \times \mathbb{R}^{q}} \{ E(h^{\rho_{k}}_{{\ensuremath{X}}}(t) \overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)})-E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))E(\overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)})\} w(t, s)dtds,\\
dVar^{\rho_k}({\ensuremath{X}}) &= dCov^{\rho_k, \rho_k}({\ensuremath{X}},{\ensuremath{X}}), \nonumber \\
dVar^{\rho_l}({\ensuremath{Y}}) &= dCov^{\rho_l, \rho_l}({\ensuremath{Y}},{\ensuremath{Y}}), \nonumber \\
dCorr^{\rho_k,\rho_l}({\ensuremath{X}},{\ensuremath{Y}}) &= \frac{dCov^{\rho_k,\rho_l}({\ensuremath{X}},{\ensuremath{Y}})}{\sqrt{dVar^{\rho_k}({\ensuremath{X}}) \cdot dVar^{\rho_l}({\ensuremath{Y}})}},\end{aligned}$$ where we limit the domain of population local correlation to $$\begin{aligned}
\mathcal{S}_{\epsilon}=\big\{(\rho_{k},\rho_{l}) \in [0,1] \times [0,1] \mbox{ that satisfies } \min\{dVar^{\rho_k}(X),dVar^{\rho_l}(Y)\} \geq \epsilon\big\}\end{aligned}$$ for a small positive $\epsilon$ that is no larger than $\min\{dVar(X),dVar(Y)\}$.

The domain of local correlation needs to be limited so the population version is well-behaved. For example, when ${\ensuremath{X}}$ is a constant or $\rho_{k}=0$, $dVar^{\rho_k}(X)$ equals $0$ and the corresponding local correlation is not well-defined. All subsequent analysis for the population local correlations is based on the domain $\mathcal{S}_{\epsilon}$, which is non-empty and compact as shown in Theorem \[thmMax\]. In practice, it suffices to set $\epsilon$ as any small positive number, see the sample version in Section \[sec:main2\]. Also note that in either indicator function, the two random variables and the distribution $dF$ are independent, e.g., at any realization $(x, x')$ of $({\ensuremath{X}},{\ensuremath{X}}')$, the first indicator equals ${\boldsymbol{I}}(\int_{B(x,\|x'-x\|)} dF_{\ensuremath{X}}(u) \leq \rho_k)$, and its expectation is taken with respect to $({\ensuremath{X}},{\ensuremath{X}}')$.

The above definition makes use of the characteristic functions, which is akin to the original definition of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} and easier to show consistency. Alternatively, the local covariance can be equivalently defined via the pairwise Euclidean distances. The alternative definition better motivates the sample version in Section \[sec:main2\], is often handy for understanding and proving theoretical properties, and suggests that local covariance is always a real number, which is not directly obvious from Equation \[eq:dcov1\].

\[thm1\] Suppose $({\ensuremath{X}},{\ensuremath{Y}}),({\ensuremath{X}}',{\ensuremath{Y}}'),({\ensuremath{X}}'',{\ensuremath{Y}}''),({\ensuremath{X}}''',{\ensuremath{Y}}''')$ are *iid* as $F_{{\ensuremath{X}}{\ensuremath{Y}}}$, and define $$\begin{aligned}
d^{\rho_{k}}_{{\ensuremath{X}}} &=(\| {\ensuremath{X}}-{\ensuremath{X}}' \| - \|{\ensuremath{X}}-{\ensuremath{X}}''\|) {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} \\
d^{\rho_{l}}_{{\ensuremath{Y}}'} &=(\| {\ensuremath{Y}}'-{\ensuremath{Y}}\| - \|{\ensuremath{Y}}'-{\ensuremath{Y}}'''\|) {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}\end{aligned}$$

The local covariance in Equation \[eq:dcov1\] can be equally defined as $$\begin{aligned}
\label{eq:dcov2}
dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) = E(d^{\rho_{k}}_{{\ensuremath{X}}} d^{\rho_{l}}_{{\ensuremath{Y}}'}) - E(d^{\rho_{k}}_{{\ensuremath{X}}}) E(d^{\rho_{l}}_{{\ensuremath{Y}}'}),\end{aligned}$$ which shows that local covariance, variance, correlation are always real numbers.

Each local covariance is essentially a local version of distance covariance that truncates large distances at each point in the support, where the neighborhood size is determined by $(\rho_{k},\rho_{l})$. In particular, distance correlation equals the local correlation at the maximal scale, which will ensure the consistency of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}.

\[thm2\] At any $(\rho_k,\rho_l)\in \mathcal{S}_{\epsilon}$, $dCov^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}}) = 0$ when ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent. Moreover, at $(\rho_k,\rho_l)=(1,1)$, $dCov^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}}) = dCov({\ensuremath{X}},{\ensuremath{Y}})$. They also hold for the correlations by replacing all the $dCov$ by $dCorr$.

Population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and Optimal Scale
------------------------------------------------------------------------------------------

The population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be naturally defined as the maximum local correlation within the domain, i.e., $$\begin{aligned}
\label{eq:pmgc}
{c}^{*}({\ensuremath{X}},{\ensuremath{Y}})=\max_{ (\rho_k,\rho_l) \in \mathcal{S}_{\epsilon}} \{dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})\},\end{aligned}$$ and the scale that attains the maximum is named the optimal scale $$\begin{aligned}
\label{eq:pmgc1}
(\rho_k,\rho_l)^{*}=\arg\max_{ (\rho_k,\rho_l) \in \mathcal{S}_{\epsilon}} \{dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})\}.\end{aligned}$$ The next theorem states the continuity of the local covariance, variance, correlation, and thus the existence of population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}.

\[thmMax\] Given two continuous random variables $({\ensuremath{X}},{\ensuremath{Y}})$,

(a)

:   The local covariance is a continuous function with respect to $(\rho_{k},\rho_{l}) \in [0,1]^2$, so is local variance in $[0,1]$ and local correlation in $\mathcal{S}_{\epsilon}$.

(b)

:   The set $\mathcal{S}_{\epsilon}$ is always non-empty unless either random variable is a constant.

(c)

:   Excluding the trivial case in (b), the set $\{dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}}), (\rho_k,\rho_l) \in \mathcal{S}_{\epsilon}\}$ is always non-empty and compact, so an optimal scale $(\rho_k,\rho_l)^{*}$ and ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$ exist.

Therefore, population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and the optimal scale exist, are distribution dependent, and may not be unique. Without loss of generality, the optimal scale is assumed unique for presentation purpose. The population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is always no smaller than [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in magnitude, and equals $0$ if and only if independence, a property inherited from [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}.

\[thm3\] When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})=dCorr({\ensuremath{X}},{\ensuremath{Y}})=0$; when ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are not independent, ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}}) \geq dCorr({\ensuremath{X}},{\ensuremath{Y}})>0$.

Sample Local Correlations {#sec:main2}
=========================

Sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} can be easily calculated via properly centering the Euclidean distance matrices, and is shown to converge to the population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} [@SzekelyRizzoBakirov2007; @SzekelyRizzo2013a; @SzekelyRizzo2014]. Similarly, we show that the sample local correlation can be calculated via the Euclidean distance matrices upon truncating large distances for each sample observation, and the sample version converges to the respective population local correlation.

Definition {#sec:defi}
----------

Given pairs of observations $(x_{i},y_{i}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{q}$ for $i=1,\ldots,n$, denote $\mathcal{X}_{n}=[x_{1},\ldots,x_{n}]$ as the data matrix with each column representing one sample observation, and similarly $\mathcal{Y}_{n}$. Let $\tilde{A}$ and $\tilde{B}$ be the $n \times n$ Euclidean distance matrices of $\mathcal{X}_{n}=\{x_{i}\}$ and $\mathcal{Y}_{n}=\{y_{i}\}$ respectively, i.e., $\tilde{A}_{ij}=\|x_{i}-x_{j}\|$. Then we compute two column-centered matrices $A$ and $B$ with the diagonals excluded, i.e., $\tilde{A}$ and $\tilde{B}$ are centered within each column such that $$\label{localCoef2}
    A_{ij}=
    \begin{cases}
      \tilde{A}_{ij}-\frac{1}{n-1}\sum_{s=1}^{n} \tilde{A}_{sj}, & \text{if $i \neq j$}, \\    
      0, & \text{if $i=j$};
    \end{cases} \qquad \qquad
    B_{ij}=
    \begin{cases}
      \tilde{B}_{ij}-\frac{1}{n-1}\sum_{s=1}^{n} \tilde{B}_{sj}, & \text{if $i \neq j$}, \\    
      0, & \text{if $i=j$};
    \end{cases}$$

Next we define $\{R^{A}_{ij}\}$ as the “rank” of $x_i$ relative to $x_j$, that is, $R^{A}_{ij}=k$ if $x_i$ is the $k^{th}$ closest point (or “neighbor”) to $x_j$, as determined by ranking the set $\{\tilde{A}_{1j},\tilde{A}_{2j},\ldots,\tilde{A}_{nj}\}$ by ascending order. Similarly define $R^{B}_{ij}$ for the $y$’s. As we assumed $({\ensuremath{X}},{\ensuremath{Y}})$ are continuous, with probability $1$ there is no repeating observation and the ranks always take value in $\{1,\ldots,n\}$. In practice ties may occur, and we recommend either using minimal rank to keep the ties or jittering to break the ties, which is discussed at the end of this section.

For any $(k,l) \in [n]^2=\{1,\ldots,n\} \times \{1,\ldots,n\}$, we define the rank truncated matrices $A^{k}$ and $B^{l}$ as $$\begin{aligned}
A_{ij}^{k} &=A_{ij} {\boldsymbol{I}}(R^{A}_{ij} \leq k), \\
B_{ij}^{l} &=B_{ij} {\boldsymbol{I}}(R^{B}_{ij} \leq l).\end{aligned}$$ Let $\circ$ denote the entry-wise product, ${\hat{E}}(\cdot)=\frac{1}{n(n-1)}\sum_{i \neq j}^{n} (\cdot)$ denote the diagonal-excluded sample mean of a square matrix, then the sample local covariance, variance, and correlation are defined as: $$\begin{aligned}
dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) &= {\hat{E}}(A^{k} \circ B^{l'})- {\hat{E}}(A^{k}){\hat{E}}(B^{l}),\\
dVar^{k}(\mathcal{X}_{n}) &={\hat{E}}(A^{k} \circ A^{k'})- {\hat{E}}^2(A^{k}), \\
dVar^{l}(\mathcal{Y}_{n}) &={\hat{E}}(B^{l} \circ B^{l'})- {\hat{E}}^2(B^{l}), \\
dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) &=dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) / \sqrt{dVar^{k}(\mathcal{X}_{n}) \cdot dVar^{l}(\mathcal{Y}_{n})}.\end{aligned}$$ If either local variance is smaller than a preset $\epsilon > 0$ (e.g., the smallest positive local variance among all), then we set the corresponding $dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})=0$ instead. Note that once the rank is known, sample local correlations can be iteratively computed in $\mathcal{O}(n^2)$ rather than a naive implementation of $\mathcal{O}(n^3)$. A detailed running time comparison is presented in Section \[sec:exp\].

In case of ties, minimal rank offers a consecutive indexing of sample local correlations, e.g., if ${\ensuremath{Y}}$ only takes two values, $R^{B}_{ij}$ takes value in $\{1,2\}$ under minimal rank, but maximal rank yields $\{\frac{n}{2},n\}$. The sample local correlations are not affected by the tie scheme, but minimal rank is more convenient to work with for implementation purposes. Alternatively, one can break ties deterministically or randomly, e.g., apply jittering to break all ties. For example, in the Bernoulli relationship of Figure \[f:dependencies\], there are only three points for computing sample local correlations and the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equals $0.9$. If white noise of variance $0.01$ were added to the data, we break all ties and obtain a much larger number of sample local correlations. The resulting Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is $0.8$, which is slightly smaller but still much larger than $0$ and implies a strong dependency.

Whether the random variable is continuous or discrete, and whether the ties in sample data are broken or not, does not affect the theoretical results except in certain theorem statements. For example, in Theorem \[thm4\], the convergence still holds for discrete random variables, but the index pair $(k,l)$ does not necessarily correspond to the population version at $(\rho_{k},\rho_{l})=(\frac{k-1}{n-1}, \frac{l-1}{n-1})$, e.g., when ${\ensuremath{X}}$ is Bernoulli with probability $0.8$ and minimal rank is used, $k=1$ corresponds to $\rho_{k}=0.8$ instead of $\rho_{k}=\frac{k-1}{n-1}$. Nevertheless, Theorem \[thm4\] and all results in the paper hold regardless of continuous or discrete random variables, but the presentation is more elegant for the continuous case.

Convergence Property
--------------------

The sample local covariance, variance, correlation are designed to converge to the respective population versions. Moreover, the expectation of sample local covariance equals the population counterpart up to a difference of $\mathcal{O}(\frac{1}{n})$, and the variance diminishes at the rate of $\mathcal{O}(\frac{1}{n})$.

\[thm4\] Suppose each column of $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}$ are generated *iid* from $({\ensuremath{X}},{\ensuremath{Y}}) \sim F_{{\ensuremath{X}}{\ensuremath{Y}}}$. The sample local covariance satisfies $$\begin{aligned}
E(dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})) &= dCov^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}}) +\mathcal{O}(1/n) \\
Var(dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})) &= \mathcal{O}(1/n)\\
dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) &\stackrel{n \rightarrow \infty}{\rightarrow} dCov^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}}),\end{aligned}$$ where $\rho_{k}=\frac{k-1}{n-1}$ and $\rho_{l}=\frac{l-1}{n-1}$. In particular, the convergence is uniform and also holds for the local correlation, i.e., for any $\epsilon$ there exists $n_{\epsilon}$ such that for all $n > n_{\epsilon}$, $$\begin{aligned}
|dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) -dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})|< \epsilon\end{aligned}$$ for any pair of $(\rho_{k},\rho_{l}) \in \mathcal{S}_{\epsilon}$.

The convergence property ensures that Theorem \[thm2\] holds asymptotically for the sample version.

\[thm5\] For any $(k,l)$, $dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow 0$ when ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent. In particular, $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow dCorr({\ensuremath{X}},{\ensuremath{Y}})$.

Moreover, one can show that $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n}) \approx dCorr(\mathcal{X}_{n},\mathcal{Y}_{n})$ for the unbiased sample distance correlation in [@SzekelyRizzo2014] up-to a small difference of $\mathcal{O}(\frac{1}{n})$, which can be verified by comparing Equation \[localCoef2\] to Equation 3.1 in [@SzekelyRizzo2014].

Centering and Ranking
---------------------

To combine distance testing with the locality principle, other than the procedure proposed in Equation \[eq:dcov2\], there are a number of alternative options to center and rank the distance matrices. For example, letting $$\begin{aligned}
d^{\rho_{k}}_{{\ensuremath{X}}} &=(\| {\ensuremath{X}}-{\ensuremath{X}}' \| - \|{\ensuremath{X}}-{\ensuremath{X}}''\|-\|{\ensuremath{X}}'-{\ensuremath{X}}''\|+\|{\ensuremath{X}}''-{\ensuremath{X}}'''\|) {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}, \\
d^{\rho_{l}}_{{\ensuremath{Y}}'} &=(\| {\ensuremath{Y}}'-{\ensuremath{Y}}\| - \|{\ensuremath{Y}}'-{\ensuremath{Y}}''\|-\|{\ensuremath{Y}}-{\ensuremath{Y}}''\|+\|{\ensuremath{Y}}''-{\ensuremath{Y}}'''\|) {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}\end{aligned}$$ still guarantees the resulting local correlation at maximal scale equals the distance correlation; and letting $$\begin{aligned}
d^{\rho_{k}}_{{\ensuremath{X}}} &=\| {\ensuremath{X}}-{\ensuremath{X}}' \|  {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}, \\
d^{\rho_{l}}_{{\ensuremath{Y}}'} &=\| {\ensuremath{Y}}'-{\ensuremath{Y}}\|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}},\end{aligned}$$ makes the resulting local correlation at maximal scale equal the [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} coefficient, the earliest distance-based correlation coefficient.

Nevertheless, the centering and ranking strategy proposed in Equation \[eq:dcov2\] is more faithful to k-nearest neighbor graph: the indicator ${\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}$ equals $1$ if and only if $\int_{B({\ensuremath{X}},\|{\ensuremath{X}}'-{\ensuremath{X}}\|)} dF_{\ensuremath{X}}(u) \leq \rho_k$, which happens with probability $\rho_{k}$. Viewed another way, when conditioned on $({\ensuremath{X}},{\ensuremath{X}}')=(x,x')$, the indicator equals $1$ if and only if $Prob(\|x'-x\|<\|{\ensuremath{X}}''-x\|) \leq \rho_{k}$, thus matching the column ranking scheme in Equation \[localCoef2\]. Indeed, the locality principle used in [@TenenbaumSilvaLangford2000; @SaulRoweis2000; @BelkinNiyogi2003] considers the k-nearest neighbors of each sample point in local computation, an essential step to yield better nonlinear embeddings. On the centering side, the [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} test appears to be an attractive option due to its simplicity in centering. All the [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{} have their theoretical consistency, while the [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} coefficient does not, despite it being merely a different centering of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}. An investigation of the population form of [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} yields some additional insights:

Given $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}$, the [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} coefficient on sample data is computed as $$\begin{aligned}
M(\mathcal{X}_{n},\mathcal{Y}_{n}) &= {\hat{E}}(\tilde{A} \circ \tilde{B})-{\hat{E}}(\tilde{A}){\hat{E}}(\tilde{B}) \\
Mantel(\mathcal{X}_{n},\mathcal{Y}_{n}) &= \frac{M(\mathcal{X}_{n},\mathcal{Y}_{n})}{\sqrt{M(\mathcal{X}_{n},\mathcal{X}_{n}) M(\mathcal{Y}_{n},\mathcal{Y}_{n})}},\end{aligned}$$ where $\tilde{A}_{ij}$ and $\tilde{B}_{ij}$ are the pairwise Euclidean distance, and ${\hat{E}}(\cdot)=\frac{1}{n(n-1)}\sum_{i \neq j}^{n} (\cdot)$ is the diagonal-excluded sample mean of a square matrix.

\[cor1\] Suppose each column of $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}$ are *iid* as $F_{{\ensuremath{X}}{\ensuremath{Y}}}$, and $({\ensuremath{X}},{\ensuremath{Y}}), ({\ensuremath{X}}',{\ensuremath{Y}}')$ are also *iid* as $F_{{\ensuremath{X}}{\ensuremath{Y}}}$. Then $$\begin{aligned}
Mantel(\mathcal{X}_{n},\mathcal{Y}_{n}) &\rightarrow Mantel({\ensuremath{X}},{\ensuremath{Y}}) = \frac{M({\ensuremath{X}},{\ensuremath{Y}})}{\sqrt{M({\ensuremath{X}},{\ensuremath{X}}) M({\ensuremath{Y}},{\ensuremath{Y}})}},\end{aligned}$$ where $$\begin{aligned}
M({\ensuremath{X}},{\ensuremath{Y}}) &= \int_{\mathbb{R}^{p} \times \mathbb{R}^{q}} \{|E(g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s))|^2-|E(g_{{\ensuremath{X}}}(t))E(g_{{\ensuremath{Y}}}(s))|^{2}\} w(t, s)dtds \\
&= E(\| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|) - E(\|{\ensuremath{X}}-{\ensuremath{X}}'\|)E(\|{\ensuremath{Y}}-{\ensuremath{Y}}'\|)) \\
&= Cov(\| {\ensuremath{X}}-{\ensuremath{X}}' \|, \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|).\end{aligned}$$

Corollary \[cor1\] suggests that [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} is actually a two-sided test based on the absolute difference of characteristic functions: under certain dependency structure, the [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} coefficient can be negative and still imply dependency (i.e., $|E(g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s))| < |E(g_{{\ensuremath{X}}}(t))E(g_{{\ensuremath{Y}}}(s))|$); whereas population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} are always no smaller than $0$, and any negativity of the sample version does not imply dependency. Therefore, [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} is only appropriate as a two-sided test, which is evaluated in Section \[sec:exp\]. Another insight is that [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{}, unlike [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, is not universally consistent: due to the integral $w$, one can construct a joint distribution such that the population [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} equals $0$ under dependence (see Remark 3.13 in [@Lyons2013] for an example of dependent random variables with uncorrelated distances). However, empirically, simple centering is still effective in a number of common dependencies (like two parabolas and diamond in Figure \[f:1DAll\]).

Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and Estimated Optimal Scale {#sec:main3}
================================================================================================

A naive sample version of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be defined as the maximum of all sample local correlations $$\begin{aligned}
\max_{(k,l) \in [n]^2}\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\}.\end{aligned}$$ Although the convergence to population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be guaranteed, the sample maximum is a biased estimator of the population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} in Equation \[eq:pmgc\]. For example, under independence, population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equals $0$, while the maximum sample local correlation has expectation larger than $0$, which may negate the advantage of searching locally and hurt the testing power.

This motivates us to compute Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} as a smoothed maximum within the largest connected region of thresholded local correlations. The purpose is to mitigate the bias of a direct maximum, while maintaining its advantage over [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in the test statistic. The idea is that in case of dependence, local correlations on the grid near the optimal scale shall all have large correlations; while in case of independence, a few local correlations may happen to be large, but most nearby local correlations shall still be small. The idea can be similarly adapted whenever there are multiple correlated test statistics or multiple models available, for which taking a direct maximum may yield too much bias [@mgc3]. From another perspective, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is like taking a regularized maximum.

Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}
--------------------------------------------------------------------

The procedure is as follows:

Input: 

:   A pair of datasets $(\mathcal{X}_{n}, \mathcal{Y}_{n})$.

Compute the Local Correlation Map:

:   Compute all local correlations:\
    $\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}), (k,l) \in [n]^2\}$.

Thresholding: 

:   Pick a threshold $\tau_n \geq 0$, denote $LC(\cdot)$ as the operation of taking the largest connected component, and compute the largest region $R$ of thresholded local correlations: $$\begin{aligned}
    \label{eq:region}
    &R=LC(\{(k,l) \mbox{ such that } dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})>\max\{\tau_n, dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n})\} \}). \end{aligned}$$ Within the region $R$, set $$\begin{aligned}
    \label{eq:smgc1}
    {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})&=\max_{ (k,l) \in R} \{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\}\\
    (k_{n},l_{n})^{*}&=\arg\max_{ (k,l) \in R} \{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\}\end{aligned}$$ as the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and the estimated optimal scale. If the number of elements in $R$ is less than $2n$, or the above thresholded maximum is no more than $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n})$, we instead set ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})=dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n})$ and $(k_{n},l_{n})^{*}=(n,n)$.

Output: 

:   Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} ${c}^{*}(\mathcal{X}_{n}, \mathcal{Y}_{n})$ and the estimated optimal scale $(k_{n},l_{n})^{*}$.

If there are multiple largest regions, e.g., $R_{1}$ and $R_{2}$ where their number of elements are more than $2n$ and coincide with each other, then it suffices to let $R = R_{1} \displaystyle \cup R_{2}$ and locate the [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} statistic within the union. The selection of at least $2n$ elements for $R$ is an empirical choice, which balances the bias-variance trade-off well in practice. The parameter can be any positive integer without affecting the validity and consistency of the test. But if the parameter is too large, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} tends to be more conservative and is unable to detect signals in strongly nonlinear relationships (e.g., trigonometric functions), and performs closer and closer to [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}; if the parameter is set to a very small fixed number, the bias is inflated so [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} tends to perform similarly as directly maximizing all local correlations.

Convergence and Consistency
---------------------------

The proposed Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is algorithmically enforced to be no less than the local correlation at the maximal scale, and also no more than the maximum local correlation. It also ensures in Theorem \[thm3\] to hold for the sample version.

\[thm6\] Regardless of the threshold $\tau_n$, the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} statistic ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})$ satisfies

(a)

:   It always holds that $$\begin{aligned}
    \max_{(k,l) \in [n]^2}\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\} \geq {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}) \geq dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n}).\end{aligned}$$

(b)

:   When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow 0$; when ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are not independent, ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow$ a positive constant.

The next theorem states that if the threshold $\tau_n$ converges to $0$, then whenever population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is larger than population [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is also larger than sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} asymptotically; otherwise if the threshold does not converge to $0$, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} may equal sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} despite of the first moment advantage in population. Moreover, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} indeed converges to population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} when the optimal scale is in the largest thresholded region $R$. The empirical advantage of Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is illustrated in Figure \[f:dependencies\].

\[thm7\] Suppose each column of $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}$ are *iid* as continuous $({\ensuremath{X}}, {\ensuremath{Y}}) \sim F_{{\ensuremath{X}}{\ensuremath{Y}}}$, and the threshold choice $\tau_n \rightarrow 0$ as $n \rightarrow \infty$.

(a)

:   Assume that ${c}^{*}({\ensuremath{X}}, {\ensuremath{Y}}) > Dcorr({\ensuremath{X}}, {\ensuremath{Y}})$ under the joint distribution. Then ${c}^{*}(\mathcal{X}_{n}, \mathcal{Y}_{n}) > Dcorr(\mathcal{X}_{n}, \mathcal{Y}_{n})$ for $n$ sufficiently large.

(b)

:   Assume there exists an element within the the largest connected area of $\{(\rho_k,\rho_l) \in \mathcal{S}_{\epsilon}$ with $dCorr^{\rho_k,\rho_l}({\ensuremath{X}},{\ensuremath{Y}})> dCorr({\ensuremath{X}},{\ensuremath{Y}}) \}$, such that the the local correlation of that element equals ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$. Then ${c}^{*}(\mathcal{X}_{n}, \mathcal{Y}_{n}) \rightarrow {c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$.

Alternatively, Theorem \[thm7\](b) can be stated that the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} always converges to the maximal population local correlation within the largest connected area of thresholded local correlations. Therefore, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} converges either to [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} (when the area is empty) or something larger, thus improving over [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} statistic in first moment.

Choice of Threshold
-------------------

The choice of threshold $\tau_n$ is imperative for Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} to enjoy a good finite-sample performance, especially at small sample size. According to Theorem \[thm7\], the threshold shall converge to $0$ for Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} to prevail sample [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}.

A model-free threshold $\tau_n$ was previously used in [@ShenEtAl2016]: for the following set $$\begin{aligned}
\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) \mbox{ s.t. }  dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})<0\},\end{aligned}$$ let $\sigma^2$ be the sum of all its elements squared, and set $\tau_n=5\sigma$ as the threshold; if there is no negative local correlation and the set is empty, use $\tau_n=0.05$. Although the previous threshold is a data-adaptive choice that works pretty well empirically and does not affect the consistency of Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} in Theorem \[thm8\], it does not converge to $0$. The following finite-sample theorem from [@SzekelyRizzo2013a] motivates an improved threshold choice here:

Under independence of $({\ensuremath{X}},{\ensuremath{Y}})$, assume the dimensions of ${\ensuremath{X}}$ are exchangeable with finite variance, and so are the dimensions of ${\ensuremath{Y}}$. Then for any $n \geq 4$ and $v=\frac{n(n-3)}{2}$, as $p,q$ increase the limiting distribution of $(dCorr^{n,n}(\mathcal{X}_{n}, \mathcal{Y}_{n})+1)/2$ equals the symmetric Beta distribution with shape parameter $\frac{v-1}{2}$.

The above theorem leads to the new threshold choice:

\[cor2\] Denote $v=\frac{n(n-3)}{2}$, $z \sim Beta(\frac{v-1}{2})$, $F^{-1}_{z}(\cdot)$ as the inverse cumulative distribution function. The threshold choice $$\begin{aligned}
\tau_n= 2F^{-1}_{z} \Big(1-\frac{0.02}{n}\Big)-1 \end{aligned}$$ converges to $0$ as $n \rightarrow \infty$.

The limiting null distribution of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is still a good approximation even when $p,q$ are not large, thus provides a reliable bound for eliminating local correlations that are larger than [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} by chance or by noise. The intuition is that Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is mostly useful when it is much larger than [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in magnitude, which is often the case in non-monotone relationships as shown in Section \[sec:exp\] Figure \[f:dependencies\]. Alternatively, directly setting $\tau_{n}=0$ also guarantees the theoretical properties and works equally well when the sample size $n$ is moderately large.

Permutation Test {#sec:permutation}
----------------

To test independence on a pair of sample data $(\mathcal{X}_{n},\mathcal{Y}_{n})$, the random permutation test has been the popular choice [@GoodPermutationBook] for almost all methods introduced, as the null distribution of the test statistic can be easily approximated by randomly permuting one data set. We discuss the computation procedure, prove the testing consistency of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, and analyze the running time.

To compute the p-value of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} from the permutation test, first compute the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} statistic ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})$ on the observed data pair. Then the [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} statistic is repeatedly computed on the permuted data pair, e.g. $\mathcal{Y}_{n}=[y_{1},\ldots,y_{n}]$ is permuted into $\mathcal{Y}_{n}^{\pi}=[y_{\pi(1)},\ldots,y_{\pi(n)}]$ for a random permutation $\pi$ of size $n$, and compute ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi})$. The permutation procedure is repeated for $r$ times to estimate the probability $Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}))$, and the estimated probability is taken as the p-value of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}. The independence hypothesis is rejected if the p-value is smaller than a pre-set critical level, say $0.05$ or $0.01$. The following theorem states that [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} via the permutation test is consistent and valid.

\[thm8\] Suppose each column of $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}$ are generated *iid* from $F_{{\ensuremath{X}}{\ensuremath{Y}}}$. At any type $1$ error level $\alpha>0$, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is a valid test statistic that is consistent against all possible alternatives under the permutation test.

Miscellaneous Properties {#sec:misc}
------------------------

In this subsection, we first show a useful lemma expressing sample local covariance in Section \[sec:defi\] by matrix trace and eigenvalues, then list a number of fundamental and desirable properties for the local variance, local correlation, and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, akin to these of Pearson’s correlation and distance correlation as shown in [@SzekelyRizzoBakirov2007; @SzekelyRizzo2009].

\[lem1\] Denote $tr(\cdot)$ as the matrix trace, $\lambda_{i} [\cdot]$ as the $i$th eigenvalue of a matrix, and $J$ as the matrix of ones of size $n$. Then the sample covariance equals $$\begin{aligned}
dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) &= tr (A^{k}B^{l})- tr (A^{k}J)tr(B^{l}J)\\
&=tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)] \\
& = \sum_{i=1}^{n} \lambda_{i} [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)].\end{aligned}$$

\[thm:dvar\] For any random variable ${\ensuremath{X}}\sim F_{{\ensuremath{X}}} \in \mathbb{R}^{p}$, and any $\mathcal{X}_{n} \in \mathbb{R}^{p \times n}$ with each column *iid* as $F_{{\ensuremath{X}}}$,

(a)

:   Population and sample local variances are always non-negative, i.e., $$\begin{aligned}
    dVar^{\rho_{k}}({\ensuremath{X}}) \geq 0\\
    dVar^{k}(\mathcal{X}_{n}) \geq 0\end{aligned}$$ at any $\rho_{k} \in [0,1]$ and any $k \in [n]$.

(b)

:   $dVar^{\rho_{k}}({\ensuremath{X}}) =0$ if and only if either $\rho_k =0$ or $F_{{\ensuremath{X}}}$ is a degenerate distribution;

    $dVar^{k}(\mathcal{X}_{n}) =0$ if and only if either $k=1$ or $F_{{\ensuremath{X}}}$ is a degenerate distribution.

(c)

:   For two constants $v \in \mathbb{R}^{p},u \in \mathbb{R}$, and an orthonormal matrix $Q \in \mathbb{R}^{p \times p}$, $$\begin{aligned}
    dVar^{\rho_{k}}(v+uQ{\ensuremath{X}}) &=u^2 \cdot dVar^{\rho_{k}}({\ensuremath{X}})\\
    dVar^{k}(v^{T} J +u\mathcal{X}_{n} Q) &=u^2 \cdot dVar^{k}(\mathcal{X}_{n} ).\end{aligned}$$

Therefore, the local variances end up having properties similar to the distance variance in [@SzekelyRizzoBakirov2007], except the distance variance definition there takes a square root.

\[thm:dcor\] For any pair of random variable $({\ensuremath{X}},{\ensuremath{Y}}) \sim F_{{\ensuremath{X}}{\ensuremath{Y}}} \in \mathbb{R}^{p} \times \mathbb{R}^{q}$, and any $(\mathcal{X}_{n},\mathcal{Y}_{n}) \in \mathbb{R}^{p \times n} \times \mathbb{R}^{q \times n}$ with each column *iid* as $F_{{\ensuremath{X}}{\ensuremath{Y}}}$,

(a)

:   Symmetric and Boundedness: $$\begin{aligned}
    dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})=dCorr^{\rho_{l},\rho_{k}}({\ensuremath{Y}},{\ensuremath{X}}) \in [-1,1]\\
    dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})=dCorr^{l,k}(\mathcal{Y}_{n},\mathcal{X}_{n}) \in [-1,1]\end{aligned}$$ at any $(\rho_k,\rho_l) \in (0,1]^2$ and any $(k,l) \in [2,\ldots,n]^2$.

(b)

:   Assume $F_{\ensuremath{X}}$ is non-degenerate. Then at any $\rho_{k} > 0$, $dCorr^{\rho_{k},\rho_k}({\ensuremath{X}},{\ensuremath{Y}})=1$ if and only if $({\ensuremath{X}}, u {\ensuremath{Y}})$ are dependent via an isometry for some non-zero constant $u \in \mathbb{R}$.

    Assume $F_{\ensuremath{X}}$ is non-degenerate. Then at any $k > 1$, $dCorr^{k,k}(\mathcal{X}_{n},\mathcal{Y}_{n})=1$ if and only if $({\ensuremath{X}}, u {\ensuremath{Y}})$ are dependent via an isometry for some non-zero constant $u \in \mathbb{R}$.

(c)

:   Both population and Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} are symmetric and bounded: $$\begin{aligned}
    {c}^{*}({\ensuremath{X}},{\ensuremath{Y}})={c}^{*}({\ensuremath{Y}},{\ensuremath{X}}) \in [-1,1] \\
    {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})={c}^{*}(\mathcal{Y}_{n},\mathcal{X}_{n}) \in [-1,1].\end{aligned}$$

(d)

:   Assume $F_{\ensuremath{X}}$ is non-degenerate. Then ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})=1$ if and only if $({\ensuremath{X}}, u {\ensuremath{Y}})$ are dependent via an isometry for some non-zero constant $u \in \mathbb{R}$.

    Assume $F_{\ensuremath{X}}$ is non-degenerate. Then ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})=1$ if and only if $({\ensuremath{X}}, u {\ensuremath{Y}})$ are dependent via an isometry for some non-zero constant $u \in \mathbb{R}$.

The proof of Theorem \[thm:dcor\](b)(d) also shows that the local correlations and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} cannot be $-1$.

Experiments {#sec:exp}
===========

In the experiments, we compare Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} with [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{}, and [[[<span style="font-variant:small-caps;">Copula</span>]{}]{}]{} test on $20$ different simulation settings based on a combination of simulations used in previous works [@SzekelyRizzoBakirov2007; @SimonTibshirani2012; @GorfineHellerHeller2012]. Among the $20$ settings, the first $5$ are monotonic relationships (and several of them are linear or nearly so), the last simulation is an independent relationship, and the remaining settings consist of common non-monotonic and strongly nonlinear relationships. The exact distributions are shown in Appendix.

The Sample Statistics {#the-sample-statistics .unnumbered}
---------------------

Figure \[f:dependencies\] shows the sample statistics of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, and [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{} for each of the $20$ simulations in a univariate setting. For each simulation, we generate sample data $(\mathcal{X}_{n},\mathcal{Y}_{n})$ at $p=q=1$ and $n=100$ without any noise, then compute the sample statistics. From type $1-5$, the test statistics for both [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} are remarkably greater than $0$ and almost identical to each other. For the nonlinear relationships (type $6-19$), [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} benefits from searching locally and achieves a larger test statistic than [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}’s, which can be very small in these nonlinear relationships. For type $20$, the test statistics for both [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} are almost $0$ as expected. On the other hand, [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{}’s test statistic is large whenever there exists certain linear association, and almost $0$ otherwise. The comparison of sample statistics indicate that [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} may have inferior finite-sample testing power in nonlinear relationships, but a strong dependency signal is actually hidden in a local structure that [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} may recover.

Finite-Sample Testing Power {#finite-sample-testing-power .unnumbered}
---------------------------

Figure \[f:noise\] shows the finite-sample testing power of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, and [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{} for a linear and a quadratic relationship at $n=20$ and $p=q=1$ with white noise (controlled by a constant). The testing power of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is estimated as follows: we first generate dependent sample data $(\mathcal{X}_{n},\mathcal{Y}_{n})$ for $r=10,000$ replicates, compute Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} for each replicate to estimate the alternative distribution of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}. Then we generate independent sample data $(\mathcal{X}_{n},\mathcal{Y}_{n})$ using the same marginal distributions for $r=10,000$ replicates, compute Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} to estimate the null distribution, and estimate the testing power at type $1$ error level $\alpha=0.05$. The testing power of [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is estimated in the same manner, while the testing power of [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{} is directly computed via the t-test. [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} has the best power in the quadratic relationship, while being almost identical to [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Pearson</span>]{}]{}]{} in the linear relationship.

The same phenomenon holds throughout all the simulations we considered, i.e., [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} achieves almost the same power as [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} in monotonic relationships, while being able to improve the power in monotonic and strongly nonlinear relationships. The testing power of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} versus all other methods are shown in Figure \[f:1DAll\] for the univariate settings, and we plot the power versus the sample size from $5$ to $100$ for each simulation. Note that the noise level is tuned for each dependency for illustration purposes.

Figure \[f:nDAll\] compares the testing performance for the same $20$ simulations with a fixed sample size $n=100$ and increasing dimensionality. The relative powers in the univariate and multivariate settings are then summarized in Figure \[f:Summary2\]. [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is overall the most powerful method, followed by [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}. Since non-monotone relationships are prevalent among the $20$ settings, it is not a surprise that [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} is overall worse than [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}, both of which also excel at nonlinear relationships.

Note that the same $20$ simulations were also used in [@ShenEtAl2016] for evaluation purposes. The main difference is that the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} algorithm is now based on the improved threshold with theoretical guarantee. Comparing to the previous algorithm, the new threshold slightly improves the testing power in monotonic relationships (the first $5$ simulations).

Running Time {#sec:time .unnumbered}
------------

Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be computed and tested in the same running time complexity as distance correlation: Assume $p$ is the maximum feature dimension of the two datasets, distance computation and centering takes ${\mathcal{O}}(n^2 p)$, the ranking process takes ${\mathcal{O}}(n^2 \log n)$, all local covariances and correlations can be incrementally computed in $O(n^2)$ (the pseudo-code is shown in [@ShenEtAl2016]), the thresholding step of Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} takes $O(n^2)$ as well. Overall, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be computed in ${\mathcal{O}}(n^2 \max\{\log n,p\})$. In comparison, the [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{} statistic requires the same complexity as [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, while distance correlation saves on the $\log n$ term.

As the only part of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} that has the additional $\log n$ term is the column-wise ranking process, a multi-core architecture can reduce the running time to ${\mathcal{O}}(n^2 \max\{\log n,p\}/T)$. By making $T=\log(n)$ ($T$ is no more than $30$ at $1$ billion samples), [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} effectively runs in ${\mathcal{O}}(n^2 p)$ and is of the same complexity as [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}. The permutation test multiplies another $r$ to all terms except the distance computation, so overall the [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} testing procedure requires ${\mathcal{O}}(n^2 \max\{r,p\})$, which is the same as [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, [[[<span style="font-variant:small-caps;">Hhg</span>]{}]{}]{}, and [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}. Figure \[f:time\] shows that [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} has approximately the same complexity as [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{}, and is slower by a constant in the actual running time.

Conclusion {#sec:dis}
==========

In this paper, we formalize the population version of local correlation and [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, connect them to the sample counterparts, prove the convergence and almost unbiasedness from the sample version to the population version, as well as a number of desirable properties for a well-defined correlation measure. In particular, population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equals $0$ and the sample version converges to $0$ if and only if independence, making Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} valid and consistent under the permutation test. Moreover, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is designed in a computationally efficient manner, and the new threshold choice achieves both theoretical and empirical improvements. The numerical experiments confirm the empirical advantages of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} in a wide range of linear, nonlinear, high-dimensional dependencies.

There are many potential future avenues to pursue. Theoretically, proving when and how one method dominates another in testing power is highly desirable. As the methods in comparison have distinct formulations and different properties, it is often difficult to compare them directly. However, a relative efficiency analysis may be viable when limited to methods of similar properties, such as [[[<span style="font-variant:small-caps;">Dcorr</span>]{}]{}]{} and [[[<span style="font-variant:small-caps;">Hsic</span>]{}]{}]{}, or local statistic and global statistic. In terms of the locality principle, the geometric meaning of the local scale in [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is intriguing — for example, does the family of local correlations fully characterize the joint distribution, and what is the relationship between the optimal local scale and the dependency geometry — answering these questions may lead to further improvement of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, and potentially make the family of local correlations a valuable tool beyond testing.

Method-wise, there are a number of alternative implementations that may be pursued. For example, the sample local correlations can be defined via $\epsilon$ ball instead of nearest neighbor graphs, i.e., truncate large distances based on absolute magnitude instead of the nearest neighbor graph. The maximization and thresholding mechanism may be further improved, e.g., thresholding based on the covariance instead of correlation, or design a better regularization scheme. There are many alternative approaches that can maintain consistency in this framework, and it will be interesting to investigate a better algorithm. In particular, we name our method as “multiscale graph correlation" because the local correlations are computed via the k-nearest neighbor graphs, which is one way to generalize the distance correlation.

Application-wise, the [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} method can directly facilitate new discoveries in many kinds of scientific fields, especially data of limited sample size and high-dimensionality such as in neuroscience and omics [@ShenEtAl2016]. Within the domain of statistics and machine learning, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} can be a very competitive candidate in any methodology that requires a well-defined dependency measure, e.g., variable selection [@LiZhongZhu2012], time series [@Zhou2012], etc. Moreover, the very idea of locality may improve other types of distance-based tests, such as the energy distance for K-sample testing [@SzekelyRizzo2013b].

1

[1]{}

Acknowledgment {#acknowledgment .unnumbered}
==============

This work was partially supported by the National Science Foundation award DMS-1712947, and the Defense Advanced Research Projects Agency’s (DARPA) SIMPLEX program through SPAWAR contract N66001-15-C-4041. The authors are grateful to the anonymous reviewers for the invaluable feedback leading to significant improvement of the manuscript, and thank Dr. Minh Tang and Dr. Shangsi Wang for useful discussions and suggestions.

[**APPENDIX**]{}

Proofs {#sec:proofs}
======

Theorem \[thm1\] {#theoremthm1 .unnumbered}
----------------

Equation \[eq:dcov1\] defines the local covariance as $$\begin{aligned}
dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) = \int_{\mathbb{R}^{p}\times \mathbb{R}^{q}} E(h^{\rho_{k}}_{{\ensuremath{X}}}(t) \overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)})-E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))E(h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)) w(t, s)dtds.\end{aligned}$$ Expanding the first integral term yields $$\begin{aligned}
&\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t) \overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)}) w(t, s)dtds \\
=&\ E(\int (g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}'}(t)}-g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}''}(t)}) (\overline{g_{{\ensuremath{Y}}'}(s)}g_{{\ensuremath{Y}}}(s)-\overline{g_{{\ensuremath{Y}}'}(s)}g_{{\ensuremath{Y}}'''}(s)) w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}) \\
=&\ E(\int g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s) \overline{g_{{\ensuremath{X}}' {\ensuremath{Y}}'}(t,s)} w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}) \\
& -E( \int g_{{\ensuremath{X}}{\ensuremath{Y}}}(t,s) \overline{g_{{\ensuremath{X}}''}(t)g_{{\ensuremath{Y}}'}(s)} w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}) \\
& - E( \int \overline{g_{{\ensuremath{X}}' {\ensuremath{Y}}'}(t,s)} g_{{\ensuremath{X}}}(t)g_{{\ensuremath{Y}}'''}(s) w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}})\\
& + E( \int g_{{\ensuremath{X}}}(t)g_{{\ensuremath{Y}}'''}(s) \overline{ g_{{\ensuremath{X}}''}(t) g_{{\ensuremath{Y}}'}(s) } w(t, s)dtds \cdot  {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}) \\
=&\ E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} )  - E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ) \\
& - E( \| {\ensuremath{X}}'-{\ensuremath{X}}\| \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ) + E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ) \\
=&\ E(d^{\rho_{k}}_{{\ensuremath{X}}} d^{\rho_{l}}_{{\ensuremath{Y}}'}).\end{aligned}$$ Every other step being routine, the third equality transforms the $w(t,s)$ integral to Euclidean distances via the same technique employed in Remark 1 and the proof of Theorem 8 in [@SzekelyRizzo2009]. Also note that all four expectations are finite. For example, the first expectation in the third equality is finite, because $\| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|$ is always non-negative, and $E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|)$ is non-negative and finite by the finite second moments assumption on [$X$]{} and [$Y$]{}, such that $$\begin{aligned}
0 \leq E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ) \leq E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|),\end{aligned}$$ which can be similarly established for the other three expectations.

The second integral term can be decomposed into $$\begin{aligned}
\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))E(h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)) w(t, s)dtds = \int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))w(t, s)dtds \cdot \int E(h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)) w(t, s)dtds,\end{aligned}$$ because the first expectation only has $t$ and the second expectation only has $s$, and $w(t,s)$ is a product of $t$ and $s$. Then $$\begin{aligned}
&\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t)) w(t, s)dtds = E(\int g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}'}(t)}-g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}''}(t)} w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) \\
=& E(\int g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}'}(t)}  w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) - E(\int g_{{\ensuremath{X}}}(t)\overline{g_{{\ensuremath{X}}''}(t)}  w(t, s)dtds \cdot {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) \\
=& E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}})  - E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) \\
=& E(d^{\rho_{k}}_{{\ensuremath{X}}}),\end{aligned}$$ where the two expectations involved are also finite. Similarly $\int E(\overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)}) w(t, s)dtds = E( \| {\ensuremath{Y}}'-{\ensuremath{Y}}\| {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}})  - E( \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \| {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}})=E(d^{\rho_{l}}_{{\ensuremath{Y}}'})$. Thus $$\begin{aligned}
\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))E(h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)) w(t, s)dtds &= E(d^{\rho_{k}}_{{\ensuremath{X}}})E(d^{\rho_{l}}_{{\ensuremath{Y}}'}).\end{aligned}$$

Combining the results verifies that Equation \[eq:dcov2\] equals Equation \[eq:dcov1\]. Moreover, as every term in Equation \[eq:dcov2\] is of real-value, local covariance, variance, correlation are all real numbers.

Theorem \[thm2\] {#theoremthm2 .unnumbered}
----------------

When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, $$\begin{aligned}
\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t)\overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)})w(t, s)dtds =\int E(h^{\rho_{k}}_{{\ensuremath{X}}}(t))E(\overline{h^{\rho_{l}}_{{\ensuremath{Y}}'}(s)})w(t, s)dtds, \end{aligned}$$ thus $dCov^{\rho_{k}, \rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})=0$ at any $(\rho_k,\rho_l)$. So is the local correlation at any $(\rho_k,\rho_l) \in \mathcal{S}_{\epsilon}$.

To show the local covariance at the maximal scale $(\rho_k,\rho_l)=(1,1)$ equals the distance covariance, we proceed via the alternative definition in Theorem \[thm1\]: $$\begin{aligned}
&dCov^{\rho_{k}=1, \rho_{l}=1}({\ensuremath{X}},{\ensuremath{Y}}) =E(d^{\rho_{k}}_{{\ensuremath{X}}} d^{\rho_{l}}_{{\ensuremath{Y}}'}) \\
=&\ E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| )  - E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|) \\
&- E( \| {\ensuremath{X}}'-{\ensuremath{X}}\| \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \| ) + E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| ) E(\| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \|  ) \\
=&\ E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| )  - E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \|) \\
&- E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}'' \| ) + E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| ) E(\| {\ensuremath{Y}}-{\ensuremath{Y}}'' \|  ) \\
=&\ dCov({\ensuremath{X}},{\ensuremath{Y}}),\end{aligned}$$ where the first equality follows by noting that $E(d^{\rho_{k}}_{{\ensuremath{X}}})=E(d^{\rho_{l}}_{{\ensuremath{Y}}'})=0$ at $\rho_{k}=\rho_{l}=1$, the second equality holds by switching the random variable notations within each expectation, and the last equality is the alternative definition of distance covariance in Theorem 8 of [@SzekelyRizzo2009]. It follows that $dVar^{\rho_{k}=1}({\ensuremath{X}})=dVar({\ensuremath{X}})$, $dVar^{\rho_{l}=1}({\ensuremath{Y}})=dVar({\ensuremath{Y}})$, and $dCorr^{\rho_{k}=1, \rho_{l}=1}({\ensuremath{X}},{\ensuremath{Y}})=dCorr({\ensuremath{X}},{\ensuremath{Y}})$.

Theorem \[thmMax\] {#theoremthmmax .unnumbered}
------------------

Given two continuous random variables $({\ensuremath{X}},{\ensuremath{Y}})$, we first illustrate the continuity of local covariance with respect to $\rho_{k}$ at fixed $\rho_{l}$: For any $\delta$ with the understanding that $\rho_{k} \pm \delta \in [0,1]$, we have $$\begin{aligned}
dCov^{\rho_k+\delta, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) -dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) = E((d^{\rho_{k}+\delta}_{{\ensuremath{X}}}-d^{\rho_{k}}_{{\ensuremath{X}}}) d^{\rho_{l}}_{{\ensuremath{Y}}'}) - E(d^{\rho_{k}+\delta}_{{\ensuremath{X}}}-d^{\rho_{k}}_{{\ensuremath{X}}}) E(d^{\rho_{l}}_{{\ensuremath{Y}}'}),\end{aligned}$$ where the expectation is taken with respect to all random variables inside, and $$\begin{aligned}
d^{\rho_{k}+\delta}_{{\ensuremath{X}}}&=(\|{\ensuremath{X}}-{\ensuremath{X}}'\|-\|{\ensuremath{X}}-{\ensuremath{X}}''\|) {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}+\delta}\\
d^{\rho_{k}}_{{\ensuremath{X}}}&=(\|{\ensuremath{X}}-{\ensuremath{X}}'\|-\|{\ensuremath{X}}-{\ensuremath{X}}''\|) {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}\end{aligned}$$ Then Cauchy-Schwarz and finite second moment of ${\ensuremath{X}}$ yield that $$\begin{aligned}
&\lim_{\delta \rightarrow 0} |E(d^{\rho_{k}+\delta}_{{\ensuremath{X}}}-d^{\rho_{k}}_{{\ensuremath{X}}})|^2 \\
\leq &\ E\{(\|{\ensuremath{X}}-{\ensuremath{X}}'\|-\|{\ensuremath{X}}-{\ensuremath{X}}''\|)^2\} \lim_{\delta \rightarrow 0} E(|{\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}+\delta}-{\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}|^2)\\
=&\ 0.\end{aligned}$$ Moreover, the finite second moment of [$Y$]{} guarantees finiteness of $E( d^{\rho_{l}}_{{\ensuremath{Y}}'})$ and $$\begin{aligned}
&\lim_{\delta \rightarrow 0} |E((d^{\rho_{k}+\delta}_{{\ensuremath{X}}}-d^{\rho_{k}}_{{\ensuremath{X}}}) d^{\rho_{l}}_{{\ensuremath{Y}}'})|^2 \\
\leq &\ E\{(\|{\ensuremath{X}}-{\ensuremath{X}}'\|-\|{\ensuremath{X}}-{\ensuremath{X}}''\|)^2 {d^{\rho_{l}}_{{\ensuremath{Y}}'}}^2\} \lim_{\delta \rightarrow 0} E(|{\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}+\delta}-{\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}|^2)\\
=&\ 0,\end{aligned}$$ which leads to the continuity of local covariance with respect to $\rho_{k}$: $$\begin{aligned}
\lim_{\delta \rightarrow 0} dCov^{\rho_k+\delta, \rho_l}({\ensuremath{X}},{\ensuremath{Y}}) - dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}})=0.\end{aligned}$$

The same holds for fixed $\rho_{k}$ such that $$\begin{aligned}
\lim_{\delta \rightarrow 0} dCov^{\rho_k, \rho_l+\delta}({\ensuremath{X}},{\ensuremath{Y}}) - dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}})=0.\end{aligned}$$ Applying the above yields that $$\begin{aligned}
&dCov^{\rho_k+\delta_{1}, \rho_l+\delta_{2}}({\ensuremath{X}},{\ensuremath{Y}}) -dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}})\\
= &\ dCov^{\rho_k+\delta_{1}, \rho_l+\delta_{2}}({\ensuremath{X}},{\ensuremath{Y}}) -dCov^{\rho_k, \rho_l+\delta_{2}}({\ensuremath{X}},{\ensuremath{Y}}) + dCov^{\rho_k, \rho_l+\delta_{2}}({\ensuremath{X}},{\ensuremath{Y}}) - dCov^{\rho_k, \rho_l}({\ensuremath{X}},{\ensuremath{Y}})\\
\rightarrow & \  0 \mbox { for any $\delta_{1}$ and $\delta_{2}$ satisfying $|\delta_{1}+\delta_{2}| \rightarrow 0$.}\end{aligned}$$ So the local covariance is continuous with respect to $(\rho_{k},\rho_{l}) \in [0,1]\times[0,1]$. The continuity of the local variance can be shown similarly, and it follows that the local correlation is continuous in $\mathcal{S}_{\epsilon}$.

At $\rho_{k}=1$, $dVar^{\rho_k}({\ensuremath{X}})=dVar({\ensuremath{X}}) \geq 0$ with equality if and only if ${\ensuremath{X}}$ is a constant, and $\mathcal{S}_{\epsilon}$ is empty in the trivial case. Otherwise by the continuity of local variance, for any $\epsilon < dVar({\ensuremath{X}})$ there exists $\epsilon_{k}$ such that for all $\rho_{k} \in [\epsilon_{k},1]$, $dVar^{\rho_k}({\ensuremath{X}}) \geq \epsilon$. Same for $dVar^{\rho_l}({\ensuremath{Y}})$, thus $\mathcal{S}_{\epsilon}$ is non-empty except when either random variable is a constant. It follows that the local correlation is continuous within the non-empty and compact domain $\mathcal{S}_{\epsilon}$, and extreme value theorem ensures the existence of population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} and the optimal scale.

Theorem \[thm3\] {#theoremthm3 .unnumbered}
----------------

By Theorem \[thm2\] and definition of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}, it holds that $$\begin{aligned}
{c}^{*}({\ensuremath{X}},{\ensuremath{Y}}) \geq dCorr^{\rho_{k}=\rho_{l}=1}({\ensuremath{X}},{\ensuremath{Y}})=dCorr({\ensuremath{X}},{\ensuremath{Y}}).\end{aligned}$$ When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, all local correlations are $0$ by Theorem \[thm2\], so ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})=0$ as well. When dependent, distance correlation is larger than $0$, and it follows that ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}}) \geq dCorr({\ensuremath{X}},{\ensuremath{Y}})>0$. Therefore, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equals $0$ if and only if independence, just like the distance correlation.

Theorem \[thm4\] {#theoremthm4 .unnumbered}
----------------

We prove this theorem by three steps: **(i)**, the expectation of the sample local covariance is shown to equal the population local covariance; **(ii)**, the variance of the sample statistic is of $\mathcal{O}(\frac{1}{n})$; **(iii)**, sample local covariance is shown to convergence to the population counterpart uniformly. Then the convergence trivially extends to the sample local variance and correlation.

**(i)**: Expanding the first and second term of population local covariance in Equation \[eq:dcov2\], we have $E(d^{\rho_{k}}_{{\ensuremath{X}}} d^{\rho_{l}}_{{\ensuremath{Y}}'})=\alpha_{1}-\alpha_{2}-\alpha_{3}+\alpha_{4}$ with $$\begin{aligned}
\alpha_{1}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ),\\
\alpha_{2}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| \| {\ensuremath{Y}}-{\ensuremath{Y}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ),\\
\alpha_{3}&=E( \| {\ensuremath{X}}'-{\ensuremath{X}}\| \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}} ),\\
\alpha_{4}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| \| {\ensuremath{Y}}'-{\ensuremath{Y}}''' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}} {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}),\end{aligned}$$ and $E(d^{\rho_{k}}_{{\ensuremath{X}}})E(d^{\rho_{l}}_{{\ensuremath{Y}}'})=\alpha_{5}-\alpha_{6}-\alpha_{7}+\alpha_{8}$ with $$\begin{aligned}
\alpha_{5}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) E(\| {\ensuremath{Y}}-{\ensuremath{Y}}' \|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}),\\
\alpha_{6}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) E(\| {\ensuremath{Y}}''-{\ensuremath{Y}}' \|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}),\\
\alpha_{7}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) E(\| {\ensuremath{Y}}-{\ensuremath{Y}}' \|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}),\\
\alpha_{8}&=E( \| {\ensuremath{X}}-{\ensuremath{X}}'' \| {\boldsymbol{I}}_{{\ensuremath{X}},{\ensuremath{X}}'}^{\rho_{k}}) E(\| {\ensuremath{Y}}''-{\ensuremath{Y}}' \|  {\boldsymbol{I}}_{{\ensuremath{Y}}',{\ensuremath{Y}}}^{\rho_{l}}).\end{aligned}$$ All the $\alpha$’s are bounded due to the finite first moment assumption on $({\ensuremath{X}},{\ensuremath{Y}})$. Note that for distance covariance, one can go through the same proof with only three terms – $\alpha_{1}, \alpha_{2}, \alpha_{5}$ – while the local version involves eight terms, due to the additional random variables for local scales.

For the sample local covariance, the expectation of the first term can be expanded as $$\begin{aligned}
&\ \frac{1}{n(n-1)}\sum_{i \neq j}^{n}E(A_{ij}B_{ji}{\boldsymbol{I}}(R^{A}_{ij} \leq k){\boldsymbol{I}}(R^{B}_{ji} \leq l))\\
=&\ E( (\frac{n-2}{n-1}\tilde{A}_{ij}-\frac{1}{n-1}\sum_{s \neq i,j} \tilde{A}_{sj}) \\
&\ \cdot (\frac{n-2}{n-1}\tilde{B}_{ji}-\frac{1}{n-1}\sum_{s \neq i,j} \tilde{B}_{si}) {\boldsymbol{I}}(R^{A}_{ij} \leq k){\boldsymbol{I}}(R^{B}_{ji} \leq l))\\
=&\ \frac{(n-2)^2}{(n-1)^2} (\alpha_{1}-\alpha_{2}-\alpha_{3})+\frac{(n-2)(n-3)}{(n-1)^2}\alpha_{4}+\mathcal{O}(\frac{1}{n}) \\
=&\ \alpha_{1}-\alpha_{2}-\alpha_{3}+\alpha_{4}+\mathcal{O}(\frac{1}{n}).\end{aligned}$$ The expectation of the second term can be similarly expanded as $$\begin{aligned}
& E(\frac{1}{n(n-1)}\sum_{i \neq j}^{n}A^{k}_{ij} \frac{1}{n(n-1)}\sum_{i \neq j}^{n}B^{l}_{ji}) \\
=&\ \frac{1}{n^2(n-1)^2} \sum_{u\neq v}^{n} E(A_{uv}{\boldsymbol{I}}(R^{A}_{uv} \leq k) \sum_{i \neq j}^{n}B_{ji}{\boldsymbol{I}}(R^{B}_{ji} \leq l))\\
=&\ \frac{1}{n(n-1)} E( (\frac{n-2}{n-1}\tilde{A}_{uv}-\frac{1}{n-1}\sum_{s \neq u,v} \tilde{A}_{sv}){\boldsymbol{I}}(R^{A}_{uv} \leq k)  \\
&\cdot \sum_{i \neq j}^{n} (\frac{n-2}{n-1}\tilde{B}_{ji}-\frac{1}{n-1}\sum_{s \neq i,j} \tilde{B}_{si}) {\boldsymbol{I}}(R^{B}_{ji} \leq l)\\
=&\ \alpha_{5}-\alpha_{6}-\alpha_{7}+\alpha_{8}+\mathcal{O}(\frac{1}{n}).\end{aligned}$$ Combining the results yields that $E(dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}))=dCov^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})+\mathcal{O}(\frac{1}{n})$.

**(ii)**: The variance of sample local covariance is computed as $$\begin{aligned}
&Var( {\hat{E}}( A^{k}- {\hat{E}}(A^{k})) (B^{l'}- {\hat{E}}(B^{l'})))\\
=&\ \frac{1}{n^2(n-1)^2} Var(\sum_{i \neq j}^{n} (A^{k}_{ij}- {\hat{E}}(A^{k})) (B^{l}_{ji}- {\hat{E}}(B^{l})))\\
=&\ \frac{n^4}{n^2(n-1)^2} \mathcal{O}(\frac{1}{n}) + \frac{n^3}{n^2(n-1)^2} \mathcal{O}(1).\end{aligned}$$ The last equality follows because: there are $n^4$ covariance terms in the numerator of $\mathcal{O}(\frac{1}{n})$, because $Cov( (A^{k}_{ij}- {\hat{E}}(A^{k})) (B^{l}_{ji}- {\hat{E}}(B^{l})), (A^{k}_{uv}- {\hat{E}}(A^{k})) (B^{l}_{vu}- {\hat{E}}(B^{l})))$ are only related via the column centering when $(i,j)$ does not equal $(u,v)$; and there remains $n^3$ covariance terms of at most $\mathcal{O}(1)$. Note that the finite second moment assumption of $({\ensuremath{X}},{\ensuremath{Y}})$ is required for the big $\mathcal{O}$ notation to have a bounding constant. Therefore, the variance of sample local covariance is of $\mathcal{O}(\frac{1}{n})$.

**(iii)**: $dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})$ converges to the population local covariance by applying the strong law of large numbers on U-statistics [@KoroljukBook]. Namely, the first term of sample local covariance satisfies $$\begin{aligned}
&\frac{1}{n(n-1)}\sum_{i \neq j}^{n}A_{ij}B_{ji}{\boldsymbol{I}}(R^{A}_{ij} \leq k){\boldsymbol{I}}(R^{B}_{ji} \leq l) \\
=&\ \frac{1}{n}\sum_{i=1}^{n}(\frac{1}{n-1}\sum_{j \neq i}^{n} (\frac{n-2}{n-1}\tilde{A}_{ij}-\frac{1}{n-1}\sum_{s \neq i,j} \tilde{A}_{sj}) \\
& \cdot (\frac{n-2}{n-1}\tilde{B}_{ji}-\frac{1}{n-1}\sum_{s \neq i,j} \tilde{B}_{si}) {\boldsymbol{I}}(R^{A}_{ij} \leq k){\boldsymbol{I}}(R^{B}_{ji} \leq l) ) \\
\rightarrow &\ \frac{1}{n}\sum_{i=1}^{n} (\alpha_{1|(x_i,y_i)}-\alpha_{2|(x_i,y_i)}-\alpha_{3|(x_i,y_i)}+\alpha_{4|(x_i,y_i)})\\
\rightarrow &\ \alpha_{1}-\alpha_{2}-\alpha_{3}+\alpha_{4},\end{aligned}$$ where the second line applies law of large numbers at each $i$ by conditioning on $({\ensuremath{X}},{\ensuremath{Y}})=(x_i,y_i)$ for each $\alpha$’s, and the last line follows by applying law of large numbers to the independently distributed conditioned $\alpha$’s. Similarly, the second term of sample local covariance can be shown to converge to the second term in population local covariance. The convergence is also uniform: each local covariance are dependent with each other, and actually repeats the summands with each other. Thus there exists a scale $(k,l)$ such that $dCor^{k,l}$ has the largest deviation from the mean than all other local covariances, and one can find a suitable $\epsilon$ to bound the maximum deviation for all $dCor^{k,l}$.

Alternatively, convergence in probability can be directly established from (i) and (ii) by applying the Chebyshev’s inequality; the almost sure convergence can also be proved via the integral definition using almost the same steps as in Theorems 1 and 2 from [@SzekelyRizzoBakirov2007], i.e., first define the empirical characteristic function via the $w$ integral for the sample local covariance, and show it converges to the population local covariance in Equation \[eq:dcov1\] by the law of large numbers on U-statistics.

Corollary \[thm5\] {#corollarythm5 .unnumbered}
------------------

It follows directly from Theorem \[thm2\], Theorem \[thm4\], and the convergence of sample distance correlation to the population [@SzekelyRizzoBakirov2007].

Corollary \[cor1\] {#corollarycor1 .unnumbered}
------------------

The population [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} and its equivalence to expectation of Euclidean distances can be established via the same steps as in Theorem \[thm1\]. The convergence of sample [[[<span style="font-variant:small-caps;">Mantel</span>]{}]{}]{} to its population version can be derived based on either the same procedure in Theorem \[thm4\], or Theorems 1 and 2 from [@SzekelyRizzoBakirov2007] with minimal notational changes.

Theorem \[thm6\] {#theoremthm6 .unnumbered}
----------------

**(a)**: Regardless of the threshold choice, the algorithm enforces Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} to be always no less than $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n})$, and no more than $\max\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\}$.

**(b)**: By Corollary \[thm5\], $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow dCorr({\ensuremath{X}},{\ensuremath{Y}})$, then the uniform convergence by Theorem \[thm4\] ensures that $\max\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\} \rightarrow {c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$. When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, $dCorr({\ensuremath{X}},{\ensuremath{Y}})$ and ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$ are both $0$, to which Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} must converge; when dependent, $dCorr^{n,n}(\mathcal{X}_{n},\mathcal{Y}_{n})$ converges to a positive constant, so Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} must converge to a constant that is either the same or larger.

Theorem \[thm7\] {#theoremthm7 .unnumbered}
----------------

**(a)**: Given ${c}^{*}({\ensuremath{X}},{\ensuremath{Y}})>dCorr({\ensuremath{X}},{\ensuremath{Y}})$, by the continuity of local correlations with respect to $(\rho_{k},\rho_{l})$, there always exists a non-empty connected area $\mathcal{R} \in \mathcal{S}_{\epsilon}$ such that $dCorr^{\rho_k,\rho_l}({\ensuremath{X}},{\ensuremath{Y}})>dCorr({\ensuremath{X}},{\ensuremath{Y}})$ for all $(\rho_{k},\rho_{l}) \in \mathcal{R}$. Among all possible areas we take the one with largest area.

As $n$ increases to infinity, the set $\{(\frac{k-1}{n-1},\frac{l-1}{n-1}) \ | \ (k, l) \in [n]^2\}$ is a dense subset of $[0,1] \times [0,1]$, and $\{dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})\}$ is also a dense subset of $\{dCorr^{\rho_{k},\rho_{l}}({\ensuremath{X}},{\ensuremath{Y}})\}$. Thus for $n$ sufficiently large, the area $\mathcal{R}$ can always be approximated via the largest connected component $R$ by the Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} algorithm. As all sample local correlations within the region $R$ are larger than the sample distance correlation, so is the smoothed maximum. Note that if the threshold $\tau_n$ does not converge to $0$, e.g., if $\tau_n$ is a positive constant like $0.05$, Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} will fail to identify a region $R$ when $0.05 > {c}^{*}({\ensuremath{X}},{\ensuremath{Y}})$.

**(b)**: Following (a), if optimal scale of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is in the largest area $\mathcal{R}$, the sample maximum within $R$ converges to the true maximum within $\mathcal{R}$, i.e., Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} converges to the population [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{}.

Corollary \[cor2\] {#corollarycor2 .unnumbered}
------------------

For $v=\frac{n(n-3)}{2}$, $z \sim Beta(\frac{v-1}{2})$, the convergence of $\tau_n= 2F^{-1}_{z}(1-\frac{0.02}{n}) -1$ can be shown as follows: by computing the variance of the Beta distribution and using Chebyschev’s inequality, it follows that $$\begin{aligned}
& \frac{0.04}{n} = Prob(|z-0.5| \geq \tau_n /2) \leq \mathcal{O}(\frac{1}{n^2 \tau_{n}^{2}})\\
\Rightarrow &\ \tau_{n}=\mathcal{O}(\frac{1}{\sqrt{n}}) \rightarrow 0.\end{aligned}$$ The equation also implies that the percentile choice can be either fixed or anything no larger than $1-\frac{c}{n^2}$ for some constant $c$, beyond which the convergence of $\tau_n$ to $0$ will be broken.

Theorem \[thm8\] {#theoremthm8 .unnumbered}
----------------

To prove consistency under the permutation test, it suffices to show that at any type $1$ error level $\alpha$, the p-value of [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is asymptotically less than $\alpha$. The p-value can be expressed by: $$\begin{aligned}
& Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})) \\
= &\ \sum_{j=0}^{n} Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}) | \pi \mbox{ is a partial derangement of size $j$})\\
& \times Prob(\mbox{partial derangement of size $j$})  \end{aligned}$$ by conditioning on the permutation being a partial derangement of size $j$, e.g., $j=0$ means $\pi$ is a derangement, while $j=n$ means $\pi$ does not permute any position.

As $n \rightarrow \infty$, we always have $$\begin{aligned}
&Prob(\mbox{partial derangement of size $j$}) \rightarrow e^{-1} / j!, \\
&{c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}) \rightarrow \epsilon >0 \ \mbox{ under dependence.}\end{aligned}$$ Thus it suffices to show that for any $\epsilon >0$, $$\begin{aligned}
\label{eq:permconsistency}
\lim_{n\rightarrow \infty} e^{-1}\sum_{j=0}^{n} Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > \epsilon | \mbox{ partial derangement of size $j$}) / j! \rightarrow 0.\end{aligned}$$ Then we decompose the above summations into two different cases. The first case is when $j$ is a fixed size, $\mathcal{X}_{n}$ and $\mathcal{Y}_{n}^{\pi}$ are asymptotically independent (due to the *iid* assumption), thus ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi})$ converges to $0$. The other case is the remaining partial derangements $\pi$ of size $\mathcal{O}(n)$, but these partial derangements occur with probability converging to $0$, i.e., for any $\alpha > 0$, there exists $N_{1}$ such that $$\begin{aligned}
e^{-1} \sum_{j=N_{1}+1}^{+\infty} 1/j! < \alpha / 2,\end{aligned}$$ as $\sum\limits_{j=0}^{n} 1/j!$ is bounded above and converges to $e$. Then back to the first case, there further exists $N_{2}>N_{1}$ such that for any $j\leq N_{1}$ and all $n > N_{2}$ $$\begin{aligned}
Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > \epsilon)| \mbox{ partial derangement of size $j$}) < \alpha / 2.\end{aligned}$$ It follows that for all $n > N_{2}$, $$\begin{aligned}
& e^{-1} \sum_{j=0}^{n} Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > \epsilon | \mbox{ partial derangement of size $j$}) / j!\\
< &\ e^{-1} \sum_{j=0}^{N_{1}} \alpha / 2 j! + e^{-1} \sum_{j=N_{1}+1}^{n} 1 / j!\\
< &\ \alpha.\end{aligned}$$ Thus the convergence in Equation \[eq:permconsistency\] holds.

Therefore, at any type $1$ error level $\alpha>0$, the p-value of Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} under the permutation test will eventually be less than $\alpha$ as $n$ increases, such that Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} always successfully detects any dependency. Thus Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is consistent against all dependencies with finite second moments. When ${\ensuremath{X}}$ and ${\ensuremath{Y}}$ are independent, each column of $\mathcal{X}_{n}$ and the corresponding column of $\mathcal{Y}_{n}$ are independent for any permutation. Therefore, ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi})$ distributes the same as ${c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n})$ for any random permutation $\pi$, and $Prob({c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}^{\pi}) > {c}^{*}(\mathcal{X}_{n},\mathcal{Y}_{n}))$ is uniformly distributed in $[0,1]$. Thus Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is valid.

Lemma \[lem1\] {#lemmalem1 .unnumbered}
--------------

$$\begin{aligned}
dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) &= {\hat{E}}(A^{k} \circ B^{l'})- {\hat{E}}(A^{k} \circ J){\hat{E}}(B^{l} \circ J) \\
&= tr (A^{k}B^{l})- tr (A^{k}J)tr(B^{l}J)\\
&=tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)]\\
& = \sum_{i=1}^{n} \lambda_{i} [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)],\end{aligned}$$

where the first line is the definition, the second line follows by noting that ${\hat{E}}(A \circ B^{'})=tr(AB)$ and ${\hat{E}}(A)={\hat{E}}(A \circ J)=tr(AJ)$ for any two matrices $A$ and $B$, and the last two lines follow from basic properties of matrix trace.

Theorem \[thm:dvar\] {#theoremthmdvar .unnumbered}
--------------------

For all these properties, it suffices to prove them on the sample local variance $dVar^{k}(\mathcal{X}_{n})$ first. Then the population version follows by the convergence property in Theorem \[thm4\].

**(a)**: Based on Lemma \[lem1\] it holds that $$\begin{aligned}
dVar^{k}(\mathcal{X}_{n}) 
& = \sum_{i=1}^{n} \lambda^{2}_{i}[A^{k}- tr (A^{k}J)J] \geq 0.\end{aligned}$$

**(b)**: Following part (a), we have $$\begin{aligned}
&\ dVar^{k}(\mathcal{X}_{n})=0\\
\Leftrightarrow & \ \lambda_{i}[A^{k}- tr (A^{k}J)J] =0, \ \forall i \\
\Leftrightarrow & \ A^{k}- tr (A^{k}J)J = 0_{n \times n} \\
\Leftrightarrow & \ A^{k}_{ij} = tr (A^{k}J), \ \forall i,j =1,\ldots, n\\
\Leftrightarrow & \ A^{k}_{ij}=tr (A^{k}J)=0, \ \forall i,j =1,\ldots, n,\end{aligned}$$ where the last line follows by observing that $A^{k}_{ii}=0$ by Equation \[localCoef2\]. Therefore, distance variance equals $0$ if and only if $A^{k}$ is the zero matrix.

A trivial case is $k=0$, which corresponds to $\rho_k =0$ asymptotically. Otherwise $A^{k}$ is a zero matrix if and only if for all $(i,j)$ satisfying ${\boldsymbol{I}}(R^{A}_{ij} \leq k)=1$, $$\begin{aligned}
\tilde{A}_{ij}=\frac{1}{n-1}\sum_{s=1}^{n} \tilde{A}_{sj}. \end{aligned}$$ Namely, for each point $x_j$, its $k$ smallest distance entries all equal the mean distances with respect to $x_j$, which can only happen when $\tilde{A}_{ij}$ is a constant for all $i \neq j$ at a fixed $j$. Due to the symmetry of the distance matrix, all the off-diagonal entries of $\tilde{A}$ are the same, i.e., $\tilde{A}=u (J - I)$ for some constant $u \geq 0$.

When $u=0$, all observations are the same, so ${\ensuremath{X}}$ is a constant. Otherwise all observations are equally distanced from each other by a distance of $u>0$, which occurs with probability $0$ under the *iid* assumption. This is because when ${\ensuremath{X}}^{'}$ and ${\ensuremath{X}}^{''}$ are independent, one cannot have $\|{\ensuremath{X}}^{''}-{\ensuremath{X}}\|=\|{\ensuremath{X}}^{'}-{\ensuremath{X}}\|$ almost surely unless they are degenerate.

From another point of view, for given sample data that happens to be equally distanced, e.g., $n$ points in $n-1$ dimensions, sample variances can still be $0$. But this scenario occurs with probability $0$ when each observation is assumed *iid*.

**(c)**: This follows trivially from the definition, because upon the transformation the distance matrix is unchanged up-to a factor of $u$.

Theorem \[thm:dcor\] {#theoremthmdcor .unnumbered}
--------------------

Similar as in Theorem \[thm:dvar\], it suffices to prove (a) and (b) for the sample local correlation, then they automatically hold for the population version by convergence.

**(a)**: The symmetric part is trivial: for any $(\rho_k,\rho_l) \in [0,1] \times [0,1]$, by Lemma \[lem1\] $$\begin{aligned}
dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) 
&=tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)] \\
&=tr [ (B^{l}-tr(B^{l}J)J)(A^{k}- tr (A^{k}J)J)] \\
&=dCov^{l,k}(\mathcal{Y}_{n},\mathcal{X}_{n}).\end{aligned}$$ Then by the Cauchy-Schwarz inequality on the trace, $$\begin{aligned}
|dCov^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) |
=& \ |tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)]|\\
= & \ \sqrt{tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)]}\\
& \ \times  \sqrt{tr [ (A^{k}- tr (A^{k}J)J) (B^{l}-tr(B^{l}J)J)]} \\
\leq & \ \sqrt{ dVar^{k}(\mathcal{X}_{n}) dVar^{l}(\mathcal{Y}_{n})}.\end{aligned}$$ Thus $dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n}) =dCorr^{l,k}(\mathcal{Y}_{n},\mathcal{X}_{n}) \in [-1,1]$.

**(b)**: The if direction is clear: under isometry, $\tilde{A}=|u| \tilde{B}$, both share the same k-nearest-neighbor graph, so that $A^{k} = |u| \cdot B^{k}$. Thus $dCov^{k,k}(\mathcal{X}_{n},\mathcal{Y}_{n}) = \frac{1}{u^2} dVar^{k}(\mathcal{X}_{n})= u^2 \cdot dVar^{k}(\mathcal{Y}_{n})$, and $dCorr^{k,k}(\mathcal{X}_{n},\mathcal{Y}_{n})=1$. For the only if direction: by part (a), the local correlation can be $\pm 1$ if and only if $(A^{k}- tr (A^{k}J)J)$ is a scalar multiple of $(B^{l}-tr(B^{l}J)J)$, say some constant $u$.

First we argue that the non-zero entries in $A^{k}$ must match the non-zero entries in $B^{l}$. Namely, the k-nearest neighbor graph is the same between $\tilde{A}$ and $\tilde{B}$. As $A^{k}_{ii}=B^{l}_{ii}=0$, $-tr (A^{k}J)$ must be a scalar multiple of $-tr (B^{l}J)$. Then if there exists $i \neq j$ such that $A^{k}_{ij}=0$ while $B^{l}_{ij} \neq 0$, $-tr (A^{k}J)$ must be the same scalar multiple of $B^{l}_{ij}-tr (B^{l}J)$, which is not possible unless $B^{l}_{ij}=0$. Thus $k=l$ and ${\boldsymbol{I}}(R^{A}_{ij} \leq k)= {\boldsymbol{I}}(R^{B}_{ij} \leq k)$ for all $(i,j)$.

Next we show the scalar multiple must be positive, i.e., the local correlation cannot be $-1$. Assuming it can be $-1$, then $$\begin{aligned}
& \ A^{k}- tr (A^{k}J)J=-|u|(B^{k}-tr(B^{k}J)J) \\
\Leftrightarrow & \ A^{k}+|u|B^{k}=(tr (A^{k}J)+|u|tr(B^{k}J))J \\
\Leftrightarrow & \ A^{k}+|u|B^{k}=0_{n \times n} \\
\Leftrightarrow & \ A+|u|B=0_{n \times n}, \end{aligned}$$ where the second to last line follows because the diagonal entries of $A^{k}+|u|B^{k}$ are $0$ by definition, and the last line follows by observing that $tr (A^{k}J)$ and $tr (B^{k}J)$ are both negative unless $k=n$ (e.g., $A$ is always centered to have zero matrix mean, while $A^{k}$ keeps the $k$ smallest entries per column so its matrix mean is negative til $k=n$). However, if the last line is true, then the original distance correlation shall be $-1$, which cannot happen under the *iid* assumption as shown in [@SzekelyRizzoBakirov2007]. Note that the derivation also shows that the local correlations can be $-1$ for general dissimilarity matrices without the *iid* assumption, i.e., when $\tilde{A}+|u| \tilde{B}=v(J-I)$ for some constant $v$.

Therefore, the scalar multiple must be positive, and $A^{k}- tr (A^{k}J)J = |u| (B^{k}-tr (B^{k}J)J)$. As the diagonals satisfy $A^{k}_{ii}=B^{k}_{ii}=0$, it holds that $tr (A^{k}J)= |u| tr (B^{k}J)$ and $A^{k}=|u| B^{k}$. Thus for each $(i, j)$ satisfying ${\boldsymbol{I}}(R^{A}_{ij} \leq k)=1$: $$\begin{aligned}
& \ \tilde{A}_{ij}-\frac{1}{n-1}\sum_{s=1}^{n} \tilde{A}_{sj}=|u|(\tilde{B}_{ij}-\frac{1}{n-1}\sum_{s=1}^{n} \tilde{B}_{sj}) \\
\Leftrightarrow & \ \tilde{A}_{ij}-|u|\tilde{B}_{ij} = \frac{1}{n-1}\sum_{s=1}^{n} \tilde{A}_{sj} -\frac{|u|}{n-1}\sum_{s=1}^{n} \tilde{B}_{sj}\\
\Leftrightarrow & \ \tilde{A}_{ij}-|u|\tilde{B}_{ij} = v.\end{aligned}$$ We argue that if $\tilde{A}_{ij}=|u|\tilde{B}_{ij}+v$ for each $(i, j)$ satisfying ${\boldsymbol{I}}(R^{A}_{ij} \leq k)=1$, it also holds for all $(i,j)$. Suppose there exists $(s,j)$ with ${\boldsymbol{I}}(R^{A}_{sj} \leq k)=0$ and $\tilde{A}_{sj} = |u|\tilde{B}_{sj}+v+w$ for some $w \neq 0$. Without loss of generality, there must exist one more index $t$ such that ${\boldsymbol{I}}(R^{A}_{tj} \leq k)=0$ and $\tilde{A}_{tj} = |u|\tilde{B}_{tj}+v-w$ to maintain the mean (or multiple indices in a similar manner). This requires $\|{\ensuremath{X}}^{''}-{\ensuremath{X}}\|-|u|\|{\ensuremath{Y}}^{''}-{\ensuremath{Y}}\|=\|{\ensuremath{X}}^{'}-{\ensuremath{X}}\|-|u|\|{\ensuremath{Y}}^{'}-{\ensuremath{Y}}\|+2w$, so $({\ensuremath{X}}^{''},{\ensuremath{Y}}^{''})$ and $({\ensuremath{X}}^{'},{\ensuremath{Y}}^{'})$ are related by $w$ when conditioning on $({\ensuremath{X}},{\ensuremath{Y}})$. Thus it imposes a dependency structure and violates the *iid* assumption.

Therefore $\tilde{A}-|u|\tilde{B}=v(J-I)$. When $v=0$, $\tilde{A}=|u|\tilde{B}$ is equivalent to that $({\ensuremath{X}}, u {\ensuremath{Y}}$) are related by an isometry. When $v \neq 0$, it requires each distance entries to be added by the same constant, which occurs with probability $0$ under the *iid* assumption. Namely, if $\|{\ensuremath{X}}^{'}-{\ensuremath{X}}\|-|u| \|{\ensuremath{Y}}^{'}-{\ensuremath{Y}}\|=\|{\ensuremath{X}}^{''}-{\ensuremath{X}}\|-|u| \|{\ensuremath{Y}}^{''}-{\ensuremath{Y}}\|=v \neq 0$ almost surely, then $({\ensuremath{X}}^{''},{\ensuremath{Y}}^{''})$ and $({\ensuremath{X}}^{'},{\ensuremath{Y}}^{'})$ are related by $v$ when conditioning on $({\ensuremath{X}},{\ensuremath{Y}})$, in which case these two pairs become dependent and the *iid* assumption is violated.

**(c)**: As each local correlation is symmetric and bounded for either population or sample case, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} is symmetric and within $[-1,1]$ by part (a).

**(d)**: If ${\ensuremath{X}}$ and $u {\ensuremath{Y}}$ are related by an isometry, the distance correlation (or the local correlation at the largest scale) equals $1$. For population, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} takes the maximum local correlation; for sample, [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} cannot be smaller than the local correlation at the largest scale. In both cases population and Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equal $1$.

When population or Sample [[[<span style="font-variant:small-caps;">MGC</span>]{}]{}]{} equal $1$, there exists at least one local correlation that equals $1$, i.e., $dCorr^{k,l}(\mathcal{X}_{n},\mathcal{Y}_{n})=1$. From the inequality in part (a), $k$ must equal $l$ for the equality to hold. Otherwise the number of non-zero entries does not match between $A^{k}$ and $B^{l}$, and $A^{k}$ cannot be a scalar multiple of $B^{l}$. Thus there exists $k$ such that $dCorr^{k,k}(\mathcal{X}_{n},\mathcal{Y}_{n})=1$, and the conclusion follows from part (b).

Simulation Dependence Functions {#appen:function}
===============================

This section presents the $20$ simulations used in the experiment section, which is mostly based on a combination of simulations from previous works [@SzekelyRizzoBakirov2007; @SimonTibshirani2012; @GorfineHellerHeller2012]. We only made changes to add noise and a weight vector for higher dimensions, thereby making them more difficult and easier to compare all methods throughout different dimensions and sample sizes. For the random variable ${\ensuremath{X}}\in {\mathbb{R}}^{p}$, we denote ${\ensuremath{X}}_{[d]}, d=1,\ldots,p$ as the $d^{th}$ dimension of [$X$]{}. For the purpose of high-dimensional simulations, $w \in {\mathbb{R}}^{p}$ is a decaying vector with $w_{[d]}=1/d$ for each $d$, such that $w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}$ is a weighted summation of all dimensions of [$X$]{}. Furthermore, ${\mathcal{U}}(a,b)$ denotes the uniform distribution on the interval $(a,b)$, ${\mathcal{B}}(p)$ denotes the Bernoulli distribution with probability $p$, ${\mathcal{N}}(\mu,{\Sigma})$ denotes the normal distribution with mean ${\mu}$ and covariance ${\Sigma}$, $U$ and $V$ represent some auxiliary random variables, $\kappa$ is a scalar constant to control the noise level (which equals $1$ for one-dimensional simulations and $0$ otherwise), and $\epsilon$ is sampled from an independent standard normal distribution unless mentioned otherwise.

Linear $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&=w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}+\kappa\epsilon.\end{aligned}$$

Exponential $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(0,3)^{p}, \\
{\ensuremath{Y}}&=exp(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}})+10\kappa\epsilon.\end{aligned}$$

Cubic $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p}, \\
{\ensuremath{Y}}&=128(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}-\tfrac{1}{3})^3+48(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}-\tfrac{1}{3})^2-12(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}-\tfrac{1}{3})+80\kappa\epsilon.\end{aligned}$$

Joint normal $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: Let $\rho=1/2p$, $I_{p}$ be the identity matrix of size $p \times p$, $J_{p}$ be the matrix of ones of size $p \times p$, and $\Sigma = \begin{bmatrix} I_{p}&\rho J_{p}\\ \rho J_{p}& (1+0.5\kappa) I_{p} \end{bmatrix}$. Then $$\begin{aligned}
({\ensuremath{X}}, {\ensuremath{Y}}) &\sim {\mathcal{N}}(0, \Sigma).\end{aligned}$$

Step Function $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&={\boldsymbol{I}}(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}>0)+\epsilon,\end{aligned}$$ where ${\boldsymbol{I}}$ is the indicator function, that is ${\boldsymbol{I}}(z)$ is unity whenever $z$ true, and zero otherwise.

Quadratic $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&=(w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}})^2+0.5\kappa\epsilon.\end{aligned}$$

W Shape $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $U \sim {\mathcal{U}}(-1,1)^{p}$, $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&=4\left[ \left( (w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}})^2 - \tfrac{1}{2} \right)^2 + w{^{\ensuremath{\mathsf{T}}}}U/500 \right]+0.5\kappa\epsilon.\end{aligned}$$

Spiral $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $U \sim {\mathcal{U}}(0,5)$, $\epsilon \sim {\mathcal{N}}(0, 1)$, $$\begin{aligned}
{\ensuremath{X}}_{[d]}&=U \sin(\pi U)  \cos^{d}(\pi U) \mbox{ for $d=1,\ldots,p-1$},\\
{\ensuremath{X}}_{[p]}&=U \cos^{p}(\pi U),\\
{\ensuremath{Y}}&= U \sin(\pi U) +0.4 p\epsilon.\end{aligned}$$

Uncorrelated Bernoulli $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $U \sim {\mathcal{B}}(0.5)$, $\epsilon_{1} \sim {\mathcal{N}}(0, I_{p})$, $\epsilon_{2} \sim {\mathcal{N}}(0, 1)$, $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{B}}(0.5)^{p}+0.5\epsilon_{1},\\
{\ensuremath{Y}}&=(2U-1)w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}+0.5\epsilon_{2}.\end{aligned}$$

Logarithmic $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: $\epsilon \sim {\mathcal{N}}(0, I_{p})$ $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{N}}(0, I_{p}),\\
{\ensuremath{Y}}_{[d]}&=2\log_{2}(|{\ensuremath{X}}_{[d]}|)+3\kappa\epsilon_{[d]} \mbox{ for $d=1,\ldots,p$.}\end{aligned}$$

Fourth Root $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&=|w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}}|^\frac{1}{4}+\frac{\kappa}{4}\epsilon.\end{aligned}$$

Sine Period $4\pi$ $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: $U \sim {\mathcal{U}}(-1,1)$, $V \sim {\mathcal{N}}(0,1)^{p}$, $\theta=4\pi$, $$\begin{aligned}
{\ensuremath{X}}_{[d]}&=U+0.02 p V_{[d]} \mbox{ for $d=1,\ldots,p$}, \\
{\ensuremath{Y}}&=\sin ( \theta {\ensuremath{X}})+\kappa\epsilon.\end{aligned}$$

Sine Period $16\pi$ $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: Same as above except $\theta=16\pi$ and the noise on ${\ensuremath{Y}}$ is changed to $0.5\kappa\epsilon$.

Square $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: Let $U \sim {\mathcal{U}}(-1,1)$, $V \sim {\mathcal{U}}(-1,1)$, $\epsilon \sim {\mathcal{N}}(0,1)^{p}$, $\theta=-\frac{\pi}{8}$. Then $$\begin{aligned}
{\ensuremath{X}}_{[d]}&=U \cos\theta + V \sin\theta + 0.05 p\epsilon_{[d]},\\
{\ensuremath{Y}}_{[d]}&=-U \sin\theta + V \cos\theta,\end{aligned}$$ for $d=1,\ldots,p$.

Two Parabolas $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $\epsilon \sim {\mathcal{U}}(0,1)$, $U \sim {\mathcal{B}}(0.5)$, $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{U}}(-1,1)^{p},\\
{\ensuremath{Y}}&=\left( (w{^{\ensuremath{\mathsf{T}}}}{\ensuremath{X}})^2  + 2\kappa\epsilon\right) \cdot (U-\tfrac{1}{2}).\end{aligned}$$

Circle $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: $U \sim {\mathcal{U}}(-1,1)^{p}$, $\epsilon \sim {\mathcal{N}}(0, I_{p})$, $r=1$, $$\begin{aligned}
{\ensuremath{X}}_{[d]}&=r \left(\sin(\pi U_{[d+1]})  \prod_{j=1}^{d} \cos(\pi U_{[j]})+0.4 \epsilon_{[d]}\right) \mbox{ for $d=1,\ldots,p-1$},\\
{\ensuremath{X}}_{[p]}&=r \left(\prod_{j=1}^{p} \cos(\pi U_{[j]})+0.4 \epsilon_{[p]}\right),\\
{\ensuremath{Y}}&= \sin(\pi U_{[1]}).\end{aligned}$$

Ellipse $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}$: Same as above except $r=5$.

Diamond $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: Same as “Square” except $\theta=-\frac{\pi}{4}$.

Multiplicative Noise $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: $U \sim {\mathcal{N}}(0, I_{p})$, $$\begin{aligned}
{\ensuremath{X}}&\sim {\mathcal{N}}(0, I_{p}),\\
{\ensuremath{Y}}_{[d]}&=U_{[d]}{\ensuremath{X}}_{[d]} \mbox{  for $d=1,\ldots,p$.}\end{aligned}$$

Multimodal Independence $({\ensuremath{X}},{\ensuremath{Y}}) \in {\mathbb{R}}^{p} \times {\mathbb{R}}^{p}$: Let $U \sim {\mathcal{N}}(0,I_{p})$, $V \sim {\mathcal{N}}(0,I_{p})$, $U' \sim {\mathcal{B}}(0.5)^{p}$, $V' \sim {\mathcal{B}}(0.5)^{p}$. Then $$\begin{aligned}
{\ensuremath{X}}&=U/3+2U'-1,\\
{\ensuremath{Y}}&=V/3+2V'-1.\end{aligned}$$

For the increasing dimension simulations in the main paper, we always set $\kappa=0$ and $n=100$, with $p$ increasing. For types $4,10,12,13,14,18,19,20$, $q=p$ such that $q$ increases as well; otherwise $q=1$. The decaying vector $w$ is utilized for $p>1$ to make the high-dimensional relationships more difficult (otherwise, additional dimensions only add more signal). For the one-dimensional simulations, we always set $p=q=1$, $\kappa=1$ and $n=100$.

[^1]: shenc@udel.edu

[^2]: cep@jhu.edu

[^3]: jovo@jhu.edu

[^4]: <https://github.com/neurodata/mgc-matlab>

[^5]: <https://CRAN.R-project.org/package=mgc>
