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

TL;DR

This paper introduces a proximal MCMC framework for efficiently sampling high-dimensional spherical inverse problems, enabling full uncertainty quantification in seismology and cosmology applications.

Contribution

It develops a modified proximal MCMC algorithm tailored for spherical problems with wavelet priors, including computational strategies for spherical harmonic transforms.

Findings

Effective in moderate-resolution cosmological mass-mapping

Applicable to global seismic tomography problems

Provides full uncertainty quantification

Abstract

Inverse problems defined on the sphere arise in many fields, including seismology and cosmology where problems are defined on the globe and the cosmic sphere. These are generally high-dimensional and computationally very complex and, as a result, sampling the posterior of spherical inverse problems is a challenging task. In this work, we describe a framework that leverages a proximal Markov chain Monte Carlo (MCMC) algorithm to efficiently sample the high-dimensional space of spherical inverse problems 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 spherical harmonic transforms that are computationally expensive. 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 on a common problem in global seismic tomography. We find that our approach is potentially useful at moderate resolutions, such as those of interest in seismology. Our framework is generally limited by resolution requirements, such as those required for astrophysical applications, due to the poor scaling of the complexity of spherical harmonic transforms with resolution. A new Python package, pxmcmc, containing the proximal MCMC sampler, measurement operators, wavelet transforms and sparse priors is made publicly available.