Diffeomorphic random sampling using optimal information transport

TL;DR

This paper introduces a diffeomorphic sampling algorithm on Riemannian manifolds based on optimal information transport, offering a semi-explicit alternative to optimal mass transport and MCMC methods.

Contribution

It develops a novel sampling method leveraging the Fisher-Rao metric and information geometry, avoiding complex nonlinear equations of traditional OMT.

Findings

Semi-explicit formulation of the sampling algorithm.

Potential advantages over MCMC for large sample sizes.

Connections between Fisher-Rao metric and diffeomorphism groups.

Abstract

In this article we explore an algorithm for diffeomorphic random sampling of nonuniform probability distributions on Riemannian manifolds. The algorithm is based on optimal information transport (OIT)---an analogue of optimal mass transport (OMT). Our framework uses the deep geometric connections between the Fisher-Rao metric on the space of probability densities and the right-invariant information metric on the group of diffeomorphisms. The resulting sampling algorithm is a promising alternative to OMT, in particular as our formulation is semi-explicit, free of the nonlinear Monge--Ampere equation. Compared to Markov Chain Monte Carlo methods, we expect our algorithm to stand up well when a large number of samples from a low dimensional nonuniform distribution is needed.

Figures

Equations

$${\phi }_{\ast }{\mu }_0=\mu .$$

$${\phi }_{\ast }{\mu }_0=|D\phi |{\mu }_0\circ \phi ,$$

$$\mathrm{Prob}(M)=\{\mu \in {\Omega }^d(M)\mid \mu >0,{\int }_M\mu =1\},$$

$$\mathrm{Diff}(M)\times \mathrm{Prob}(M)$$

$$(\phi ,\mu )$$

$$\mathrm{SDiff}(M)=\{\phi \in \mathrm{Diff}(M)\mid {\phi }_{\ast }{\mu }_0={\mu }_0\}.$$

$$E(\phi )={\mathrm{dist}}^2(\mathrm{id},\phi )$$

$$G_{\phi }(u\circ \phi ,v\circ \phi )={\int }_M⟨u\circ \phi ,v\circ \phi ⟩{\mu }_0$$

$$G_{\phi }^I(u\circ \phi ,v\circ \phi )=-{\int }_M⟨\Delta u,v⟩{\mu }_0+\lambda \sum _{i=1}^k{\int }_M⟨u,{\xi }_i⟩{\mu }_0{\int }_M⟨v,{\xi }_i⟩{\mu }_0,$$

$$\begin{array}{rl}\Delta f(t)& =\frac{\dot{\mu }(t)}{\mu (t)}\circ \phi (t),\\ v(t)& =\nabla (f(t)),\\ \frac{d}{dt}\phi (t{)}^{-1}& =v(t)\circ \phi (t{)}^{-1},\phi (0)=\mathrm{id}\end{array}$$

$$\mu (t)={\left(\frac{\sin \left((1-t)\theta \right)}{\sin \theta }+\frac{\sin \left(t\theta \right)}{\sin \theta }\sqrt{\frac{\mu }{{\mu }_0}}\right)}^2{\mu }_0,\cos \theta ={\int }_M\sqrt{\frac{\mu }{{\mu }_0}}{\mu }_0.$$

$$\Delta f_k=s_k.$$

$${\psi }_k=\mathrm{id}-\epsilon v_k.$$

$${\phi }_{k+1}={\phi }_k\circ {\psi }_k.$$

$$\mu \sim 3\exp (-x^2-10(y-x^2/2+1{)}^2)+2\exp (-(x+1{)}^2-y^2)+1/10,$$

Full text

Keywords: density matching, information geometry, Fisher–Rao metric, optimal transport, image registration, diffeomorphism groups, random sampling

MSC2010: 58E50, 49Q10, 58E10

1 Introduction

We construct algorithms for random sampling, addressing the following problem.

Problem 1

