A fast algorithm for computing distance correlation

1 Introduction

Detecting dependencies between two random vectors $X$ and $Y$

is a fundamental problem in statistics and machine learning. Dependence measures such as Pearson’s correlation, Spearman’s

ρ

, and Kendall’s

τ

are used in almost all quantitative areas; example areas are bioinformatics (Guo et al., 2014; Sferra et al., 2017) and time-series (Zhou, 2012). However, those classical dependence measures are usually designed to detect one specific dependence structure such as a monotonic or linear structure. It is easy to construct highly dependent

X

and

Y

whose dependence cannot be detected by classical dependence measures. To overcome these limitations, Székely et al. (2007); Székely and Rizzo (2009) proposed distance covariance as a weighted

L_{2}

distance between the joint characteristic function and the product of marginal characteristic distributions. The distance covariance is

0

if and only if two random vectors

X

and

Y

are independent. A closely related measure is the Hilbert-Schmidt independence measure (HSIC). HSIC has been extensively studied in machine learning literature (Gretton et al., 2005, 2008; Pfister et al., 2018). Sejdinovic et al. (2013) established equivalence of distance covariance with HSIC.

Despite the power of sample distance covariance to detect a dependence structure, its use for large sample sizes, is inhibited by the high computational cost required. The sample distance covariance and HSIC computation typically requires $O (n^{2})$ pairwise distance (kernel) calculations and $O (n^{2})$ memory for storing them. This is undesirable and greatly limits the application of distance correlation to large datasets. In the era of big data, it is not rare to see data that consists of millions of observations. For such data, an $O (n^{2})$ algorithm is almost impossible to run on a personal computer. To approximate the distance covariance or HSIC for large data, Nyström approach or the random Fourier feature method is often adopted. However, the use of these approximations leads to a reduction in power (Zhang et al., 2018). In this article, we describe an exact method to compute the sample distance covariance between two univariate random variables with computational cost $O (n log (n))$ and memory cost $O (n)$ . Our proposed method essentially consists of just two sorting steps, which makes it easy to implement.

A closely related $O (n log (n))$ algorithm for sample distance covariance was proposed by Huo and Székely (2016). Our algorithm differs from Huo and Székely (2016) in the following ways: First, they implicitly assume that there are no ties in the data (see Algorithm 1 and proof in Huo and Székely (2016)), whereas our proposed method is valid for any pair of real-valued univariate variables. In practice, it is common to see datasets with ties, especially for discrete variables or bootstrap sample. Second, we use a merge sort instead of an AVL tree-type implementation to compute the Frobenius inner product of the distance matrices of $x$ and $y$ . Empirical results show that our proposed method is significantly faster; for example, for one million observations our MATLAB implementation runs 10 times faster (finishing in 4 seconds) on our desktop, whereas the implementation in Huo and Székely (2016) requires 40 seconds. Because our implementation consists only of MATLAB code while the key step in the Huo and Székely (2016) routine is implemented in C, even greater speed increases are possible by rewriting the critical parts of our implementation in C.

The rest of paper is organized as follows. In Section 2

, we briefly introduce the definition of distance covariance and its sample estimate. In Section

3, we describe the proposed

O (n log (n))

algorithm for sample distance covariance. In Section 4, experiment results are presented. Finally, conclusions and remarks are made in Section 5.

2 Some Preliminaries

Denote the joint characteristic function of $X \in R^{p}$ and $Y \in R^{q}$ as $f_{X, Y} (t, s)$ , and denote the marginal characteristic functions of $X$ and $Y$ as $f_{X} (t)$ and $f_{Y} (s)$ , respectively. Denote $| \cdot |_{k}$ as the Euclidean norm in $R^{k}$ . The squared distance covariance is defined as the weighted $L_{2}$ distance between $f_{X, Y} (t, s)$ and $f_{X} (t) \cdot f_{Y} (s)$ ,

V^{2} (X, Y) = \int_{R^{p + q}} | f_{X, Y} (t, s) - f_{X} (t) f_{Y} (s) |^{2} w (t, s) d t d s,

