Energy Distance

Energy distance compares two distributions using nothing but expected pairwise distances between draws.

\[ D_E(P, Q) = 2\,\mathbb{E}\|X - Y\| - \mathbb{E}\|X - X’\| - \mathbb{E}\|Y - Y’\| \]

Here \(X,X^{\prime}\) are independent draws from P and \(Y,Y^{\prime}\) from Q. With finite first moments, the displayed Euclidean quantity is nonnegative and zero exactly for equal distributions. This article uses the squared-energy convention; SciPy’s energy_distance returns its square root in one dimension. Distinguish that convention when comparing units or values with MMD.

It is an MMD with the kernel already chosen

Sejdinovic et al. proved that for a semimetric of negative type there is always a kernel that generates it, and that \( D_{E,\rho}(P,Q) = 2\gamma_k^2(P,Q) \) — the energy distance is exactly twice the squared MMD of that kernel. For the Euclidean distance the kernel, centered at the origin, is \( k(x,y) = \tfrac{1}{2}(\|x\| + \|y\| - \|x - y\|) \).

import numpy as np
from scipy.stats import energy_distance


def pdist(A, B):
    return np.sqrt(np.maximum(((A[:, None, :] - B[None, :, :]) ** 2).sum(-1), 0.0))


def energy(X, Y):
    return 2 * pdist(X, Y).mean() - pdist(X, X).mean() - pdist(Y, Y).mean()


def mmd2_distance_kernel(X, Y):
    Z = np.vstack([X, Y])
    nz = np.sqrt((Z ** 2).sum(1))
    K = 0.5 * (nz[:, None] + nz[None, :] - pdist(Z, Z))
    m = len(X)
    return K[:m, :m].mean() + K[m:, m:].mean() - 2 * K[:m, m:].mean()


rng = np.random.default_rng(77)
print("%-14s %16s %16s %16s" % ("data", "energy distance", "2 x MMD^2", "ratio"))
for label, X, Y in (
        ("1-D normals", rng.normal(size=(400, 1)), rng.normal(0.4, 1.3, (500, 1))),
        ("4-D normals", rng.normal(size=(400, 4)), rng.normal(0.3, 1.0, (500, 4))),
        ("mixed scales", rng.normal(0, 20, (300, 3)), rng.normal(2, 25, (350, 3)))):
    e, k = energy(X, Y), 2 * mmd2_distance_kernel(X, Y)
    print("%-14s %16.10f %16.10f %16.13f" % (label, e, k, e / k))
x = rng.normal(size=600)
y = rng.normal(0.4, 1.3, 700)
e = energy(x[:, None], y[:, None])
print("sqrt(energy)                %.10f" % np.sqrt(e))
print("scipy.stats.energy_distance %.10f" % energy_distance(x, y))

data            energy distance        2 x MMD^2            ratio
1-D normals        0.1317430682     0.1317430682  1.0000000000000
4-D normals        0.0898299056     0.0898299056  1.0000000000000
mixed scales       0.8346635368     0.8346635368  1.0000000000000
sqrt(energy)                0.3353322152
scipy.stats.energy_distance 0.3353322152

The two routes agree to the last digit the format prints, in one dimension, in four, and on data whose scale is twenty times larger. The claim that energy distance needs no kernel is really a claim that its kernel is fixed rather than absent, and that is the honest summary of the trade: no bandwidth to tune, and no bandwidth to tune toward the alternative you care about.

It carries the units of the data

\( \|x - y\| \) has the dimensions of \( x \). Multiply every value by a constant c and the whole expression multiplies by |c|, which makes an absolute cut-off a statement about the column’s unit rather than about the customers in it. The next block records the same simulated customer ages in years, months, and days. It takes the energy distance from energy above applied to one-dimensional columns, the KS statistic from scipy.stats.ks_2samp, and a p-value from a permutation test: shuffle the pooled ages, split them back into two samples, and count how often the shuffled energy distance reaches the observed one, with the observed arrangement counted once. It continues the same session.

from scipy.stats import ks_2samp


def energy_1d(x, y):
    return energy(x[:, None], y[:, None])


