Posterior sampling for inverse imaging problems on the sphere in seismology and cosmology

Abstract

In this work, we describe a framework for solving spherical inverse imaging problems using posterior sampling for full uncertainty quantification. Inverse imaging problems defined on the sphere arise in many fields, including seismology and cosmology where images are defined on the globe and the cosmic sphere, and are generally high-dimensional and computationally expensive. As a result, sampling the posterior distribution of spherical imaging problems is a challenging task. Our framework leverages a proximal Markov chain Monte Carlo (MCMC) algorithm to efficiently sample the high-dimensional space of spherical images with a sparsity-promoting wavelet prior. We detail the modifications needed for the algorithm to be applied to spherical problems, and give special consideration to the crucial forward modelling step which contains computationally expensive spherical harmonic transforms. By sampling the posterior, our framework allows for full and flexible uncertainty quantification, something which is not possible with other methods based on, for example, convex optimisation. We demonstrate our framework in practice on full-sky cosmological mass-mapping and to the construction of phase velocity maps in global seismic tomography. We find that our approach is potentially useful at moderate resolutions, such as those of interest in seismology. However at high resolutions, such as those required for astrophysical applications, the poor scaling of the complexity of spherical harmonic transforms severely limits our method, which may be resolved with future GPU implementations. A new Python package, pxmcmc, containing the proximal MCMC sampler, measurement operators, wavelet transforms and sparse priors is made publicly available.

PxMCMC: A Python package for proximal Markov Chain Monte Carlo

“…The use of proximal operators in MCMC was first proposed by Pereyra (2016), modifying the gradient-based Langevin MCMC, and has since been used in astrophysical and geophysical applications (e.g. Cai et al, 2018;Marignier et al, 2023;Price et al, 2020).…”

“…The class-based API abstracts out the main components of MCMC into interoperable classes, thereby allowing users to implement their own forward models (physical model) and priors, and even their own MCMC sampler if desired. Originally developed to solve inverse imaging problems defined on spherical domains Marignier et al, 2023), the package provides priors to promote sparsity in a spherical wavelet domain using transforms from the S2LET package (Leistedt et al, 2013). Examples provided in the package include a common problem in global seismic tomography and a full-sky cosmological mass-mapping problem, the details of which can be found in Marignier et al (2023) and .…”

“…Originally developed to solve inverse imaging problems defined on spherical domains Marignier et al, 2023), the package provides priors to promote sparsity in a spherical wavelet domain using transforms from the S2LET package (Leistedt et al, 2013). Examples provided in the package include a common problem in global seismic tomography and a full-sky cosmological mass-mapping problem, the details of which can be found in Marignier et al (2023) and .…”