where $w (t, s) = (c_{p} c_{q} | t |_{p}^{1 + p} | s |_{q}^{1 + q})^{- 1}$ , $c_{p}$ and $c_{q}$ are constants. It is obvious that $V^{2} (X, Y) = 0$ if and only if $X$ and $Y$ are independent. Under some mild conditions, the squared distance covariance can be defined equivalently as the expectation of Euclidean distances

	$V^{2} (X, Y) =$	$E (\| X - X^{'} \|_{p} \| Y - Y^{'} \|_{q}) - 2 E (\| X - X^{'} \|_{p} \| Y - Y^{''} \|_{q})$
		$+ E (\| X - X^{'} \|_{p}) E (\| Y - Y^{'} \|_{q}),$		(1)

where $(X, Y), (X^{'}, Y^{'}),$ and $(X^{''}, Y^{''})$

are identical independent copies from the joint distribution of

(X, Y)

The squared distance correlation is defined by

R^{2} (X, Y) = ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} \frac{V^{2} (X, Y)}{\sqrt{V^{2} (X, X) V^{2} (Y, Y)}}, & if V^{2} (X, X) V^{2} (Y, Y) > 0 0 & otherwise . \end{matrix}

(2)

Let $X = (x_{1}, \dots, x_{n})^{⊺}$ and $Y = (y_{1}, \dots, y_{n})^{⊺}$ be the sample collected. Define

	$a_{i j} = \| x_{i} - x_{j} \|_{p},$		$b_{i j} = \| y_{i} - y_{j} \|_{q}$
	$a_{i .} = n \sum j = 1 a_{i j},$		$b_{i .} = n \sum j = 1 b_{i j}$
	$a_{. .} = n \sum i = 1 a_{i .},$		$b_{. .} = n \sum i = 1 b_{i .}$
	$D = \sum 1 \leq i, j \leq n a_{i j} b_{i j}$

The squared sample distance covariance between $X$ and $Y$ is

\begin{matrix} V_{n}^{2} (X, Y) = \frac{D}{} n^{2} - \frac{2}{n^{3}} n \sum i = 1 a_{i .} b_{i .} + \frac{a_{. .} b_{. .}}{n^{4}}, \end{matrix}

(3)

which is similar in form to (2). The squared sample distance correlation is given by

R_{n}^{2} (X, Y) = ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} \frac{V_{n}^{2} (X, Y)}{\sqrt{V_{n}^{2} (X, X) V_{n}^{2} (Y, Y)}}, & if V_{n}^{2} (X, X) V_{n}^{2} (Y, Y) > 0 0 & otherwise . \end{matrix}

(4)

From (3), it is easy to see a $O (n^{2})$ brute force algorithm exists for distance covariance. However, the brute force implementation is difficult to handle large datasets. Moreover, the p-value of distance covariance or correlation is typically calculated by using permutation test, which makes it more computationally intensive.

If we can compute $D$ and all $a_{i .}, b_{i .}$ for $1 \leq i \leq n$ and $D$ in $O (n log n)$ steps, then we can also compute $V_{n}^{2} (X, Y)$ in $O (n log n)$ steps.

In this paper, we consider the case where $X$ and $Y$ are univariate random variables; that is, $p = q = 1$ . For the rest of this document we assume that $x_{1} \leq x_{2} \leq \dots \leq x_{n}$ (because after an $O (n log n)$ sort step, we can ensure that $x_{1} \leq x_{2} \leq \dots \leq x_{n}$ ).

3 Fast Algorithm for Distance Covariance

Define the function $I (x)$ as

