Brittany D. Froese

Yunan Yang, Björn Engquist, Junzhe Sun et al.

Conventional full-waveform inversion (FWI) using the least-squares norm ($L^2$) as a misfit function is known to suffer from cycle skipping. This increases the risk of computing a local rather than the global minimum of the misfit. In our previous work, we proposed the quadratic Wasserstein metric ($W_2$) as a new misfit function for FWI. The $W_2$ metric has been proved to have many ideal properties with regards to convexity and insensitivity to noise. When the observed and predicted seismic data are regarded as two density functions, the quadratic Wasserstein metric corresponds to the optimal cost of rearranging one density into the other, where the transportation cost is quadratic in distance. The difficulty of transforming seismic signals into nonnegative density functions is discussed. Unlike the $L^2$ norm, $W_2$ measures not only amplitude differences, but also global phase shifts, which helps to avoid cycle skipping issues. In this work, we build on our earlier method to cover more realistic high-resolution applications by embedding the $W_2$ technique into the framework of the adjoint-state method and applying it to seismic relevant 2D examples: the Camembert, the Marmousi, and the 2004 BP models. We propose a new way of using the $W_2$ metric trace-by-trace in FWI and compare it to global $W_2$ via the solution of the Monge-Ampère equation. With corresponding adjoint source, the velocity model can be updated using the l-BFGS method. Numerical results show the effectiveness of $W_2$ for alleviating cycle skipping issues and sensitivity to noise. Both mathematical theory and numerical examples demonstrate that the quadratic Wasserstein metric is a good candidate for a misfit function in seismic inversion.

Jean-David Benamou, Brittany D. Froese, Adam M. Oberman

A numerical method for the solution of the elliptic Monge-Ampere Partial Differential Equation, with boundary conditions corresponding to the Optimal Transportation (OT) problem is presented. A local representation of the OT boundary conditions is combined with a finite difference scheme for the Monge-Ampere equation. Newton's method is implemented leading to a fast solver, comparable to solving the Laplace equation on the same grid several times. Theoretical justification for the method is given by a convergence proof in the companion paper (Benamou et al., 2012). In this paper, the algorithm is modified to a simpler compact stencil implementation and details of the implementation are given. Solutions are computed with densities supported on non-convex and disconnected domains. Computational examples demonstrate robust performance on singular solutions and fast computational times.

Bjorn Engquist, Brittany D. Froese, Yunan Yang

Full waveform inversion is a successful procedure for determining properties of the earth from surface measurements in seismology. This inverse problem is solved by a PDE constrained optimization where unknown coefficients in a computed wavefield are adjusted to minimize the mismatch with the measured data. We propose using the Wasserstein metric, which is related to optimal transport, for measuring this mismatch. Several advantageous properties are proved with regards to convexity of the objective function and robustness with respect to noise. The Wasserstein metric is computed by solving a Monge-Ampere equation. We describe an algorithm for computing its Frechet gradient for use in the optimization. Numerical examples are given.

Brittany D. Froese, Adam M. Oberman

The elliptic Monge-Ampère equation is a fully nonlinear Partial Differential Equation that originated in geometric surface theory and has been applied in dynamic meteorology, elasticity, geometric optics, image processing and image registration. Solutions can be singular, in which case standard numerical approaches fail. Novel solution methods are required for stability and convergence to the weak (viscosity) solution. In this article we build a wide stencil finite difference discretization for the \MA equation. The scheme is monotone, so the Barles-Souganidis theory allows us to prove that the solution of the scheme converges to the unique viscosity solution of the equation. Solutions of the scheme are found using a damped Newton's method. We prove convergence of Newton's method and provide a systematic method to determine a starting point for the Newton iteration. Computational results are presented in two and three dimensions, which demonstrates the speed and accuracy of the method on a number of exact solutions, which range in regularity from smooth to non-differentiable.

Brittany D. Froese, Adam M. Oberman

