Approximate Earth Mover’s Distance (EMD)

• Approximate Earth Mover’s Distance (EMD) is a family of algorithms that estimate transportation cost between distributions without solving the full optimal transport problem.
• The methods include (1+ε)-approximations for geometric inputs and additive approximations for arbitrary cost matrices, balancing accuracy and computational efficiency.
• These approaches are applied in high-dimensional search, streaming models, indexing, and differentiable loss functions for machine learning.

Searching arXiv for the papers on arXiv to ground the article. Approximate Earth Mover’s Distance (EMD) comprises algorithmic, geometric, streaming, indexing, and learned methods for estimating the transportation cost between two distributions without solving the full optimal transport problem exactly. In its standard form, EMD is the minimum cost of moving mass from one distribution to another under a ground distance, and exact computation is often computationally intensive or memory intensive. Approximation results therefore target different regimes: $(1+\varepsilon)$-approximation for discrete or geometric inputs, additive approximation in truly-subquadratic time for arbitrary cost matrices, lower-bound relaxations for search and indexing, and differentiable surrogates for learning pipelines (Kreveld et al., 2021, Beretta et al., 2023).

1. Formal definition and approximation regimes

For two distributions of equal total mass, the Earth Mover’s Distance measures the cost of transforming one distribution into the other, where the cost of moving a unit of mass is equal to the distance over which it is moved. In the continuous geometric setting, a transport plan is a nonnegative measurable function $\lambda(p,s)$ describing how much infinitesimal mass at $p$ is sent to $s$, and the EMD is

$$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$

where $H$ is the set of all joint measures on $P\times S$ whose marginals match the given distributions (Kreveld et al., 2021).

In the discrete setting, EMD is the standard minimum-cost flow or matching problem. For weighted points one recovers

$$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$

subject to the supply and demand constraints. For histograms $P,Q$ with cost matrix $(c_{i,j})$, the EMD formulation used in additive-approximation work is

$\lambda(p,s)$0

with the usual row and column marginal constraints (Beretta et al., 2023).

Approximation objectives differ substantially across problem classes. The geometric-object algorithms give a $\lambda(p,s)$1-approximation when one side is a set of weighted points and the other side is a set of line segments, triangles or $\lambda(p,s)$2-dimensional simplices, and a $\lambda(p,s)$3-approximation with a small additive term when both sides are continuous objects such as segments, triangles, or simplices (Kreveld et al., 2021). For arbitrary $\lambda(p,s)$4 cost matrices, an additive approximation scheme estimates EMD up to a $\lambda(p,s)$5-additive approximation in truly-subquadratic time $\lambda(p,s)$6 and makes no assumptions about triangle inequality (Beretta et al., 2023). In data streams, approximation guarantees include a $\lambda(p,s)$7-approximation in the insertion-only model and an $\lambda(p,s)$8-approximation in the turnstile model under a small-support assumption (Yousefi et al., 2014).

These regimes are not interchangeable. When neither side is a point set in the continuous geometric setting, the optimal transport cost can be arbitrarily small, so purely multiplicative approximation is impossible (Kreveld et al., 2021). In high-dimensional point-set settings, fine-grained hardness results show that significant improvements in exact and approximate algorithms would contradict the Orthogonal Vectors Conjecture or the Hitting Set Conjecture (Rohatgi, 2019). This suggests that “approximate EMD” is not a single algorithmic notion but a family of problem-specific compromises between multiplicative error, additive error, dimensional dependence, and output format.

2. Approximation for continuous geometric objects

For weighted geometric objects, approximation proceeds by reducing continuous transport to discrete transport through subdivision and discretization. When $\lambda(p,s)$9 is a set of weighted points and $p$0 is a set of line segments, the algorithm sets $p$1, subdivides each segment so that distances from any fixed point vary by at most a $p$2 factor on each subsegment unless the subsegment lies within $p$3 of the point, replaces each subsegment by a representative point of equal weight, invokes a black-box $p$4-approximation for discrete EMD, and then spreads the resulting flow evenly over the original subsegments (Kreveld et al., 2021).

The key analytical structure is a lower bound by discrete flow plus a bounded discretization loss. In the plane with $p$5 distance, the minimum-cost flow on the bipartite graph between points and subdivided subsegments has cost $p$6, and the final transport plan satisfies

$p$7

Setting $p$8 yields a $p$9-approximation. The overall running time is $s$0 (Kreveld et al., 2021).

The same scheme extends to triangles and $s$1-simplices by replacing segment subdivisions with square-/cube-quadtree subdivisions. In 2D one shows $s$2, and in $s$3 dimensions

