Why "Moving Dirt" Is One of the Hardest Problems in Math

Here is a problem that sounds almost too simple to be interesting.

You have a pile of dirt. Somewhere nearby, there is a hole. The pile and the hole have the same total volume, but different shapes. Your task: move the dirt from the pile into the hole, using as little total effort as possible.

Not just “shovel it over there.” Optimally. Every grain of dirt should travel the minimum distance necessary to fill the hole, given all the other grains are doing the same. If you move one grain too far, another grain moves too little, and the total effort increases.

This problem — now called the Optimal Transport problem, or the Monge problem after the French mathematician Gaspard Monge who first posed it in 1781 — turns out to be extraordinarily difficult. It took over 150 years to solve rigorously. Fields Medalists have worked on it. It connects measure theory, probability, geometry, and computation in ways that are still being explored today.

It also, as I will explain, has a great deal to do with imaging the inside of the Earth.

What Makes It Hard

The difficulty is not immediately obvious. It seems like the kind of problem a clever algorithm should dispatch in seconds.

The trouble is scale and coupling. Every grain of dirt is not independent — where you send one grain affects the optimal destination for every other grain. In the continuous mathematical setting, you are not moving individual particles but redistributing an entire probability distribution into another shape, while minimizing a global cost function that depends on all movement decisions simultaneously.

Monge himself could not solve the general problem. It was Leonid Kantorovich — a Soviet mathematician who would later win the Nobel Prize in Economics — who found the path to a rigorous solution in the 1940s. His key insight was subtle but transformative: he allowed mass to split. In Monge’s original formulation, each grain of dirt must go to exactly one destination. Kantorovich relaxed this constraint, permitting a single grain to be distributed across multiple destinations — converting the problem into a linear program that could be solved rigorously. The result is the Monge-Kantorovich problem, which underlies an entire field: optimal transport theory.

From Monge's problem to Kantorovich's reformulation — From Monge’s intuition to Kantorovich’s reformulation: the critical step was allowing mass to split across destinations, transforming an intractable combinatorial problem into a linear program — and ultimately yielding the Wasserstein distance as a principled measure of distributional difference.

The mathematical object at the center of this field is the Wasserstein distance — a way of measuring how different two probability distributions are, not by comparing them point-by-point, but by computing the minimum cost of transforming one into the other.

This distinction matters enormously.

Two Ways to Compare Two Things

Suppose you want to measure the difference between two shapes — say, two seismic waveforms recorded at the same location, one from a real earthquake and one from a computer simulation.

The standard approach, used throughout science and engineering for decades, is the L2 norm: square the difference at each point, sum everything up, take the square root. It is simple, mathematically convenient, and computationally cheap.

It is also, in certain situations, deeply misleading.

Here is the problem. Imagine two waveforms that are identical in shape but offset in time — one is simply a shifted version of the other. Visually, they look almost the same. A human looking at both would immediately say: these are nearly identical, just slightly misaligned.

The L2 norm disagrees. Because it compares corresponding time points directly, a shifted waveform registers as very different from the original, even when the underlying physics is nearly identical. It has no concept of “these two signals describe the same event at slightly different times.” It only knows: at this exact moment, the values differ by this amount. The consequence, visible in the misfit landscape, is a jagged surface full of local minima — many wrong answers that look, locally, like right answers.

The Wasserstein distance handles this correctly. Because it measures the cost of transforming one distribution into another, it naturally accounts for shifts, stretches, and warps. Two similar-shaped distributions that are slightly offset will be close in Wasserstein distance. The misfit landscape becomes smoother and more convex — meaning the optimization algorithm has fewer traps to fall into.

The analogy I find useful: L2 is like comparing two pieces of music by looking at their sheet music note-by-note and counting mismatches. Wasserstein distance is like comparing the melodic contour — the overall shape of the music as it moves through time. The second approach is more like how humans actually perceive similarity.