I (x) = {\begin{matrix} 1 & i f $ x > 0 $ 0 & o t h e r w i s e . \end{matrix}

(5)

For any two real $x$ and $y$ we have

| x - y | = (x - y) (2 I (x - y) - 1) .

(6)

We use (6) extensively in the rest of paper.

3.1 Fast computation of the $a_{i .}$ and $b_{i .}$

Define $s_{i} = \sum_{j = 1}^{i} x_{i}$ for $1 \leq i \leq n$ and note that $s_{1}, \dots, s_{n}$ can be computed in $O (n)$ time.

Since $x_{1} \leq x_{2} \leq \dots \leq x_{n}$ we have

	$a_{i .}$	$= \sum j < i (x_{i} - x_{j}) + \sum j > i (x_{j} - x_{i})$
		$= (2 i - n) x_{i} + (s_{n} - 2 s_{i}) .$		(7)

So $a_{1.}, \dots, a_{n .}$ can be computed in $O (n)$ time.

We can use an $O (n log n)$ sorting algorithm to determine a permutation $π (1),$ $π (2),$ $\dots,$ $π (n)$ of $1, 2, \dots, n$ such that $y_{π (1)} \leq y_{π (2)} \leq \dots \leq y_{π (n)}$ . Therefore as in (7), $b_{π (1)}, \dots, b_{π (n)}$ can be computed in $O (n)$ time after $y_{1}, \dots, y_{n}$ is sorted.

3.2 Fast computation of D

In this subsection, we describe an $O (n log (n))$ algorithm for computing $D$ . First, we have

	$D$	$= n \sum i = 1 n \sum j = 1 \| x_{i} - x_{j} \| \| y_{i} - y_{j} \|$
		$= 2 n \sum i = 1 \sum 1 \leq j < i \| x_{i} - x_{j} \| \| y_{i} - y_{j} \| .$		(8)

In (3.2) note that $| x_{i} - x_{j} | = x_{i} - x_{j}$ if $1 \leq j \leq i$ , thus showing that

$\frac{D}{2} =$	$n \sum i = 1 \sum 1 \leq j < i (x_{i} - x_{j}) (y_{i} - y_{j}) (2 I (y_{i} - y_{j}) - 1)$
$=$	$2 n \sum i = 1 \sum 1 \leq j < i (x_{i} - x_{j}) (y_{i} - y_{j}) I (y_{i} - y_{j})$
	$- n \sum i = 1 \sum 1 \leq j < i (x_{i} - x_{j}) (y_{i} - y_{j}) .$	(9)

Let $m_{x} = \sum_{i = 1}^{n} x_{i} / n$ and $m_{y} = \sum_{i = 1}^{n} y_{i} / n$ . Now the second term in (3.2) is

	$n \sum i = 1 \sum 1 \leq j < i (x_{i} - x_{j}) (y_{i} - y_{j})$	$= \frac{1}{2} \sum 1 \leq i, j \leq n (x_{i} - x_{j}) (y_{i} - y_{j})$
		$= n n \sum i = 1 (x_{i} - m_{x}) (y_{i} - m_{y}) .$

Therefore it can be computed in $O (n)$ steps.

Define $U_{i} = {j : 1 \leq j < i, y_{j} < y_{i}}$ for $1 \leq i \leq n .$ The first term in (3.2) can be easily expressed in terms of $U_{i}$ . Note that

	$n \sum i = 1 \sum 1 \leq j < i (x_{i} - x_{j}) (y_{i} - y_{j}) I (y_{i} - y_{j})$
	$= n \sum i = 1 x_{i} y_{i} \sum j \in U_{i} 1 - n \sum i = 1 x_{i} \sum j \in U_{i} y_{j} - n \sum i = 1 y_{i} \sum j \in U_{i} x_{j} + n \sum i = 1 \sum j \in U_{i} x_{j} y_{j} .$		(10)

Thus the first term in (3.2) is expanded into a sum of four terms, each of which is of the form

\begin{matrix} n \sum i = 1 s_{i} \sum j \in U_{i} t_{j} . \end{matrix}

(11)

For any $t_{1}, \dots, t_{n}$ , define $d_{i} = \sum_{j \in U_{i}} t_{j} .$ If it can be shown that all of $d_{1}, \dots, d_{n}$ can be computed in $O (n log n)$ time, then the sample distance covariance can be computed in $O (n log n)$ time. We will show this in the next subsection. The preceding arguments lead to the following theorem.

Theorem 3.1

For any real-valued univariate variables with sample $x = (x_{1}, \dots, x_{n})^{⊺}$ and $y = (y_{1}, \dots, y_{n})^{⊺}$ , the sample distance covariance $V_{n} (x, y)$ can be computed in $O (n log n)$ time.

3.3 Fast computation of $d_{i}$ where $d_{i} = \sum_{j \in U_{i}} t_{j}$

We are given a series $y_{1}, \dots, y_{n}$ along with weights $t_{1}, \dots, t_{n}$ . Define $U_{i} = {j : 1 \leq j < i, y_{j} < y_{i}}$ . Our objective is to compute $d_{1}, d_{2}, \dots, d_{n}$ in $O (n log n)$ steps where

d_{i} = \sum \begin{matrix} j < i y_{j} < y_{i} \end{matrix} t_{j} = \sum j \in U_{i} t_{j}

for $1 \leq i \leq n .$

It is well known that the number of inversions in a permutation can be obtained by a merge sort (Ginat, 2004). We use a similar strategy to compute the $d_{i}$ while performing a merge sort on $y_{1}, \dots, y_{n}$ to sort the $y_{i}$ in decreasing order.

Merge sort works by successively merging sorted subarrays until the final array is sorted. Assume that we also keep an auxiliary array for storing the original indices of each element of the array, another auxiliary array for storing partially computed $d_{i}$ (say $d [1], \dots, d [n]$ ), and third one for storing the partial sums of the intermediate results.

A merge sort makes $⌈ {log}_{2} n ⌉$ passes over the data. At the beginning of the $k^{th}$ pass, the array consists of $⌈ \frac{n}{2^{k - 1}} ⌉$ contiguous subarrays, each of which is sorted in decreasing order. We merge two consecutive subarrays $⌊ \frac{n}{2^{k}} ⌋$ times.

Let $A = [y_{n_{1}}, y_{n_{2}}, \dots, y_{n_{r}}]$ and $B = [y_{m_{1}}, y_{m_{2}}, \dots, y_{m_{s}}]$ be two such consecutive subarrays with $y_{n_{1}} \geq y_{n_{2}} \geq \dots \geq y_{n_{r}}$ and $y_{m_{1}} \geq \dots \geq y_{m_{s}}$ . We can also assume $n_{i} < m_{j}$ for all meaningful $i, j$ (here $m_{i}$ and $n_{j}$ refer to the original indices of these elements). During the merge step, if we notice that $y_{m_{α}} > y_{n_{β}}$ for some $α, β$ , then $y_{n_{β}}, y_{n_{β + 1}}, \dots, y_{n_{r}}$ are all the terms in $A$ that are less than $y_{m_{α}}$ and we increment $d [m_{α}]$ by $t_{n_{β}} + \dots + t_{n_{r}}$ . Note that if we also store running sums as we output the results, $d [m_{α}]$ can be computed using just one difference.

The extra computation of maintaining the additional auxiliary arrays does not increase the order of computation. The detailed algorithm is presented in Algorithm 1. For a better understanding of the proposed algorithm, we also include MATLAB code in the Appendix.

Input: $Y = (y_{1}, y_{2}, \dots, y_{n})$ and $T = (t_{1}, t_{2}, \dots, t_{n})$
Output: $d = (d_{1}, \dots, d_{n})$

1: Initialize

i d x

to be a

2 \times n

matrix whose first row is

(1, 2, \dots, n)

and whose second row is

(0, 0, \dots, 0)

i = 1; r = 1; s = 2;

3: while

i < n

g a p = 2 * i; k = 0; i d x_r = i d x [r,];

c s u m T = c u s u m (T [i d x_r]); c s u m T = (0, c s u m T); / * cum sum * /

j = 1;

7: while

j < n

s t 1 = j; e 1 = min (s t 1 + i - 1, n);

s t 2 = j + i; e 2 = m i n (s t 2 + i - 1, n);

10: while

s t 1 <= e 1

and

s t 2 <= e 2

11:

k = k + 1

12:

i d x 1 = i d x_r [s t 1]; i d x 2 = i d x_r [s t 2];

13: if

Y [i d x 1] >= Y [i d x 2]

then

14:

i d x [s, k] = i d x 1;

s t 1 = s t 1 + 1;

15: else

16:

i d x [s, k] = i d x 2; s t 2 = s t 2 + 1;

17:

d [i d x 2] = d [i d x 2] + (c s u m T [e 1 + 1] - c s u m T [s t 1]);

18: end if

19: end while

20: if

s t 1 <= e 1

then

21:

k f = k + e 1 - s t 1 + 1;

22:

i d x [s, (k + 1) : k f] = i d x_r [, s t 1 : e 1];

23:

k = k f;

24: else if

s t 2 <= e 2

then

25:

k f = k + e 2 - s t 2 + 1;

26:

i d x [s, (k + 1) : k f] = i d x_r [, s t 2 : e 2];

27:

k = k f;

28: end if

29:

j = j + g a p;

30: end while

31:

i = g a p;

32:

r = 3 - r; s = 3 - s;

33: end while

34: return d

Algorithm 1 Fast computation of

d_{1}, \dots, d_{n}

	index	1	2	3	4	5	6	7	8
	y	3	5	7	3	8	4	6	7
Iter 0	t	1	5	3	2	4	6	7	5
	cusumT	1	6	9	11	15	21	28	33
	d	0	0	0	0	0	0	0	0
	index	2	1	3	4	5	6	8	7
	y	5	3	7	3	8	4	7	6
Iter 1	t	5	1	3	2	4	6	5	7
	cusumT	5	6	9	11	15	21	26	33
	d	1	0	0	0	0	0	7	0
	index	3	2	1	4	5	8	7	6
	y	7	5	3	3	8	7	6	4
Iter 2	t	3	5	1	2	4	5	7	6
	cusumT	3	8	9	11	15	20	27	33
	d	6	1	0	0	0	13	6	0
	index	5	3	8	7	2	6	1	4
Iter 3	y	8	7	7	6	5	4	3	3
	d	11	6	21	14	1	3	0	0

Table 1: Illustration of fast algorithm for series

Y = (3, 5, 7, 3, 8, 4, 6, 7)^{⊺}

and weight

T = (1, 5, 3, 2, 4, 6, 7, 5)^{⊺}

. The index row denotes the original order of series, and the cusumT row denotes the cumulative sum of

t_{i}

. The proposed algorithm finishes in three iterations. Final results :

(d_{1}, \dots, d_{8}) = (0, 1, 6, 0, 11, 3, 14, 21)

We illustrate the proposed algorithm by using a simple series,
$Y = (3, 5, 7, 3, 8, 4, 6, 7)^{⊺}$ and weight $T = (1, 5, 3, 2, 4, 6, 7, 5)^{⊺}$ . The iteration history is presented in Table 1. The index row in Table 1 denotes the original order of the $Y$ series, and the cusumT row denotes the cumulative sum of $t_{i}$ . The algorithm finishes in ${log}_{2} (8) = 3$ iterations and outputs the results $(d_{1}, \dots, d_{8}) = (0, 1, 6, 0, 11, 3, 14, 21)$ . Consider the computation of $d_{7}$ and $d_{8}$ for example. At iteration 1, we merge $y_{7}, y_{8}$ , because $y_{7} < y_{8}$ we exchange $y_{7}$ and $y_{8}$ and set $d_{7} = 0$ and $d_{8} = 7$ . At iteration 2, we merge $y_{5}, y_{6}$ and $y_{8}, y_{7}$ ; and because $y_{8} > y_{6}$ and $y_{8} < y_{5}$ we increment $d_{8} = 7 + t_{6} = 13$ and because $y_{7} > y_{6}$ we increment $d_{7} = 0 + t_{6} = 6$ . At the final merge, we note $y_{8} > y_{2}$ so we increment $d_{8} = 13 + t_{2} + t_{1} + t_{4} = 13 + (11 - 3) = 21$ , where $t_{2} + t_{1} + t_{4}$ is calculated from the difference of cumulative sums of $t_{i}$ . Similarly for $d_{7}$ we have $d_{7} = 6 + t_{2} + t_{1} + t_{4} = 6 + (11 - 3) = 14.$

The proposed algorithm has the same order of computational cost as a merge sort. Denote the computational cost as $T (n)$ . We know $T (n) = 2 T (n / 2) + O (n)$ . It is trivial to show that the complexity is $O (n log (n))$ by using the master theorem (Cormen et al., 2001). As a byproduct of calculating all $d_{i}$ , $y_{1}, y_{2}, \dots, y_{n}$ are sorted. Therefore, calculating all $b_{i .}$ as in Section 3.1 takes only extra $O (n)$ time.

As discussed previously, the proposed algorithm essentially consists of two sorting steps. First we sort X and calculate $a_{i \cdot}$ for $i = 1, \dots, n$ . Then we sort Y and calculate $D$ and all $b_{i \cdot}$ for $i = 1, \dots, n$ . The correctness of the described algorithm can be easily concluded from the previous discussion. To verify the correctness of our implementation we matched our numbers with a simple brute force $O (n^{2})$ implementation and we confirmed that the numbers match exactly.

We also note another factor that impacts performance. A non recursive merge sort algorithm makes $⌈ {log}_{2} (n) ⌉$ passes over the data, and in the $k^{th}$ pass makes $\approx \frac{n}{2^{k}}$ merges of two subarrays of size $2^{k - 1}$ . If $n$ is large then for small $k$ we end up merging a large number of small subarrays, which has a large overhead. We have observed that the speed increases by a factor of $\approx 1.3$ if we replace these initial merges by a single insert sort step — that is, we divide the input array into $\approx \frac{n}{16}$ groups of length $16$ and sort each one of them using insert sort, store the intermediate $d$ values, and then continue with the merge steps with $k \geq 5$ . A MATLAB implementation is provided as supplemental material to this paper.

4 Experiments

In this section, we compare the speed of the proposed fast algorithm with the dyadic updating method (Huo and Székely, 2016) and also with the brute force implementation. We implemented the proposed algorithm and the brute force method in MATLAB. For the dyadic updating method, we use the authors’ MATLAB implementation, in which the key step is implemented in C. Therefore, this comparison strongly favors the dyadic updating method because more than $90 %$ of its calculations are done in C instead of MATLAB according to the MATLAB code profiler. All simulations are run on a PC with Intel^® Xeon^® Gold CPU @ 2.40GHZ processor and 16GB memory running MATLAB version 9.0.0.341360 (R2016a).

The data are generated from a simple quadratic model, $y_{i} = x_{i}^{2} + ϵ_{i}$ for $i = 1, \dots, n$ , where $x_{i}$ and $ϵ_{i}$ are i.i.d. standard normal. For each sample size $n = 2^{6}, 2^{8}, 2^{10}, 2^{12}, \dots, 2^{22}$ , average running time is calculated based on 10 generated samples. The speed comparison results are presented in Table 2. The column Merge Sort 2 in the table corresponds to proposed algorithm with initial merge steps replaced by insertion sort. For a moderately sized data (for example $n = 2^{14}$ ), the dyadic update method is about 25 times faster than the brute force implementation, whereas our proposed algorithm is about 200 times faster. Since brute force implementation requires $O (n^{2})$ space, no results are available when $n \geq 2^{16}$ . For large data ( $n \geq 2^{16}$ ), proposed method is 7 to 10 times faster than dyadic updating. Greater speed increases are expected if our proposed algorithm is implemented in C. Note that the p-value of the distance covariance is typically obtained by a permutation test, which is computationally intensive. For $1, 000$ permutations and $n = 2^{14} = 16, 384$ , the brute force implementation takes around $9.5 \times 1000$ seconds ( $158$ minutes) versus $0.048 \times 1000$ seconds (0.8 minutes) or $0.024 \times 1000$ seconds (0.4 minutes) for proposed $O (n log (n))$ implementation.

=2.5ex

n	Brute Force	Dyadic Update	Merge Sort	Merge Sort 2
$2^{6}$	0.0011(0.0052)	0.0006(0.0001)	0.0004(0.0002)	0.0001(0.0001)
$2^{8}$	0.0011(0.0003)	0.0027(0.0013)	0.0011(0.0004)	0.0004(0.0001)
$2^{10}$	0.0416(0.0034)	0.0107(0.0007)	0.0031(0.0009)	0.0013(0.0003)
$2^{12}$	0.6014(0.0160)	0.0678(0.0115)	0.0126(0.0026)	0.0059(0.0022)
$2^{14}$	9.5096(0.1134)	0.3585(0.0247)	0.0479(0.0054)	0.0241(0.0036)
$2^{16}$	-	1.7816(0.0757)	0.2511(0.0190)	0.1624(0.0261)
$2^{18}$	-	9.2188(0.2154)	1.1844(0.0552)	0.8258(0.0625)
$2^{20}$	-	45.302(0.4065)	5.4248(0.1077)	4.0098(0.1737)
$2^{22}$	-	219.04(2.7388)	25.345(0.4245)	19.8976(0.4665)

Table 2:

Average time in seconds for 10 replications. The numbers in parentheses are the estimated standard deviation of running time. No results are available for brute force implementation when

n \geq 2^{16} = 65536

because the required memory exceeds the maximum memory size. Column Merge Sort 2 corresponds to proposed algorithm with initial merge steps replaced by insertion sort.

5 Conclusions

In this paper, we presented an $O (n log (n))$ algorithm for calculating the sample distance covariance for univariate variables. For multivariate random variables, random projection approach can be adopted (Huang and Huo, 2017), which depends on calculation of distance covariance for univariate variables. The proposed algorithm is intuitive and simple to implement. Empirical results show that it outperforms existing methods in the literature. Our algorithm will speed up any computational technique that depends on the calculation of univariate distance covariance for example, feature screening (Li et al., 2012).

The proposed faster algorithm provides a tool for scientists to explore complicated dependence structures using distance covariance in larger data than what was previously possible.

References

Cormen et al. (2001) Cormen, T. H., Leiserson, C. E., Rivest, R. L., and Stein, C. (2001). Introduction to algorithms second edition.
Ginat (2004) Ginat, D. (2004). Do Senior CS Students Capitalize on Recursion? SIGCSE Bull., 36(3):82–86.
Gretton et al. (2005) Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005). Measuring statistical dependence with hilbert-schmidt norms. In International Conference on Algorithmic Learning Theory, pages 63–77. Springer.
Gretton et al. (2008) Gretton, A., Fukumizu, K., Teo, C. H., Song, L., Schölkopf, B., and Smola, A. J. (2008). A kernel statistical test of independence. In Advances in Neural Information Processing Systems, pages 585–592.
Guo et al. (2014) Guo, X., Zhang, Y., Hu, W., Tan, H., and Wang, X. (2014). Inferring nonlinear gene regulatory networks from gene expression data based on distance correlation. PloS one, 9(2):e87446.
Huang and Huo (2017) Huang, C. and Huo, X. (2017). A statistically and numerically efficient independence test based on random projections and distance covariance. arXiv preprint arXiv:1701.06054.
Huo and Székely (2016) Huo, X. and Székely, G. J. (2016). Fast computing for distance covariance. Technometrics, 58(4):435–447.
Knuth (1997) Knuth, D. E. (1997). The Art of Computer Programming: Sorting and Searching, volume 3. Pearson Education.
Li et al. (2012) Li, R., Zhong, W., and Zhu, L. (2012). Feature screening via distance correlation learning. Journal of the American Statistical Association, 107(499):1129–1139.
Pfister et al. (2018) Pfister, N., Bühlmann, P., Schölkopf, B., and Peters, J. (2018). Kernel-based tests for joint independence. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):5–31.
Sejdinovic et al. (2013) Sejdinovic, D., Sriperumbudur, B., Gretton, A., and Fukumizu, K. (2013). Equivalence of distance-based and rkhs-based statistics in hypothesis testing. The Annals of Statistics, pages 2263–2291.
Sferra et al. (2017) Sferra, G., Fratini, F., Ponzi, M., and Pizzi, E. (2017). Phylo_dcor: distance correlation as a novel metric for phylogenetic profiling. BMC bioinformatics, 18(1):396.
Székely and Rizzo (2009) Székely, G. J. and Rizzo, M. L. (2009). Brownian distance covariance. The annals of applied statistics, pages 1236–1265.
Székely et al. (2007) Székely, G. J., Rizzo, M. L., Bakirov, N. K., et al. (2007). Measuring and testing dependence by correlation of distances. The Annals of Statistics, 35(6):2769–2794.
Zhang et al. (2018) Zhang, Q., Filippi, S., Gretton, A., and Sejdinovic, D. (2018). Large-scale kernel methods for independence testing. Statistics and Computing, 28(1):113–130.
Zhou (2012) Zhou, Z. (2012). Measuring nonlinear dependence in time-series, a distance correlation approach. Journal of Time Series Analysis, 33(3):438–457.