$s$4

The same lemmas yield $s$5, so setting $s$6 gives a $s$7-approximation for points to triangles or simplices (Kreveld et al., 2021).

When both sides are continuous, the main obstacle is that the optimum can be arbitrarily small. The algorithms therefore introduce a small-distance greedy match that removes parts of $s$8 and $s$9 lying within a radius threshold, incurring only a controlled additive cost, and then apply the same subdivision-to-discretization reduction to the remainder (Kreveld et al., 2021). For segments↔segments in the plane, Theorem 3.5 gives

$$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$0

so setting $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$1 yields $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$2 multiplicative plus $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$3 additive, with running time $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$4 (Kreveld et al., 2021). For triangles↔triangles, the additive term is $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$5, and for $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$6-simplices↔$$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$7-simplices the final plan has cost at most $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$8 (Kreveld et al., 2021).

A central feature of these results is that the algorithms actually calculate the transport plan, rather than just the cost, and the paper states that these are the first combinatorial algorithms with a provable approximation ratio for the Earth Mover’s Distance when the objects are continuous rather than discrete points (Kreveld et al., 2021).

3. Discrete algorithmic paradigms

Approximation algorithms for discrete EMD span several distinct paradigms: subquadratic additive schemes, primal-dual regularization, streaming sketches, nearest-neighbor heuristics with analysis, and lower-bound relaxations for massively parallel search.

A concise comparison of representative guarantees is given below.

Setting	Guarantee	Running time / space
Arbitrary $$d_e(P,S)\;=\;\inf_{\lambda\in H}\;\int_{P}\int_{S}d(p,s)\;d\lambda(p,s),$$9 cost matrix	$H$0	$H$1 (Beretta et al., 2023)
Insertion-only stream, one side has $H$2 distinct points	$H$3-approximation	$H$4 space (Yousefi et al., 2014)
Turnstile stream, one side has $H$5 distinct points	$H$6-approximation	$H$7 space (Yousefi et al., 2014)
Flux formulation with $H$8 type regularization	approximate EMD via primal-dual updates	per iteration $H$9, memory $P\times S$0 (Li et al., 2016)
NNS-EMD	practical approximation with theoretical error bound	memory $P\times S$1 (Meng et al., 2024)
OMR / ICT / ACT-$P\times S$2	lower bounds with $P\times S$3	linear-complexity GPU variants for constant $P\times S$4 (Atasu et al., 2018)

For arbitrary costs, the truly-subquadratic additive scheme adapts a classical primal–dual algorithm for minimum-weight bipartite matching and implements its main steps in sublinear time via the oracles LargeMatchingForward and AugmentEligible. After $P\times S$5 outer iterations, there are at most $P\times S$6 unmatched vertices on each side, and discarding these as outliers yields the additive approximation guarantee (Beretta et al., 2023). The method works even when the costs do not satisfy triangle inequality, and it extends to max-cardinality bipartite matching with a knapsack constraint (Beretta et al., 2023).

A different route reformulates EMD as an $P\times S$7 type minimization over a flux field: $P\times S$8 The quadratic perturbation restores uniqueness for any $P\times S$9, and a Chambolle–Pock primal–dual iteration yields closed-form local updates based on shrinkage. Per iteration cost is $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$0, memory scales like $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$1, and the overall method is described as typically near-linear in problem size (Li et al., 2016). This is an approximation in the sense that the regularized problem converges to a particular solution of the original $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$2 problem as $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$3 (Li et al., 2016).

In streaming models, the insertion-only $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$4 algorithm uses a $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$5-coreset for $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$6-median on one multiset and stores the other multiset when the number of distinct points is at most $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$7; the final answer computes EMD between the coreset and the stored set (Yousefi et al., 2014). The turnstile algorithm instead combines multiple random grid shifts, turnstile $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$8-sketches, and an $$d_e(P,S)\;=\;\min_{f_{i,j}\ge0}\; \sum_{i=1}^n\sum_{j=1}^m f_{i,j}\,d(p_i,s_j)$$9-sketch, and returns an $P,Q$0-approximation (Yousefi et al., 2014).

For large-scale discrete histograms, NNS-EMD repeatedly matches each active consumer to its nearest active supplier, assigns flow by a Random Protocol or Greedy Protocol, and peels off exhausted weights. The constructed flow is feasible, and the paper states the error bound

$P,Q$1

Memory is $P,Q$2, and GPU vectorization reduces the per-iteration nearest-neighbor step in practice (Meng et al., 2024). Lower-bound relaxations such as OMR, ICT, and ACT-$P,Q$3 instead impose weakened in-flow constraints; the ordering