L2 norm versus Wasserstein distance waveform comparison — Two ways to measure waveform mismatch. L2 norm compares point-by-point: a slight time shift registers as a large error, and the misfit landscape is jagged with many local minima. Wasserstein distance measures the cost of transforming one signal into the other: shifted signals remain close, and the misfit landscape is smoother and more convex — a far more reliable guide for optimization.

The Earthquake Imaging Problem

Now for the connection to geophysics.

Full waveform inversion — FWI, in the shorthand of the field — is a technique for imaging the interior of the Earth using seismic waves. The principle is elegant: send seismic waves into the ground (via an earthquake, an explosion, or a controlled source), record how those waves return to the surface, and then work backwards to determine what the underground structure must look like in order to produce exactly those recorded signals.

This inverse problem is one of the most computationally demanding in applied science. The Earth’s interior is three-dimensional and extraordinarily complex. The relationship between subsurface structure and surface recordings is non-linear and non-unique. And the datasets involved — continuous recordings from networks of seismometers — are vast.

The standard approach to FWI uses the L2 norm to measure the mismatch between observed and simulated waveforms, then iteratively adjusts the subsurface model to reduce that mismatch. At each step, the algorithm asks: in which direction should I adjust the model to make the simulation more closely match the data?

This works, under favorable conditions. The conditions are: the initial model is already close to the truth, and the waveforms are well-aligned in time.

When those conditions fail — which happens often in practice — the L2-based approach falls into a trap called cycle skipping.

The Trap of Cycle Skipping

A seismic wave is periodic. It oscillates. And if the simulated wave is more than half a wavelength out of phase with the observed wave, something goes wrong.

The gradient of the L2 misfit — the signal that tells the inversion algorithm which direction to move — points in the wrong direction. The algorithm “thinks” it is getting closer to the correct model when it is actually moving away. The misfit value decreases, iteration after iteration, but the recovered model diverges from the truth. It converges confidently to a wrong answer.

Think of navigating in complete darkness with only a compass that tells you whether you are getting closer to or farther from your destination. If the compass is calibrated against a reference point that is itself wrong — if your “target” signal is actually pointing at the mirror image of where you want to go — you will walk with perfect confidence in exactly the wrong direction. The L2 misfit landscape is full of such false basins: local minima where the algorithm settles, having “improved” its way into an incorrect model.

Cycle skipping trap infographic — The cycle skipping trap: when the simulated waveform is more than half a wavelength out of phase with the observed data, the L2 gradient points away from the true model. The algorithm appears to converge — the misfit decreases — but the recovered model is wrong. The misfit landscape reveals why: many deceptive local minima surround the true global minimum, and L2-based inversion can descend into the wrong basin from the very first iteration.

Cycle skipping is the reason FWI often requires an excellent initial model. Without one, the L2 gradient is unreliable, and the inversion diverges. This is a significant practical limitation: in many real-world scenarios, such as imaging beneath complex geological structures or in regions without prior well data, a sufficiently accurate starting model simply does not exist.

Optimal Transport as a Better Misfit

This is where optimal transport enters the picture — not as a theoretical curiosity, but as a practical solution to a specific failure mode.

If the mismatch between observed and simulated waveforms is measured with Wasserstein distance rather than L2, the misfit landscape changes fundamentally. The Wasserstein misfit is more convex — meaning it has fewer local minima for the inversion algorithm to get trapped in. It is more sensitive to the overall shape and timing of waveforms than to point-by-point differences. And critically, it degrades more gracefully when the initial model is poor: even when the simulated and observed waveforms are significantly misaligned, the Wasserstein gradient still points, roughly, in the right direction.

My PhD research focused on making this idea computationally tractable. The theoretical appeal of Wasserstein distance in FWI had been recognized before my work — the challenge was that computing exact Wasserstein distances is expensive, and exact optimal transport between high-dimensional distributions is even more so.