Let $\mu$ be a probability distribution on a manifold $M$. Generate $N$ random samples from $\mu$.

The classic approach to sample from a probability distribution on a higher dimensional space is to use Markov Chain Monte Carlo (MCMC) methods, for example the Metropolis–Hastings algorithm [6]. An alternative idea is to use diffeomorphic density matching between the density $\mu$ and a standard density $\mu_{0}$ from which samples can be drawn easily. Standard samples are then transformed by the diffeomorphism to generate non-uniform samples. In Bayesian inference, for example, the distribution $\mu$ would be the posterior distribution and $\mu_{0}$ would be the prior distribution. In case the prior itself is hard to sample from the uniform distribution can be used. For $M$ being a subset of the real line, the standard approach is to use the cumulative distribution function to define the diffeomorphic transformation. If, however, the dimension of $M$ is greater then one there is no obvious change of variables to transform the samples to the distribution of the prior. We are thus led to the following matching problem.

Problem 2

Given a probability distribution $\mu$ on $M$, find a diffeomorpism $\varphi$ such that

[TABLE]

Here, $\mu_{0}$ denotes a standard distribution on $M$ from which samples can be drawn, and $\varphi_{*}$ is the the push-forward of $\varphi$ acting on densities, i.e.,

[TABLE]

where $|D\varphi|$ is the Jacobian determinant.

A benefit of transport-based methods over traditional MCMC methods is cheap computation of additional samples; it amounts to drawing uniform samples and then evaluating the transformation. On the other hand, transport-based methods scale poorly with increasing dimensionality of $M$, contrary to MCMC.

The action of the diffeomorphism group on the space of smooth probability densities is transitive (Moser’s lemma [13]), so existence of a solution to Problem 2 is guaranteed. However, if the dimension of $M$ is greater then one, there is an infinite-dimensional space of solutions. Thus, one needs to select a specific diffeomorphism within the set of all solutions. Moselhy and Marzouk [12] and Reich [15] proposed to use optimal mass transport (OMT) to construct the desired diffeomorphism $\varphi$, thereby enforcing $\varphi=\nabla c$ for some convex function $c$. The OMT approach implies solving, in one form or another, the heavily non-linear Monge–Ampere equation for $c$. A survey of the OMT approach to random sampling is given by Marzouk et. al. [9].

In this article we pursue an alternative approach for diffeomorphic based random sampling, replacing OMT by optimal information transport (OIT), which is diffeomorphic transport based on the Fisher–Rao geometry [11]. Building on deep geometric connections between the Fisher–Rao metric on the space of probability densities and the right-invariant information metric on the group of diffeomorphisms [7, 11], we developed in [3] an efficient numerical method for density matching. The efficiency stems from a solution formula for $\varphi$ that is explicit up to inversion of the Laplace operator, thus avoiding the solution of nonlinear PDE such as Monge–Ampere. In this paper we explore this method for random sampling (the initial motivation in [3] is medical imaging, although other applications, including random sampling, are also suggested). The resulting algorithm is implemented in a short MATLAB code, available under MIT license at https://github.com/kmodin/oit-random.

2 Density Transport Problems

Let $M$ be an $d$–dimensional orientable, compact manifold equipped with a Riemannian metric $g=\langle.,.\rangle$. The volume density induced by $g$ is denoted $\mu_{0}$ and without loss of generality we assume that the total volume of $M$ with respect to $\mu_{0}$ is one, i.e., $\int_{M}\mu_{0}=1$. Furthermore, the space of smooth probability densities on $M$ is given by

[TABLE]

where $\Omega^{d}(M)$ denotes the space of smooth $d$-forms. The group of smooth diffeomorphisms $\operatorname{Diff}(M)$ acts on the space of probability densities via push-forward:

[TABLE]

By a result of Moser [13] this action is transitive.

We introduce the subgroup of volume preserving diffeomorphisms

[TABLE]