MATLAB Code for Fast Distance Covariance

mygreenrgb0,0.6,0

⬇

blue function covsq = fastDcov(x,y)

mygreen%fastDcov computes distance correlation between column vectors x and y

mygreen%Sort the data by x. This does not change the answer.

n = blue length(x);

[x, Index] = blue sort(x);

y = y(Index);

mygreen%a_x is the vector of row sums of distance matrix of x

si = blue cumsum(x);

s = si(n);

a_x = (-(n-2):2:n).’␣.*x␣+␣(s␣-␣2*si);

%Compute␣Frobenius␣inner␣product␣between␣the

%distance␣matrices␣while␣finding␣indexes␣corresponding

%to␣sorted␣y␣using␣merge-sort.

%Weight␣vectors

v␣=␣[x␣y␣x.*y];

nw␣=␣size(v,␣2);

%␣The␣columns␣of␣idx␣are␣buffers␣to␣store␣sort␣indices␣and␣output␣buffer

%␣of␣merge-sort

%␣we␣alternate␣between␣them␣and␣avoid␣uneccessary␣copying

idx␣=␣zeros(n,␣2);

idx(:,1)␣=␣1:n;

%␣iv1,␣iv2,␣iv3,␣iv4␣are␣used␣to␣store␣sum␣of␣weights.