$P,Q$4

is explicit, and the linear-complexity GPU implementation LC-ACT-$P,Q$5 has overall time $P,Q$6 for constant $P,Q$7 (Atasu et al., 2018).

4. High-dimensional, search, and sketch-based approximation

High-dimensional EMD has motivated reductions, embeddings, and hashing schemes that trade exact transport for subquadratic time or efficient nearest-neighbor search. A recent reduction shows that if a $P,Q$8-approximate Closest Pair can be computed in time $P,Q$9, then a $(c_{i,j})$0 approximation to EMD can be computed in time $(c_{i,j})$1. Plugging in the fastest known Closest Pair algorithm yields a randomized $(c_{i,j})$2-approximation for EMD on high-dimensional point sets in time

$(c_{i,j})$3

improving over the prior fastest running time of $(c_{i,j})$4 (Beretta et al., 9 Aug 2025).

The technical core of that reduction is a sublinear implementation of the Multiplicative Weights Update framework for EMD. The weights are never explicitly stored; instead, the algorithm stores only the current dual variables $(c_{i,j})$5 and samples implicitly from the induced constraint distribution using a Closest Pair oracle plus rejection sampling (Beretta et al., 9 Aug 2025). This suggests that geometric structure can replace explicit quadratic-sized flow representations in high dimension.

For nearest-neighbor search under $(c_{i,j})$6, data-dependent locality sensitive hashing improves the approximation from the data-independent $(c_{i,j})$7 of Andoni–Indyk–Krauthgamer to $(c_{i,j})$8. The final theorem states that for any distribution $(c_{i,j})$9 on $\lambda(p,s)$00, there is a random hash with approximation factor

$\lambda(p,s)$01

for constant $\lambda(p,s)$02, computable in $\lambda(p,s)$03 time (Jayaram et al., 2024). The same paper states that the resulting data-dependent LSH also gives optimal (distributional) sketches for the Earth Mover’s Distance, and that the approximation is optimal up to $\lambda(p,s)$04 factors among schemes that collide close points with constant probability (Jayaram et al., 2024).

In bounded-doubling metrics, approximation can be shifted to preprocessing and query data structures. With $\lambda(p,s)$05 preprocessing time, one can build a data structure of size $\lambda(p,s)$06 that returns an $\lambda(p,s)$07-approximation in $\lambda(p,s)$08 query time, where $\lambda(p,s)$09 is the doubling dimension (Li, 2010). A weaker “encoding scheme” computes an encoding $\lambda(p,s)$10 in $\lambda(p,s)$11 time so that from $\lambda(p,s)$12 and $\lambda(p,s)$13, EMD can be $\lambda(p,s)$14-approximated in $\lambda(p,s)$15 time (Li, 2010).

Several indexing methods are lower-bound based rather than approximation-scheme based. A Hilbert-curve method maps bins to one dimension, computes the 1D cumulative-difference EMD, and repeats with multiple random shifts; overall complexity is $\lambda(p,s)$16, and the paper states a theorem giving $\lambda(p,s)$17 approximation with probability at least $\lambda(p,s)$18 when $\lambda(p,s)$19 (Jang et al., 2011). For uncertain databases, another approach projects distributions onto a vector, approximates each projected distribution by a normal distribution plus error terms, represents each normal as a point in a Hough transformed space, and uses stochastic dominance for indexing. The resulting lower bound $\lambda(p,s)$20 satisfies

$\lambda(p,s)$21

and empirical speed-ups are reported for $\lambda(p,s)$22-NN query processing (Ruttenberg et al., 2011).

5. Learned and differentiable approximations

Approximate EMD is also used as a differentiable or fast surrogate in machine learning, especially when exact optimal transport is too slow to serve as a training loss.

DeepEMD approximates EMD between point clouds by predicting the matching itself rather than regressing the scalar distance. The input is the concatenation of two point clouds with a learned group-ID embedding and a small positional embedding, processed by a Transformer encoder with $\lambda(p,s)$23 multi-head self-attention layers, $\lambda(p,s)$24 heads, and model dimension $\lambda(p,s)$25. A matching head produces an attention matrix, from which directional matchings are extracted; the training loss is the average of two directional cross-entropies against the Hungarian one-hot matching (Sinha et al., 2023). The approximate EMD is then computed from the predicted assignments, and the paper reports that on point clouds of size $\lambda(p,s)$26, DeepEMD is approximately $\lambda(p,s)$27 faster than Hungarian and approximately $\lambda(p,s)$28 faster than $\lambda(p,s)$29-iteration Sinkhorn, with Pearson $\lambda(p,s)$30, Spearman $\lambda(p,s)$31, and Kendall $\lambda(p,s)$32 on held-out pairs (Sinha et al., 2023).