Note that ${\operatorname{SDiff}}(M)$ is the isotropy group of $\mu_{0}$ with respect to the action of $\mathrm{Diff}(M)$. The spaces $\operatorname{Prob}(M)$, $\mathrm{Diff}(M)$, and ${\operatorname{SDiff}}(M)$ all have the structure of smooth, infinite dimensional Fréchet manifold. Furthermore, $\mathrm{Diff}(M)$ and ${\operatorname{SDiff}}(M)$ are infinite dimensional Fréchet Lie groups. A careful treatment of these Fréchet topologies can be found in the work by Hamilton [5].

In the following we will focus our attention on the diffeomorphic density matching problem (Problem 2). A common approach to overcome the non-uniqueness in the solution is to add a regularization term to the problem. That is, to search for a minimum energy solution that has the required matching property, for some energy functional $E$ on the diffeomorphism group. Following ideas from mathematical shape analysis [10] it is a natural approach to define this energy functional using the geodesic distance function $\operatorname{dist}$ of a Riemannian metric on the diffeomorphism group. Then the regularized diffeomorphic matching problem can be written as follows.

Problem 3

Given a probability density $\mu\in\operatorname{Prob}(M)$ we want to find the diffeomorphism $\varphi\in\operatorname{Diff}(M)$ that minimizes the energy functional

[TABLE]

over all diffeomorphisms $\varphi$ with $\varphi_{*}\mu_{0}=\mu$.

The free variable in the above matching problem is the choice of Riemannian metric—thus distance function—on the group of diffeomorphisms. Although not formulated as here, Moselhy and Marzouk [12] proposed to use the $L^{2}$ metric on $\mathrm{Diff}(M)$

[TABLE]

for $u\circ\varphi,v\circ\varphi\in T_{\varphi}\mathrm{Diff}(M)$. This corresponds to distance-squared optimal mass transport (OMT), which induces the Wasserstein $L^{2}$ distance on $\mathrm{Prob}(M)$, see, for example, [14, 8, 16].

In this article we use the right-invariant $H^{1}$-type metric

[TABLE]

where $\lambda>0$, $\Delta$ is the Laplace–de Rham operator lifted to vector fields, and $\xi_{1},\ldots,\xi_{k}$ is an orthonormal basis of the harmonic 1-forms on $M$. Because of the Hodge decomposition theorem, $G^{I}$ is independent of the choice of orthonormal basis $\xi_{1},\ldots,\xi_{k}$ for the harmonic vector fields. This construction is related to the Fisher-Rao metric on the space of probability density [4, 2], which is predominant in the field of information geometry [1]. We call $G^{I}$ the information metric. See [7, 11, 3] for more information on the underlying geometry.

The connection between the information metric and the Fisher-Rao metric allows us to construct almost explicit solutions formulas for Problem 2 using the explicit formulas for the geodesics of the Fisher-Rao metric.

Theorem 1** ([11, 3])**

Let $\mu\in\mathrm{Prob}(M)$ be a smooth probability density. The diffeomorphism $\varphi\in\mathrm{Diff}(M)$ minimizing $\operatorname{dist}_{G^{I}}(\operatorname{id},\varphi)$ under the constraint $\varphi_{*}\mu_{0}=\mu$ is given by $\varphi(1)$, where $\varphi(t)$ is obtained as the solution to the problem

[TABLE]

where $\mu(t)$ is the (unique) Fisher-Rao geodesic connecting $\mu_{0}$ and $\mu$

[TABLE]

The algorithm for diffeomorphic random sampling, described in the following section, is directly based on solving the equations (8).

3 Numerical Algorithm

In this section we explain the algorithm for random sampling using optimal information transport. It is a direct adaptation of [3, Algorithm 1].

Algorithm 1 (OIT based random sampling)

  Assume we have a numerical way to represent functions, vector fields, and diffeomorphisms on $M$, and numerical methods for

• •

composing functions and vector fields with diffeomorphisms,

• •

