The Sinkhorn divergence is not a metric, yet it now has geodesic equations: its geometric estimators can now report uncertainty

The space of probability measures now carries a Levi-Civita covariant derivative for the Sinkhorn divergence, and geodesics, exponential map, logarithmic map and parallel transport all follow from it. Wasserstein geometry already had the full second-order calculus running from covariant differentiation to parallel transport, while the Sinkhorn side had only a tangent space and a metric tensor, with the geometric objects underived. The contrast sits in the projection inside the covariant derivative: Wasserstein requires solving a weighted elliptic partial differential equation, whereas the Sinkhorn projection is a single operator on a fixed RKHS, with both projected objects living in that same space.

The geodesic and parallel transport equations each become an ordinary differential equation with smooth coefficients on a fixed reproducing kernel Hilbert space (RKHS, the function space generated by a kernel). Reaching that form requires first proving the metric tensor is Hadamard differentiable along curves of measures, meaning its directional derivative exists uniformly over all sequences converging to the point, before the covariant derivative can be assembled along Gigli's route. The route is restricted to costs of finite smoothness: the squared Euclidean cost does not qualify, because the Gaussian kernel RKHS it induces is not closed under pointwise multiplication, while the Matérn kernel does qualify and approximates the squared Euclidean cost arbitrarily well.

Linearizing the solutions of those ODEs in the initial data gives √n Gaussian limits for plug-in estimators of the exponential map, the logarithmic map, parallel transport and the squared Sinkhorn distance, converging to a centered Gaussian process around the population value. Before this, in a transport-based geometry these canonical geometric operations had neither √n limit theorems nor bootstrap consistency results. Under nonconjugacy and a unique constant-speed variational minimizer, the whole estimated parallel-imputation curve (the initial direction of one geodesic carried to a new base point) is Gaussian in the limit as well, and the nonparametric bootstrap reproduces that limit consistently in outer probability, so no separate asymptotic variance estimator is needed.

The results cover costs of finite smoothness whose induced RKHS is norm-equivalent to a Sobolev space, population measures with a Lebesgue density uniformly bounded below, a logarithmic map used only near nonconjugate pairs, and geodesics taken as the unique constant-speed variational minimizer. The discussion proposes a next step of replacing these infinite-dimensional plug-in estimators with finite-dimensional approximations whose error is $o_P(n^{-1/2})$, so that the established limit distributions and bootstrap guarantees are inherited; another direction is examining whether the density lower bound can be relaxed.

The strongest of the three assumptions is the uniform positive lower bound on the density. Solutions of the ambient geodesic equations are only guaranteed to stay on the unit sphere with unit mass, not to remain probability measures, and they coincide with the probability-valued interpolant only when the initial velocity comes from a variational minimizer satisfying that bound; the parallel-imputation curve is not guaranteed to be a probability measure even at the population level. The numerical example only shows that the minimum density along the discretized interpolant stays positive when the endpoint densities have a positive lower bound; whether the endpoint condition implies a bound along the interpolant itself is not established.