The approach I developed and studied is based on the Sinkhorn divergence — a regularized approximation to the Wasserstein distance that introduces an entropic penalty term. This regularization does two things: it makes the problem strictly convex (which is mathematically convenient), and it allows the distance to be computed using an efficient iterative algorithm — the Sinkhorn-Knopp algorithm — that scales to large problems. The work involved systematic comparison of four OT-based misfit functions — W1, W2, GHK, and Sinkhorn divergence — across synthetic and realistic inversion scenarios, studying how each behaved as a function of initial model quality, noise levels, and data complexity.

Sinkhorn divergence infographic — Sinkhorn divergence makes optimal transport practical for large-scale FWI. By introducing entropic regularization, the Sinkhorn-Knopp algorithm computes approximate Wasserstein distances efficiently — yielding a misfit landscape that is smoother, more convex, and more robust to misalignment than the standard L2 norm. The practical applications span from oil and gas exploration to earthquake hazard assessment and urban subsurface mapping.

The conclusions were not simple: no single choice dominates all scenarios, and the interaction between the misfit function and the inversion algorithm matters as much as the misfit function itself.

This kind of nuance is, I think, typical of good applied mathematics. The elegant idea — use Wasserstein distance instead of L2 — is the beginning of the work, not the end of it.

What It Means in Practice

The practical implications extend beyond academic interest.

Seismic imaging underlies a remarkable range of applications: oil and gas exploration, earthquake hazard assessment, carbon sequestration monitoring, geothermal resource characterization, and — increasingly — urban subsurface mapping for infrastructure planning. In all of these contexts, the quality of the subsurface image depends directly on the quality of the inversion algorithm.

An inversion method that is more robust to poor initial models, more tolerant of data noise, and less prone to cycle skipping is not a minor refinement. It is the difference between a technique that works in favorable laboratory conditions and one that can be deployed in the complex, data-limited, high-stakes conditions of real-world geophysical surveys.

The research implementation I developed during my PhD is publicly available. Releasing the code openly was a deliberate choice: the barrier between theoretical work and practical adoption in geophysics is often the absence of a working, documented implementation that practitioners can test and build upon. Lowering that barrier is part of what it means to do science that actually reaches the world.

The Unreasonable Usefulness of Abstract Mathematics

In 1781, Gaspard Monge was thinking about logistics — specifically, about the optimal way to move excavated earth during military fortification construction. He was not thinking about seismic waves, or probability distributions, or machine learning, all of which would use his ideas extensively two centuries later.

This is a familiar pattern in mathematics. The abstract structures that mathematicians develop for one purpose — often for purely intellectual reasons, without any application in mind — turn out to be exactly the right language for problems that haven’t been invented yet.

Optimal transport is a particularly striking example. Monge’s problem, Kantorovich’s reformulation, Wasserstein’s distance, and the computational breakthroughs of the 1990s and 2000s have collectively created a mathematical framework that now appears in image processing, generative AI, economics, biology, climate science, and — yes — earthquake imaging. At the end of that chain sits something concrete: a color-coded map of seismic velocity in the Earth’s interior, layers of rock and sediment rendered visible by mathematics that began with a pile of dirt and a hole in the ground.

The ground beneath a city holds a record of everything that has happened to it: every geological event, every structural modification, every load cycle. Reading that record requires mathematics sophisticated enough to extract signal from noise, to compare complex shapes, to navigate high-dimensional optimization landscapes without getting lost.

It turns out that the best way to move dirt — really, the mathematically optimal way — is also, in a precise sense, one of the best ways to image what lies beneath it.

Monge didn’t know it in 1781. But he was also thinking about earthquakes.

Miao YU is a geophysicist whose PhD research at Université de Paris focused on entropic optimal transport methods for seismic full waveform inversion. She subsequently led R&D in AI-driven urban infrastructure monitoring systems, working at the intersection of geophysical theory and large-scale engineering deployment.