def perm_p(x, y, reps=500):
    stat = energy_1d(x, y)
    z, m, rng, hits = np.concatenate([x, y]), len(x), np.random.default_rng(9), 0
    for _ in range(reps):
        rng.shuffle(z)
        hits += energy_1d(z[:m], z[m:]) >= stat
    return (hits + 1) / (reps + 1)


rng = np.random.default_rng(404)
base_years = 30 + 9 * rng.normal(size=500)
cur_years = 33 + 9 * rng.normal(size=500)

print("%-22s %14s %10s %10s %10s"
      % ("customer age recorded in", "energy", "sqrt", "KS D", "energy p"))
for unit, factor in (("years", 1.0), ("months", 12.0), ("days", 365.25)):
    x, y = base_years * factor, cur_years * factor
    e = energy_1d(x, y)
    print("%-22s %14.4f %10.4f %10.4f %10.4f"
          % (unit, e, np.sqrt(e), ks_2samp(x, y).statistic, perm_p(x, y)))

customer age recorded in         energy       sqrt       KS D   energy p
years                          0.8126     0.9014     0.1780     0.0020
months                         9.7510     3.1227     0.1780     0.0020
days                         296.7954    17.2278     0.1780     0.0020

The same simulated customers, the same drift, three ways of recording an age: 0.8126, 9.7510, 296.7954. The KS reading stays at 0.1780 in every row and the permutation p-value stays at 0.0020. A guide that says 0.15 means similar and 1.2 means different is calibrated to whatever unit its author happened to have, and it transfers to nothing.

The same applies across features. In a raw feature vector, columns with larger scales can dominate Euclidean distances. Choose scaling or weights to reflect the features’ units and the changes that matter; standardization is one option, and changes the geometry being tested. Assess sampling evidence with a permutation test under exchangeability rather than treating the distance alone as significance.

Area, where KS takes height

In one dimension the definition collapses into something very familiar: \( D_E = 2\int (F - G)^2 \). Under this convention it is twice the squared \( L^2 \) distance between the two cumulative curves, and the KS statistic is the \( L^\infty \) distance between the same two curves. Same object, different norm. The next block evaluates both from the same CDF gap on a fine grid over [0, 1]: the largest absolute gap is the KS statistic, and twice the sum of squared gaps times the grid step approximates the integral, to be compared with energy_1d from the previous block.

def cdf_gap(x, y, grid):
    x, y = np.sort(x), np.sort(y)
    return (np.searchsorted(x, grid, "right") / x.size
            - np.searchsorted(y, grid, "right") / y.size)


rng = np.random.default_rng(2026)
u = rng.random(3000)
g = 0.15
scenarios = {
    "narrow: 15% of the mass slid one notch": np.where((u >= 0.40) & (u < 0.40 + g), u + g, u),
    "broad: every value pulled upward": u ** (1 / 1.5),
}
grid = np.linspace(0, 1, 20001)
step = grid[1] - grid[0]

print("%-40s %8s %10s %14s" % ("drift", "KS D", "energy", "2 * integral"))
for name, cur in scenarios.items():
    gap = cdf_gap(u, cur, grid)
    print("%-40s %8.4f %10.6f %14.6f"
          % (name, np.max(np.abs(gap)), energy_1d(u, cur), 2 * np.sum(gap ** 2) * step))

drift                                        KS D     energy   2 * integral
narrow: 15% of the mass slid one notch     0.1400   0.003915       0.003915
broad: every value pulled upward           0.1617   0.024064       0.024064

The last column recomputes the energy distance as twice the integral of the squared CDF gap and lands on the same number, which is the identity rather than a coincidence. The two drifts look almost the same to KS, 0.1400 against 0.1617, and differ by a factor of six in energy distance, because one opens a narrow triangle in the CDF gap and the other lifts the whole curve.

That is the choice in one sentence. The supremum measures the largest CDF gap; the integral accumulates squared gaps across values. A narrow change and a broad change can therefore look nearly the same to the supremum and very different to the integral, as the example shows. Choose a test by evaluating power for the changes and sample sizes that matter, along with its assumptions and operational cost. Neither is a better drift monitor in the abstract.

References

• Sejdinovic, D., Sriperumbudur, B., Gretton, A., & Fukumizu, K. (2013). Equivalence of distance-based and RKHS-based statistics in hypothesis testing. The Annals of Statistics, 41(5), 2263–2291.
• SciPy documentation, scipy.stats.energy_distance