The theory of viscosity solutions has been effective for representing and approximating weak solutions to fully nonlinear Partial Differential Equations (PDEs) such as the elliptic Monge-Ampère equation. The approximation theory of Barles-Souganidis [Barles and Souganidis, Asymptotic Anal., 4 (1999) 271-283] requires that numerical schemes be monotone (or elliptic in the sense of [Oberman, SIAM J. Numer. Anal, 44 (2006) 879-895]. But such schemes have limited accuracy. In this article, we establish a convergence result for nearly monotone schemes. This allows us to construct finite difference discretizations of arbitrarily high-order. We demonstrate that the higher accuracy is achieved when solutions are sufficiently smooth. In addition, the filtered scheme provides a natural detection principle for singularities. We employ this framework to construct a formally second-order scheme for the Monge-Ampère equation and present computational results on smooth and singular solutions.

Jean-David Benamou, Brittany D. Froese, Adam M. Oberman

In this work we present a numerical method for the Optimal Mass Transportation problem. Optimal Mass Transportation (OT) is an active research field in mathematics.It has recently led to significant theoretical results as well as applications in diverse areas. Numerical solution techniques for the OT problem remain underdeveloped. The solution is obtained by solving the second boundary value problem for the MA equation, a fully nonlinear elliptic partial differential equation (PDE). Instead of standard boundary conditions the problem has global state constraints. These are reformulated as a tractable local PDE. We give a proof of convergence of the numerical method, using the theory of viscosity solutions. Details of the implementation and a fast solution method are provided in the companion paper arXiv:1208.4870.

Brittany D. Froese, Adam M. Oberman, Tiago Salvador

The elliptic 2-Hessian equation is a fully nonlinear partial differential equation (PDE) that is related to intrinsic curvature for three dimensional manifolds. We introduce two numerical methods for this PDE: the first is provably convergent to the viscosity solution, and the second is more accurate, and convergent in practice but lacks a proof. The PDE is elliptic on a restricted set of functions: a convexity type constraint is needed for the ellipticity of the PDE operator. Solutions with both discretizations are obtained using Newton's method. Computational results are presented on a number of exact solutions which range in regularity from smooth to nondifferentiable and in shape from convex to non convex.

Brittany D. Froese

We introduce meshfree finite difference methods for approximating nonlinear elliptic operators that depend on second directional derivatives or the eigenvalues of the Hessian. Approximations are defined on unstructured point clouds, which allows for very complicated domains and a non-uniform distribution of discretisation points. The schemes are monotone, which ensures that they converge to the viscosity solution of the underlying PDE as long as the equation has a comparison principle. Numerical experiments demonstrate convergence for a variety of equations including problems posed on random point clouds, complex domains, degenerate equations, and singular solutions.

Jean-David Benamou, Brittany D. Froese

We consider the Monge-Kantorovich optimal transportation problem between two measures, one of which is a weighted sum of Diracs. This problem is traditionally solved using expensive geometric methods. It can also be reformulated as an elliptic partial differential equation known as the Monge-Ampere equation. However, existing numerical methods for this non-linear PDE require the measures to have finite density. We introduce a new formulation that couples the viscosity and Aleksandrov solution definitions and show that it is equivalent to the original problem. Moreover, we describe a local reformulation of the subgradient measure at the Diracs, which makes use of one-sided directional derivatives. This leads to a consistent, monotone discretisation of the equation. Computational results demonstrate the correctness of this scheme when methods designed for conventional viscosity solutions fail.

Brittany D. Froese

We consider the numerical approximation of surfaces of prescribed Gaussian curvature via the solution of a fully nonlinear partial differential equation of Monge-Ampère type. These surfaces need not be continuous up to the boundary of the domain and the Dirichlet boundary condition must be interpreted in a weak sense. As a consequence, sub-solutions do not always lie below super-solutions, standard comparison principles fail, and existing convergence theorems break down. By relying on a geometric interpretation of weak solutions, we prove a relaxed comparison principle that applies only in the interior of the domain. We provide a general framework for proving existence and stability results for consistent, monotone finite difference approximations and modify the Barles-Souganidis convergence framework to show convergence in the interior of the domain. We describe a convergent scheme for the prescribed Gaussian curvature equation and present several challenging examples to validate these results.

Brittany D. Froese, Tiago Salvador

We introduce generalised finite difference methods for solving fully nonlinear elliptic partial differential equations. Methods are based on piecewise Cartesian meshes augmented by additional points along the boundary. This allows for adaptive meshes and complicated geometries, while still ensuring consistency, monotonicity, and convergence. We describe an algorithm for efficiently computing the non-traditional finite difference stencils. We also present a strategy for computing formally higher-order convergent methods. Computational examples demonstrate the efficiency, accuracy, and flexibility of the methods.

Bjorn Engquist, Brittany D. Froese, Yen-Hsi Richard Tsai

The idea of using fast sweeping methods for solving stationary systems of conservation laws has previously been proposed for efficiently computing solutions with sharp shocks. We further develop these methods to allow for a more challenging class of problems including problems with sonic points, shocks originating in the interior of the domain, rarefaction waves, and two-dimensional systems. We show that fast sweeping methods can produce higher-order accuracy. Computational results validate the claims of accuracy, sharp shock curves, and optimal computational efficiency.