computing the gradient of functions,

• •

computing solutions to Poisson’s equation on $M$,

• •

sampling from the standard distribution $\mu_{0}$ on $M$, and

• •

evaluating diffeomorphisms.

An OIT based algorithm for Problem 1 is then given as follows:

Choose a step size $\varepsilon=1/K$ for some positive integer $K$ and calculate the Fisher-Rao geodesic $\mu(t)$ and its derivative $\dot{\mu}(t)$ at all time points $t_{k}=\frac{k}{K}$ using equation (9). 2. 2.

Initialize $\varphi_{0}=\operatorname{id}$. Set $k\leftarrow 0$. 3. 3.

Compute $s_{k}=\frac{\dot{\mu}(t_{k})}{\mu(t_{k})}\circ\varphi_{k}$ and solve the Poisson equation

[TABLE] 4. 4.

Compute the gradient vector field $v_{k}=\nabla f_{k}$. 5. 5.

Construct approximations $\psi_{k}$ to $\exp(-\varepsilon v_{k})$, for example

[TABLE] 6. 6.

Update the diffeomorphism111If needed, one may also compute the inverse by $\varphi_{k+1}^{-1}=\varphi_{k}^{-1}+\varepsilon v\circ\varphi_{k}^{-1}$.

[TABLE] 7. 7.

Set $k\leftarrow k+1$ and continue from step 3 unless $k=K$. 8. 8.

Draw $N$ random samples $x_{1},\ldots x_{N}$ from the uniform distribution $\mu_{0}$. 9. 9.

Set $y_{n}=\varphi_{K}(x_{n})$, $n\in\{1,\ldots N\}$.

The algorithm generates $N$ random samples $y_{1},\ldots,y_{N}$ from the distribution $\mu$. One can save $\varphi_{K}$ and repeat 8-9 whenever additional samples are needed.

The computationally most intensive part of the algorithm is the solution of Poisson’s equation at each time step. Notice, however, that we do not need to solve nonlinear equations, such as Monge–Ampere, as is necessary in OMT.

4 Example

In this example we consider $M=\mathbb{T}^{2}\simeq(\mathbb{R}/2\pi\mathbb{Z})^{2}$ with distribution defined in Cartesian coordinates $x,y\in[-\pi,\pi)$ by

[TABLE]

normalized so that the ratio between the maximum and mimimum of $\mu$ is 100. The resulting density is depicted in Fig. 1 (left).

We draw $10^{5}$ samples from this distribution using a MATLAB implementation of our algorithm, available under MIT license at

https://github.com/kmodin/oit-random

The implementation can be summarized as follows. To solve the lifting equations (8) we discretize the torus by a $256\times 256$ mesh and use the fast Fourier transform (FFT) to invert the Laplacian. We use 100 time steps. The resulting diffeomorphism is shown as a mesh warp in Fig. 2. We then draw $10^{5}$ uniform samples on $[-\pi,\pi]^{2}$ and apply the diffeomorphism on each sample (applying the diffeomorphism corresponds to interpolation on the warped mesh). The resulting random samples are depicted in Fig. 1 (right). To draw new samples is very efficient. For example, another $10^{7}$ samples can be drawn in less than a second on a standard laptop.

5 Conclusions

In this paper we explore random sampling based on the optimal information transport algorithm developed in [3]. Given the semi-explicit nature of the algorithm, we expect it to be an efficient competitor to existing methods, especially for drawing a large number of samples from a low dimensional manifold. However, a detailed comparison with other methods, including MCMC methods, is outside the scope of this paper and left for future work.

We provide an example of a complicated distribution on the flat 2-torus. The method is straighforward to extended to more elaborate manifolds, e.g., by using finite element methods for the efficient solution of Poisson’s equation on manifolds. For non-compact manifolds, most importantly $\mathbb{R}^{n}$, one might use standard techniques, such as Box–Muller, to first transform the required distribution to a compact domain.