Regularity of optimal transport on oblate ellipsoids

in collaboration with Q. Mérigot



The question of regularity of optimal transport plans has received a lot of attention. Examples of discontinuous transport plans have been found by Caffarelli, between probability densities on the plane with non-convex support, and by Loeper on compact manifolds with negative curvature even for densities lower bounded by a positive constant (see e.g. Chapter 12 of the book of C. Villani). An even more striking recent result is that even on simple convex surfaces, it is possible to find examples of pairs of non-vanishing probability densities such that the optimal transport map is discontinuous (see recent papers of [B. Bonnard,J.B. Caillau,L. Rifford] and [A. Figalli,L. Rifford,C. Villani]). We study below numerically by our optimization method the case of oblate ellipsoids, which is one of the simplest cases where discontinuity can occur. More details can be found in our paper Discrete optimal transport: complexity, geometry and applications.

The family of ellipsoids $ S_\mu $ we will consider is defined by $ S_\mu := \left\{ (x,y,z)\in R^3; ~ x^2+y^2 + \left(\frac{z}{\mu}\right)^2 = 1 \right\} $ Bonnard, Caillau and Rifford have shown that for $ \mu $ smaller than $1/\sqrt{3}$, there exists a point on the equator $\{ z = 0 \}$ whose injectivity domain is non-convex. This fact implies the existence of two absolutely continuous probability measure $\mu,\nu$ with densities in $[C^{-1}, C]$ such that the optimal transport plan between $\mu$ and $\nu$ is discontinuous.

Quantization by 10000 Dirac masses of the initial and the target measures on the sphere and on an oblate ellipsoid (second row)

To obtain a qualitative description of this phenomena, we introduce two non-vanishing densities with equal total mass on $S_1$, whose formula in Cartesian coordinates is given by:

\[ f(x,y,z)= 0.05+4\exp\left(-c(1-x)^2\right) \]

and

\[ g(x,y,z)= 0.05+2\exp\left(-c(1-\frac{\sqrt{2}}{2}(z+y))^2\right) + 2\exp\left(-c(1-\frac{\sqrt{2}}{2}(z-y))^2\right) \]

with $c=500$. Previous parameters have been selected to define non-vanishing densities which are not too concentrated. Even in this regular situation, we will see that optimal transport maps may have very different qualitative behaviors depending on the eccentricity of the ellipsoid.

We construct a quantization by a sum of 10000 unit Dirac masses of the two measures of densities $f$ and $g$ with respect to the two dimensional Hausdorff measure on $S_1$ using Lloyd's algorithm. Our goal is to compare here the optimal transport between the two discrete measures on the sphere and between two analogous measures defined on the oblate ellipsoid $S_{\mu}$, with $\mu=0.3$. By analogous we mean that we deduce our two measures on $S_{\mu}$ by the naive bi-Lipschitz transformation

\[(x,y,z)\rightarrow (x,y,\mu z).\]

To illustrate the difference between the two optimal transports, we computed for every sampling point of the source densities the geodesic distance between its image and the center of the Gaussian part of the source density, i.e. the point $c:=(1,0,0)$. More formally, this is equivalent to compute the pullback of the geodesic distance function to $(1,0,0)$ by the optimal transport. These functions are rescaled to have the same infinity norm, and their values are reported in the following figures, with the same color map. One can observe that the function is smooth along the equator in the spherical case ( $\mu=1$) whereas in the oblate case ( $\mu=0.3$) some points close to the equator can be transported to points whose geodesic distance to $ c$ is very different. This qualitative behavior, observed here at a discrete level, is the first numerical illustration of the discontinuity of optimal transport plans between two non-vanishing densities supported on a convex manifold.

Pullback of the geodesic distance function to (1,0,0) by the optimal transports on the sphere and on the oblate ellipsoid respectively.

Move your mouse on the pictures to animate (slow connection be patient)

Color maps associated to the azimuthal coordinates of the image points in the case of the sphere and of an oblate ellipsoid