Distance Covariance From Euclidean Distances

Why Characteristic Functions Lead to Distances

A characteristic-function discrepancy compares the joint distribution of two random vectors with the product of their marginal distributions. Its frequency-domain form is useful for defining a broad measure of association, but an integral over frequencies can be awkward to compute directly. For Euclidean spaces, that same measure has an expression in terms of pairwise distances. This tutorial derives the connection and shows how to calculate it from paired observations.

Let $X$ take values in $\mathbb{R}^p$ and $Y$ in $\mathbb{R}^q$, where $p,q\geq 1$. Write $\varphi_{X,Y}(t,s)=E[\exp\{i(t^\top X+s^\top Y)\}]$ for their joint characteristic function, and $\varphi_X(t)$, $\varphi_Y(s)$ for their marginal characteristic functions. We use the Euclidean norm, denoted by $\lVert\cdot\rVert$. A convenient sufficient condition for the distance expectations and derivation below is $E\lVert X\rVert<\infty$, $E\lVert Y\rVert<\infty$, and $E[\lVert X\rVert\lVert Y\rVert]<\infty$.

Definition: The squared distance covariance of $X$ and $Y$ is $$ \mathcal{V}^2(X,Y)=\frac{1}{c_p c_q}\int_{\mathbb{R}^p}\int_{\mathbb{R}^q} \frac{\left|\varphi_{X,Y}(t,s)-\varphi_X(t)\varphi_Y(s)\right|^2} {\lVert t\rVert^{p+1}\lVert s\rVert^{q+1}}\,ds\,dt, \qquad c_d=\frac{\pi^{(d+1)/2}}{\Gamma((d+1)/2)}. $$ The integrand at $t=0$ or $s=0$ can be assigned any value; these sets have Lebesgue measure zero. Distance covariance itself is the nonnegative square root $\mathcal{V}(X,Y)$.

The weights emphasize low frequencies and depend on the dimensions $p$ and $q$. The constants $c_p,c_q$ are chosen so that the frequency integral translates exactly into products of ordinary Euclidean distances. The squared quantity is the one most directly expressed as an inner product of distance matrices.

The Euclidean-Distance Identity

The key Fourier identity is that a norm can be recovered by integrating $1-\cos(t^\top x)$ against the same type of weight used in the definition.

Lemma: For $x\in\mathbb{R}^d$, $$ \int_{\mathbb{R}^d}\frac{1-\cos(t^\top x)}{\lVert t\rVert^{d+1}}\,dt=c_d\lVert x\rVert. $$

Here is a derivation of the constant as well as the dependence on $x$. Rotational invariance lets us rotate $x$ onto the first coordinate axis, and scaling $t$ then shows that the integral is proportional to $\lVert x\rVert$. For a unit vector in the first coordinate direction, write $t=(u,v)$, with $v\in\mathbb{R}^{d-1}$. Integrating over $v$ gives $\int_{\mathbb{R}^{d-1}}(u^2+\lVert v\rVert^2)^{-(d+1)/2}\,dv =\pi^{(d-1)/2}|u|^{-2}/\Gamma((d+1)/2)$. The remaining integral is $\int_{\mathbb{R}}(1-\cos u)u^{-2}\,du=\pi$. Their product is $c_d$, proving the identity. In dimension one, the same calculation uses no $v$-coordinates.

From the Weighted Integral to Expected Distances

Take three independent copies $(X,Y),(X',Y'),(X'',Y'')$ of the paired vector $(X,Y)$. Marginal components taken from different copies are independent. Under the moment conditions stated above, the population quantity also has this distance representation:

