Comparing Earth Mover’s Distance for Masked Distributions

Cell-wise metrics — mean absolute error, KL divergence, chi-squared — treat two grid cells 40 km apart as exactly as different as two adjacent cells. For masked maps that is the wrong geometry, and it produces the characteristic failure where a mask that displaced everything by one cell scores worse than one that moved a tenth of the mass across the city. Earth mover’s distance is the metric that knows where the cells are.

Core Calculation # permalink

EMD — the 1-Wasserstein distance — is the minimum total work to turn one distribution into the other, where work is mass times the distance it travels:

subject to the transport plan $\gamma$ being a valid coupling:

With $d$ the ground distance in metres and $P, Q$ normalised to sum to 1, the result is a mean displacement in metres — the average distance the mask moved a unit of population. That interpretability is most of why the metric is worth the compute.

For a one-dimensional or a monotone-ordered problem, EMD collapses to the $L_1$ distance between cumulative distributions:

which is why 1-D marginals are cheap and the full 2-D grid is not.

Worked numeric example # permalink

A 3-cell toy strip, cells 500 m apart, true shares $P = (0.6, 0.3, 0.1)$.

Mask A shifts everything one cell right: $Q_A = (0.1, 0.6, 0.3)$. Mask B moves 10 % of cell 1’s mass to the far cell: $Q_B = (0.5, 0.3, 0.2)$.

Cell-wise MAE:

MAE says B is five times better. Now EMD, with the optimal plan:

EMD agrees on the ordering here but not on the magnitude — B is 2.5× better, not 5×. Now change mask B to move that 10 % across 20 km instead:

MAE is unchanged at 0.067, because the cell counts are identical; only the geography moved. EMD reports an eight-fold degradation. That divergence — same cell-wise error, wildly different spatial error — is the case EMD exists for, and it is exactly the case a masking bug produces.

Python Implementation # permalink

from __future__ import annotations

import numpy as np
from scipy.optimize import linprog
from scipy.stats import wasserstein_distance


def emd_grid(p: np.ndarray, q: np.ndarray, centroids: np.ndarray,
             normalise: bool = True) -> float:
    """Exact EMD between two grid distributions, in the units of `centroids`.

    Solves the transport LP directly. Exact and O(n^3) in the worst case, so
    it is for grids up to a few thousand cells; beyond that use the Sinkhorn
    approximation or the sliced form below.
    """
    p = np.asarray(p, dtype=float).ravel()
    q = np.asarray(q, dtype=float).ravel()
    if normalise:
        p, q = p / p.sum(), q / q.sum()
    if not np.isclose(p.sum(), q.sum()):
        raise ValueError("distributions must carry equal mass; normalise first")

    n, m = len(p), len(q)
    d = np.linalg.norm(centroids[:, None, :] - centroids[None, :, :], axis=2)

    # Marginal constraints: rows sum to p, columns sum to q.
    a_eq = np.zeros((n + m, n * m))
    for i in range(n):
        a_eq[i, i * m:(i + 1) * m] = 1.0
    for j in range(m):
        a_eq[n + j, j::m] = 1.0

    res = linprog(d.ravel(), A_eq=a_eq, b_eq=np.concatenate([p, q]),
                  bounds=(0, None), method="highs")
    if not res.success:
        raise RuntimeError(f"transport LP failed: {res.message}")
    return float(res.fun)


def emd_sliced(p: np.ndarray, q: np.ndarray, centroids: np.ndarray,
               n_projections: int = 128, seed: int = 0) -> float:
    """Sliced EMD: average 1-D Wasserstein over random directions.

    Linear in the number of cells and within a few percent of exact on
    realistic urban grids, which makes it the one to run inside a CI gate.
    """
    rng = np.random.default_rng(seed)
    p = np.asarray(p, dtype=float).ravel()
    q = np.asarray(q, dtype=float).ravel()
    p, q = p / p.sum(), q / q.sum()

    total = 0.0
    for _ in range(n_projections):
        theta = rng.uniform(0, np.pi)
        proj = centroids @ np.array([np.cos(theta), np.sin(theta)])
        total += wasserstein_distance(proj, proj, u_weights=p, v_weights=q)
    return total / n_projections