%␣On␣output,␣for␣

1 \leq j \leq n

%␣

iv1 (j) = \sum \begin{matrix} i < j, y_{i} < y_{j} \end{matrix} 1

%␣

iv2 (j) = \sum \begin{matrix} i < j, y_{i} < y_{j} \end{matrix} x_{i}

%␣

iv3 (j) = \sum \begin{matrix} i < j, y_{i} < y_{j} \end{matrix} y_{i}

%␣

iv4 (j) = \sum \begin{matrix} i < j, y_{i} < y_{j} \end{matrix} x_{i} y_{i}

iv1␣=␣zeros(n,1);iv2␣=␣zeros(n,1);iv3␣=␣zeros(n,1);iv4␣=␣zeros(n,1);

%␣The␣merge␣sort␣loop.

i␣=␣1;␣r␣=␣1;␣s␣=␣2;

while␣i␣<␣n;

␣␣␣␣gap␣=␣2*i;

␣␣␣␣k␣=␣0;

␣␣␣␣idx_r␣=␣idx(:,r);

␣␣␣␣csumv␣=␣[zeros(1,␣nw);␣cumsum(v(idx_r,:))];

␣␣␣␣for␣j␣=␣1:gap:n;

␣␣␣␣␣␣␣␣st1␣=␣j;␣e1␣=␣min(st1␣+␣i␣-␣1,n);

