Numerical methods for the elliptic Monge-Ampere equation and optimal transport

The problem of optimal transport, which involves finding the most cost-efficient way of transporting mass from one location to another, has been widely-studied, going back to the late eighteenth century. Recent years have revealed numerous applications in areas such as medical imaging, meteorology, cosmology, oceanography, and economics. Despite the importance of optimal transport, the computation of solutions remains extremely challenging. In the simplest case, where the cost function is quadratic, the problem takes on additional structure. In this setting, the constraint that mass must be conserved can be expressed as a fully non-linear partial differential equation known as the elliptic Monge-Ampere equation. The numerical solution of the Monge-Ampere equation has received a great deal of attention in recent years, yet the correct and efficient computation of solutions remains a challenge. Because of the nonlinearity of the equation, solutions can be singular and standard numerical approaches can fail. This means that novel solution techniques are needed to correctly capture the behaviour of weak solutions. We describe a monotone finite difference discretisation, which provably converges to the viscosity solution of the Monge-Ampere equation. The accuracy of the discretisation is improved by combining higher-order schemes with the monotone scheme needed to capture the correct behaviour of solutions near singularities. In doing this, we provide a general result about the convergence of higher-order finite difference methods for elliptic equations. The resulting nonlinear equations are solved efficiently using Newton's method. To ensure that mass is mapped into the desired region, the Monge-Ampere equation must be coupled to a transport boundary condition. This type of boundary condition is non-standard, and previously has been implemented only in very simple cases (such as transporting a square to a square). We propose a new method for implementing the transport condition by solving a sequence of more tractable Monge-Ampere equations with Neumann boundary conditions. To demonstrate the effectiveness and efficiency of the resulting methods, we provide computational results for a number of challenging problems including the recovery of inverse maps, mapping onto unbounded density functions, mapping from a disconnected domain, and mapping onto non-convex sets.