The sliced form is the practical default. It runs in milliseconds on a 128 × 128 grid where the exact LP takes minutes, and its bias is a consistent slight underestimate — acceptable when the gate is a threshold rather than a published figure.

Verification # permalink

Check the metric axioms hold on your implementation. Three assertions catch nearly every bug:

def check_emd(centroids: np.ndarray, seed: int = 0) -> dict:
    rng = np.random.default_rng(seed)
    n = len(centroids)
    p = rng.dirichlet(np.ones(n))
    q = rng.dirichlet(np.ones(n))
    r = rng.dirichlet(np.ones(n))

    return {
        # Identity: no work is needed to turn a distribution into itself.
        "identity": np.isclose(emd_grid(p, p, centroids), 0.0, atol=1e-9),
        # Symmetry: transport is reversible at the same cost.
        "symmetry": np.isclose(emd_grid(p, q, centroids),
                               emd_grid(q, p, centroids), rtol=1e-6),
        # Triangle inequality: a metric, not merely a divergence.
        "triangle": (emd_grid(p, r, centroids) <=
                     emd_grid(p, q, centroids) + emd_grid(q, r, centroids) + 1e-6),
        # Scale: doubling every ground distance doubles the EMD exactly.
        "scale_linear": np.isclose(emd_grid(p, q, centroids * 2),
                                   2 * emd_grid(p, q, centroids), rtol=1e-6),
    }

The scale check is the one that catches the most common real error — computing distances in degrees of latitude and longitude rather than metres, which makes the result uninterpretable and, away from the equator, anisotropic. Project to a local metric CRS before taking centroid distances.

Calibrate the threshold against a known-good release. An EMD of 340 m means nothing in isolation. Compute it for a release everyone agrees is fine, and set the gate at some multiple of that. A threshold set from first principles will be either unreachable or vacuous.

Confirm the sliced approximation against exact on a subsample. Take a 32 × 32 crop, compute both, and record the ratio. If sliced is within 5 % of exact on the crop, the full-grid sliced figure is trustworthy as a gate; if not, raise the projection count.

Edge Cases & Adjustments # permalink

• Unequal mass. Suppression removes population, so $P$ and $Q$ often do not sum to the same total. Normalising hides that loss entirely — a release that dropped 20 % of the population can have a small EMD. Report the mass difference as a separate number alongside the normalised EMD, or use partial-transport EMD with an explicit penalty for unmatched mass.
• Degrees versus metres. Always project first. At 55° latitude a degree of longitude is 64 km and a degree of latitude is 111 km; an unprojected EMD is a weighted average of two different units.
• Grid alignment. Comparing distributions on different grids requires a common support. Resample both to the coarser grid; resampling to the finer one invents structure and depresses the distance.
• The metric is not privacy. A small EMD says the map is useful, not that it is safe. Utility metrics and risk metrics are independent gates and a release must clear both — a perfect-utility release is the raw data.
• Very large grids. Beyond about 10⁴ cells the exact LP is impractical. Use the sliced form for gating and the exact form on the worst-scoring region for diagnosis.

FAQ # permalink

Why not just use KL divergence?

KL is infinite whenever the released map has a zero where the truth does not — which suppression guarantees. It is also blind to geography: moving mass one cell and moving it across the city cost the same.

What EMD value is acceptable?

Whatever your calibration says. As a rule of thumb, an EMD below one grid-cell width means the mask is doing local rounding, and an EMD above three cell widths means it is relocating population. Both statements are about your cell size, not about EMD.

Does EMD replace cell-wise metrics?

No — it complements them. Cell-wise error catches a mask that changes counts without moving them; EMD catches one that moves them without changing totals. Report both.

Is the sliced approximation good enough for a published figure?

For a gate, yes. For a number in a methodology document, compute the exact value once and cite that, noting the sliced figure is what the pipeline monitors.

Related # permalink

• Utility Preservation Metrics for Masked Maps — the metric family this sits in
• Measuring Utility Loss in DP Heatmaps — the cell-wise companion metrics
• Plotting the Privacy–Utility Frontier for Spatial Queries — putting EMD on the y-axis
• Geospatial Masking & Perturbation Techniques — the masks being measured
• Validating Synthetic Trajectories Against Real Mobility Metrics — distribution comparison in the synthetic case

← Back to Utility Preservation Metrics for Masked Maps