␣␣␣␣␣␣␣␣st2␣=␣j␣+␣i;␣e2␣=␣min(st2␣+␣i␣-␣1,n);

␣␣␣␣␣␣␣␣while␣(st1␣<=␣e1)␣&&␣(st2␣<=␣e2);

␣␣␣␣␣␣␣␣␣␣␣␣k␣=␣k␣+1;

␣␣␣␣␣␣␣␣␣␣␣␣idx1␣=␣idx_r(st1);

␣␣␣␣␣␣␣␣␣␣␣␣idx2␣=␣idx_r(st2);

␣␣␣␣␣␣␣␣␣␣␣␣if␣y(idx1)␣>=␣y(idx2);

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣idx(k,s)␣=␣idx1;

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣st1␣=␣st1␣+␣1;

␣␣␣␣␣␣␣␣␣␣␣␣else

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣idx(k,s)␣=␣idx2;

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣st2␣=␣st2␣+␣1;

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣iv1(idx2,␣1)␣=␣iv1(idx2)␣+␣e1␣-st1␣+1;

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣iv2(idx2)␣=␣iv2(idx2)␣+␣(csumv(e1+1,␣1)␣-␣csumv(st1,␣1));

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣iv3(idx2)␣=␣iv3(idx2)␣+␣(csumv(e1+1,␣2)␣-␣csumv(st1,␣2));

␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣␣iv4(idx2)␣=␣iv4(idx2)␣+␣(csumv(e1+1,␣3)␣-␣csumv(st1,␣3));

␣␣␣␣␣␣␣␣␣␣␣␣end;

␣␣␣␣␣␣␣␣end;

␣␣␣␣␣␣␣␣if␣st1␣<=␣e1;

␣␣␣␣␣␣␣␣␣␣␣␣kf␣=␣k␣+␣e1␣-␣st1␣+␣1;

␣␣␣␣␣␣␣␣␣␣␣␣idx((k+1):kf,␣s)␣=␣idx_r(st1:e1,:);

␣␣␣␣␣␣␣␣␣␣␣␣k␣=␣kf;

␣␣␣␣␣␣␣␣elseif␣st2␣<=e2;

␣␣␣␣␣␣␣␣␣␣␣␣kf␣=␣k␣+␣e2␣-␣st2␣+␣1;

␣␣␣␣␣␣␣␣␣␣␣␣idx((k+1):kf,␣s)␣=␣idx_r(st2:e2,:);

␣␣␣␣␣␣␣␣␣␣␣␣k␣=␣kf;

␣␣␣␣␣␣␣␣end;

␣␣␣␣end;

␣␣␣␣i␣=␣gap;

␣␣␣␣r␣=␣3-r;␣s␣=␣3-s;