Theorem: $$ \mathcal{V}^2(X,Y) =E[\lVert X-X'\rVert\lVert Y-Y'\rVert] +E\lVert X-X'\rVert\,E\lVert Y-Y'\rVert -2E[\lVert X-X'\rVert\lVert Y-Y''\rVert]. $$

To derive it, apply the Euclidean-distance identity separately in the $t$ and $s$ integrals. A product of distances becomes a double integral of products of terms of the form $1-\cos(t^\top x)$ and $1-\cos(s^\top y)$. Expanding those products gives three types of terms: joint differences from two paired copies, independent marginal differences, and the cross term involving $X-X'$ and $Y-Y''$. Their coefficients are respectively $1$, $1$, and $-2$, as in the displayed formula.

One can check the characteristic-function algebra that identifies those terms directly. The squared discrepancy expands to $\lvert\varphi_{X,Y}\rvert^2+\lvert\varphi_X\rvert^2\lvert\varphi_Y\rvert^2 -2\operatorname{Re}(\varphi_{X,Y}\overline{\varphi_X}\,\overline{\varphi_Y})$. The three terms correspond to the paired, independent-marginal, and cross terms just described. Since the frequency weights are even in $t$ and $s$, the cosine expansions can be averaged over sign changes without changing the integral. The distance identity then gives the stated coefficients and distances. The finite moment conditions ensure the distance expectations are finite and justify the integrations; equivalently, the calculation can first be made on truncated frequency domains and then passed to the limit.

The cross expectation uses $Y''$, not $Y'$: $X$ and $Y$ in $X-X'$ and $Y-Y''$ belong to the same copy, while the subtracted components come from independent copies. This pairing is what retains the dependence structure in the formula.

Double-Centering the Sample Distances

Now observe $n$ paired vectors $(x_i,y_i)$, with $x_i\in\mathbb{R}^p$ and $y_i\in\mathbb{R}^q$. Form the Euclidean distance matrices $a_{ij}=\lVert x_i-x_j\rVert$ and $b_{ij}=\lVert y_i-y_j\rVert$. Define their row means, column means, and grand means in the usual way. Since distance matrices are symmetric, row and column means agree, but keeping both in the definition makes the centering operation explicit.

Definition: The doubly centered distance entries are $$ A_{ij}=a_{ij}-\frac{1}{n}\sum_{k=1}^n a_{ik} -\frac{1}{n}\sum_{k=1}^n a_{kj} +\frac{1}{n^2}\sum_{k=1}^n\sum_{\ell=1}^n a_{k\ell}, $$ with $B_{ij}$ defined by the same formula using $b$. The empirical squared distance covariance, also called the biased sample distance covariance squared, is $$ \mathcal{V}_n^2=\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n A_{ij}B_{ij}. $$

This is the distance-matrix version of the weighted characteristic-function definition applied to the empirical distribution that puts mass $1/n$ on each observed pair. In particular, the formula pairs the two centered distances at the same indices $i,j$. It is not obtained by centering each coordinate of $x_i$ and $y_i$ and then taking ordinary covariance.

Theorem: For any paired observations, the doubly centered formula satisfies $$ \mathcal{V}_n^2 =\frac{1}{n^2}\sum_{i,j}a_{ij}b_{ij} +\left(\frac{1}{n^2}\sum_{i,j}a_{ij}\right) \left(\frac{1}{n^2}\sum_{i,j}b_{ij}\right) -\frac{2}{n^3}\sum_{i=1}^n \left(\sum_{j=1}^n a_{ij}\right)\left(\sum_{j=1}^n b_{ij}\right). $$

For the proof, let $\bar a_{i\cdot}=n^{-1}\sum_j a_{ij}$, $\bar a_{\cdot j}=n^{-1}\sum_i a_{ij}$, and $\bar a_{\cdot\cdot}=n^{-2}\sum_{i,j}a_{ij}$, and define the corresponding quantities for $b$. Expand $A_{ij}B_{ij}$ and sum over $i,j$. The $a_{ij}b_{ij}$ terms give the first term after division by $n^2$. Terms involving a grand mean give $\bar a_{\cdot\cdot}\bar b_{\cdot\cdot}$. The row-mean and column-mean cross terms sum to $-2n^{-1}\sum_i\bar a_{i\cdot}\bar b_{i\cdot}$, because both distance matrices are symmetric. Substituting the definitions of these means yields the displayed expression. Thus the expression is exactly the inner product of the centered matrices, not a separate approximation. $\square$

Although the expanded formula contains a negative term, the total is nonnegative: the empirical weighted characteristic-function integral is an integral of a squared modulus, and the double-centering formula equals that integral. The individual summands in the expanded expression need not be nonnegative.

Worked Calculations

Worked Example: Three Paired Scalar Observations

Take $x=(0,1,2)$ and $y=(0,1,0)$, with each coordinate pair representing one observation. The distance matrices are $$ (a_{ij})=\begin{pmatrix}0&1&2\\1&0&1\\2&1&0\end{pmatrix}, \qquad (b_{ij})=\begin{pmatrix}0&1&0\\1&0&1\\0&1&0\end{pmatrix}. $$

For $a$, the row means are $1,2/3,1$, and the grand mean is $8/9$. Double-centering gives $$ A=\frac{1}{9}\begin{pmatrix}-10&2&8\\2&-4&2\\8&2&-10\end{pmatrix}. $$ For $b$, the row means are $1/3,2/3,1/3$, and the grand mean is $4/9$, so $$ B=\frac{1}{9}\begin{pmatrix}-2&4&-2\\4&-8&4\\-2&4&-2\end{pmatrix}. $$

The sum of entrywise products of the integer matrices in these displays is $20+8-16+8+32+8-16+8+20=72$. Therefore $\sum_{i,j}A_{ij}B_{ij}=72/81=8/9$, and $\mathcal{V}_3^2=(8/9)/9=8/81$. The centering may produce negative entries, but their paired products, summed over all indices, give the nonnegative statistic.

Worked Example: Two Observations and Scale

Let $x_1=(0,0)$, $x_2=(2,0)$, $y_1=(0)$, and $y_2=(3)$. Then $$ (a_{ij})=\begin{pmatrix}0&2\\2&0\end{pmatrix}, \qquad (b_{ij})=\begin{pmatrix}0&3\\3&0\end{pmatrix}. $$ Both grand means equal their row means, so double-centering gives $A=\begin{pmatrix}-1&1\\1&-1\end{pmatrix}$ and $B=\begin{pmatrix}-3/2&3/2\\3/2&-3/2\end{pmatrix}$.

Each of the four entrywise products is $3/2$. Their sum is $6$, and division by $n^2=4$ gives $\mathcal{V}_2^2=3/2$. This also checks the units: multiplying all $x$-distances by a positive factor $r$ and all $y$-distances by $s$ multiplies squared distance covariance by $rs$. The square root $\mathcal{V}_n$ consequently scales by $\sqrt{rs}$.

Worked Example: Paired Vectors in Two Dimensions

Consider two paired observations $x_1=(0,0)$, $x_2=(3,4)$, $y_1=(1,1)$, and $y_2=(1,13)$. Their within-variable distances are $\lVert x_1-x_2\rVert=\sqrt{3^2+4^2}=5$ and $\lVert y_1-y_2\rVert=12$. Hence the distance matrices have off-diagonal entries $5$ and $12$, respectively, and zero diagonal entries.

With two observations, double-centering turns a distance matrix with off-diagonal value $d$ into $\begin{pmatrix}-d/2&d/2\\d/2&-d/2\end{pmatrix}$. Thus $A$ has entries $-5/2,5/2,5/2,-5/2$, and $B$ has entries $-6,6,6,-6$. Every paired product equals $15$, so their sum is $60$, and $\mathcal{V}_2^2=60/4=15$. The directions of the vectors do not enter separately: the construction uses their Euclidean distances.

Common Pitfalls and Why the Formulation Matters

A frequent error is to subtract only row means from each distance matrix, or to center the original observations rather than the distances. The required operation is double-centering the pairwise distances: subtract the row and column means and add the grand mean. Another error is to divide the sum of centered products by $n$ instead of $n^2$. The formula here is the empirical-distribution, or biased, version, with exactly $n^2$ ordered pairs, including the diagonal.

The distance representation matters because it turns a frequency-weighted comparison into a calculation using observed Euclidean geometry. It applies when the two variables are vectors of different dimensions, and it is unchanged by translating all observations within either variable: pairwise distances are unchanged by translation. The result is a scalar summary of paired multivariate observations; interpreting what its value implies about dependence requires care and is not supplied by the distance-matrix algebra alone.

Key takeaway: Under suitable moment conditions, the weighted characteristic-function discrepancy is a combination of expected products of Euclidean distances. For paired data, its empirical form is the $n^{-2}$-scaled inner product of the two doubly centered distance matrices.

Check Your Understanding

Use the definitions and calculations above to answer the following questions.

What role do $c_p$ and $c_q$ play in the weighted characteristic-function definition?
Why does the population distance formula require a cross term involving $Y''$ from an independent copy?
For a symmetric distance matrix, how are its row means and column means related?
Why can centered distance entries be negative even though the empirical squared distance covariance is nonnegative?
If all $x_i$ are multiplied by $r>0$ and all $y_i$ by $s>0$, by what factor does $\mathcal{V}_n^2$ change?