A CNN-based differentiable approximation is used for detector-image comparisons in high-luminosity LHC data compression. The surrogate network takes a pair of $\lambda(p,s)$33 tensors, applies four 2D-convolutional layers with $\lambda(p,s)$34 filters and kernel size $\lambda(p,s)$35, one dense layer with $\lambda(p,s)$36 neurons, and a final linear output, with symmetrization

$\lambda(p,s)$37

Ground-truth labels are computed offline via ot.emd2(...), and the surrogate is trained with mean-squared error on EMD (Shenoy et al., 2023). On the validation set, the reported mean relative error is approximately $\lambda(p,s)$38 and $\lambda(p,s)$39, while a single forward pass on a GPU takes $\lambda(p,s)$40 ms versus $\lambda(p,s)$41 ms for a standard linear-program EMD solver, giving an approximately $\lambda(p,s)$42 speed-up in inference (Shenoy et al., 2023).

These methods do not provide the same type of guarantee as a provable $\lambda(p,s)$43-approximation algorithm. DeepEMD emphasizes accurate estimation of the EMD and its gradient, explicit matching prediction, and strong out-of-distribution behavior on point clouds several times larger than those seen during training (Sinha et al., 2023). The LHC surrogate emphasizes differentiability and downstream utility as a training loss for an autoencoder-inspired neural network, where training with the EMD CNN loss improves median EMD and cluster $\lambda(p,s)$44 resolution relative to MSE-based losses (Shenoy et al., 2023). A plausible implication is that learned approximations are most useful when gradient fidelity or end-to-end optimization is more important than worst-case approximation guarantees.

6. Guarantees, lower bounds, and limitations

Approximate EMD methods are constrained not only by implementation trade-offs but also by negative results in fine-grained complexity. Under the Orthogonal Vectors Conjecture, there is a constant $\lambda(p,s)$45 such that exact EMD in dimension $\lambda(p,s)$46 cannot be computed in truly subquadratic time (Rohatgi, 2019). Under the Hitting Set Conjecture, for every constant $\lambda(p,s)$47 and any $\lambda(p,s)$48, there is no truly subquadratic-time algorithm that outputs a bijection with cost at most $\lambda(p,s)$49 (Rohatgi, 2019). The paper also shows that asymmetric EMD matching cannot be approximated in truly subquadratic time for $\lambda(p,s)$50 in $\lambda(p,s)$51 dimensions (Rohatgi, 2019).

These hardness statements matter because many nearly-linear-time or subquadratic approximation schemes exhibit exponential dependence on dimension or weaken the problem in some other way. The conditional lower bounds explicitly note that all known nearly-linear-time approximation schemes for EMD have a factor $\lambda(p,s)$52 or $\lambda(p,s)$53 dependence, and that improving the dependence on $\lambda(p,s)$54 or beating quadratic in $\lambda(p,s)$55 for nontrivial approximation would refute standard conjectures in fine-grained complexity (Rohatgi, 2019).

Limitations also arise within individual approximation paradigms. In the continuous-object setting, when neither side is a point set, purely multiplicative approximation is impossible because the optimal transport cost can be arbitrarily small (Kreveld et al., 2021). In ACT-$\lambda(p,s)$56 and related lower-bound relaxations, approximation quality depends on histogram overlap pattern, and there are no strict additive or multiplicative error bounds in high dimensions (Atasu et al., 2018). In NNS-EMD, the method is presented as achieving high accuracy, low time complexity, and high memory efficiency, but it remains an iterative approximation whose practical behavior depends on the number of active-set reductions and on the choice between the Random Protocol and Greedy Protocol (Meng et al., 2024). In learned surrogates, DeepEMD still has quadratic attention cost, and hard argmax assignments break gradient flow; the LHC surrogate gives no formal proof of convergence to the true EMD (Sinha et al., 2023, Shenoy et al., 2023).

Taken together, these results delineate the current meaning of approximate EMD. Provable $\lambda(p,s)$57 schemes exist for several geometric and discrete settings, additive approximation breaks the quadratic barrier for arbitrary costs, sketching and LSH support search and indexing, and differentiable surrogates enable gradient-based optimization. At the same time, dimensional dependence, additive terms, support assumptions, and hardness results impose clear limits on what approximation can achieve in general (Kreveld et al., 2021, Beretta et al., 2023, Rohatgi, 2019).