end;

%␣d␣is␣the␣Frobenius␣inner␣product␣of␣the␣distance␣matrices

covterm␣=␣␣n*(x␣-␣mean(x)).’␣*␣(y␣-␣mean(y));

c1␣=␣iv1.’␣*␣v(:,␣3);

c2␣=␣sum(iv4);

c3␣=␣iv2.’␣*␣y;

c4␣=␣iv3.’␣*␣x;

d␣=␣4*((c1␣+␣c2)␣-␣(c3␣+␣c4))␣-␣2*covterm;

%␣b_y␣is␣the␣vector␣of␣row␣sums␣of␣distance␣matrix␣of␣y

ySorted␣=␣y(idx(n:-1:1,␣r));

si␣=␣cumsum(ySorted);

s␣=␣si(n);

b_y␣=␣zeros(n,␣1);

b_y(idx(n:-1:1,␣r))␣=␣(-(n-2):2:n).’␣.*ySorted␣+␣(s␣-␣2*si);

%covsq␣equals␣

V_{n}^{2} (x, y)

␣the␣square␣of␣the␣distance␣covariance

%between␣x␣and␣y.

nsq␣=␣n*n;

ncb␣=␣nsq*n;

nq␣=␣ncb*n;

term1␣=␣d␣/␣nsq;

term2␣=␣2*␣(a_x.’␣*␣b_y)␣/␣ncb;

term3␣=␣sum(a_x)␣*␣sum(b_y)␣/␣nq;

covsq␣=␣(term1␣+␣term3)␣-␣term2;

A fast algorithm for computing distance correlation

Authors

Detecting independence of random vectors II. Distance multivariance and Gaussian multivariance

Distance Metrics for Measuring Joint Dependence with Application to Causal Inference

Estimating Feature-Label Dependence Using Gini Distance Statistics

Mutual Dependence: A Novel Method for Computing Dependencies Between Random Variables

Estimating The Proportion of Signal Variables Under Arbitrary Covariance Dependence

Distance-based and RKHS-based Dependence Metrics in High Dimension

Detecting independence of random vectors I. Generalized distance covariance and Gaussian covariance

1 Introduction

2 Some Preliminaries

3 Fast Algorithm for Distance Covariance

3.1 Fast computation of the $a_{i .}$ and $b_{i .}$

3.2 Fast computation of D

Theorem 3.1

3.3 Fast computation of $d_{i}$ where $d_{i} = \sum_{j \in U_{i}} t_{j}$

4 Experiments

5 Conclusions

References

MATLAB Code for Fast Distance Covariance

A fast algorithm for computing distance correlation

Authors

Related Research

Detecting independence of random vectors II. Distance multivariance and Gaussian multivariance

Distance Metrics for Measuring Joint Dependence with Application to Causal Inference

Estimating Feature-Label Dependence Using Gini Distance Statistics

Mutual Dependence: A Novel Method for Computing Dependencies Between Random Variables

Estimating The Proportion of Signal Variables Under Arbitrary Covariance Dependence

Distance-based and RKHS-based Dependence Metrics in High Dimension

Detecting independence of random vectors I. Generalized distance covariance and Gaussian covariance

This week in AI

1 Introduction

2 Some Preliminaries

3 Fast Algorithm for Distance Covariance

3.1 Fast computation of the ai. and bi.

3.2 Fast computation of D

Theorem 3.1

3.3 Fast computation of di where di=∑j∈Uitj

4 Experiments

5 Conclusions

References

MATLAB Code for Fast Distance Covariance

3.1 Fast computation of the $a_{i .}$ and $b_{i .}$

3.3 Fast computation of $d_{i}$ where $d_{i} = \sum_{j \in U_{i}} t_{j}$