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

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

Augustin Marignier  1,2 Jason D. McEwen E-mail:
<EMAIL_ADDRESS>1 Ana M. G. Ferreira  2,3 Thomas D. Kitching1
1Mullard Space Science Laboratory  University College London  Dorking RH5 6NT
UK
2Department of Earth Sciences  University College London  London WC1E 6BT  UK
3CERIS  Instituto Superior Técnico  Universidade de Lisboa  Lisbon  Portugal

(Accepted XXX. Received YYY; in original form ZZZ; 2022)

###### 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.

###### keywords:

Data Methods — Bayesian — Mass-Mapping — Seismic Tomography

##  1 Introduction

Inverse problems on the sphere are common in many fields, from astrophysics
(e.g. Planck Collaboration, 2016; Jeffrey et al., 2021), to geophysics (e.g.
Ritsema et al., 2011; Chang et al., 2015) and more. These tend to be high-
dimensional and computationally challenging imaging problems, increasingly so
at high resolutions. Spherical inverse problems are difficult to solve using
posterior sampling methods such as Markov chain Monte Carlo (MCMC), largely
due to the generally large number of parameters to sample and the high
computational cost of repeated evaluations of the forward problem on the
sphere. MCMC methods have grown in popularity in recent decades (e.g.
Mosegaard & Tarantola, 1995; Malinverno, 2002; Bodin et al., 2012; Lewis &
Bridle, 2002; Corless et al., 2009; Cai et al., 2018a), benefitting from their
ability to sample the full posterior probability density function (pdf), which
constitutes the solution of the inverse problem. This allows for flexible
calculation of any measure of uncertainty in the solution. Furthermore, they
can be used to solve non-linear inverse problems, which are commonplace in,
for example, geophysics (Malinverno, 2002; Bodin et al., 2012; Ferreira et
al., 2020). However, sampling methods come at significant computational cost,
for both solving the inverse problem and for model comparison, which is often
computationally infeasible. The simpler MCMC methods, such as the Metropolis-
Hastings (MH) algorithm, are known to struggle in high-dimensional parameter
spaces (e.g. Roberts & Rosenthal, 1998; Neal, 2012), such as those for
spherical imaging problems. The exponential increase in volume of the
parameter space makes MH unable to converge to a solution. Alternative
gradient-based methods such as the Hamiltonian Monte Carlo (HMC) (Duane et
al., 1987; Neal, 2012) and the unadjusted Langevin algorithm (ULA) (Roberts &
Tweedie, 1996) scale favourably compared to the MH algorithm as the number of
dimensions increases (Roberts & Rosenthal, 1998; Neal, 2012). The caveat here
is that these methods can only be applied to smooth distributions, limiting
the form of prior information that can be used as regularisation. Conversely,
nested sampling methods (Skilling, 2006) can be used for non-smooth problems.
However, most nested sampling approaches cannot scale to the high-dimensional
settings of inverse problems without gradient information. Proximal MCMC
methods have recently been developed (Pereyra, 2016), leveraging proximal
mappings to scale to high dimensions and allow for non-smooth distributions.
Proximal mappings have also recently been used to scale nested sampling
methods to high-dimensional imaging problems (Cai et al., 2022).

Compressed sensing (Donoho, 2006; Candès et al., 2011) has demonstrated that
sparse signals can be accurately recovered from incomplete data, recovering
both sharp and smooth image features simultaneously from underdetermined
systems (Loris et al., 2007; Wallis et al., 2017). As a result, sparse priors
have been widely adopted for solving inverse imaging problems, e.g. for radio
interferometry applications (Wiaux et al., 2009; Carrillo et al., 2014;
Pratley et al., 2017; Cai et al., 2018a, b), cosmological mass-mapping
(Lanusse et al., 2016; Price et al., 2020b), and have received some attention
in seismic tomography (Loris et al., 2007; Charléty et al., 2013). The prior
pdf (i.e. the distribution representing a prior beliefs to constrain the
inverse problem) often used to promote sparsity is the Laplace distribution,
which is non-differentiable and so cannot be used with gradient-based sampling
algorithms.

Current approaches to uncertainty quantification for high-dimensional inverse
problems with sparse a priori constraints exploit results from information
theory to derive approximate credible regions from a single point estimate
solution of the inverse problem (Pereyra, 2017; Cai et al., 2018b; Price et
al., 2020b). The point estimate solution can be obtained using convex
optimisation (Cai et al., 2018b; Price et al., 2020b), and thus this approach
is much faster than posterior sampling, particularly for high-resolution
spherical problems. However, this information theory based approach can only
approximate credible regions. A sampling method on the other hand will have
the flexibility to calculate any measure of uncertainty, as the samples
represent a full probability distribution.

In this work we provide a framework for solving inverse imaging problems on
the sphere using posterior sampling with sparsity-promoting priors. For this,
we leverage a proximal MCMC algorithm (Pereyra, 2016), which is a sampling
scheme that uses proximal mappings to efficiently sample high-dimensional
parameter spaces. Proximal mappings can be viewed as more general gradient
operators that can be used even on non-differentiable functions, allowing us
to use a Laplace prior to promote sparsity. In our framework we outline the
modifications needed to the proximal MCMC algorithm for spherical problems,
with particular consideration given to the spherical parameterisation and
forward operator. We then demonstrate our framework in practice first on a
common problem from the field of global seismic tomography using both
synthetic and real data, and then on a low-resolution full-sky cosmological
mass-mapping example from simulation data.

This paper is structured as follows. Section 2 gives the necessary
mathematical background for Bayesian inference and representations of
spherical images. Section 3 outlines our framework for posterior sampling of
inverse problems on the sphere, with details of the proximal MCMC algorithm
that we use. Sections 4 and 5 contain our illustrative examples from global
seismic tomography and cosmological mass-mapping, respectively, and we
conclude in Section 6.

##  2 Mathematical Background

In this section we provide the necessary mathematical background for this
work, including Bayesian inference for imaging, and harmonic and wavelet
representations of spherical images.

###  2.1 Bayesian inference for imaging

Consider some observed data $\bm{d}$ and model parameters $\bm{m}$ that are
related by some general, possibly non-linear, forward operator $\bm{G}$, as

$\bm{d}=\bm{G}(\bm{m})+\bm{n},$ (1)

where $\bm{n}$ represents observational noise. The aim of the inverse problem
is to infer $\bm{m}$ from $\bm{d}$. This can be formulated in a Bayesian
statistical framework by Bayes Theorem as

$p(\bm{m}{}\bm{|d})\propto p(\bm{d|m}{})p(\bm{m}{}).$ (2)

Here $p(\bm{m}{}\bm{|d})$ is known as the posterior pdf that encapsulates
knowledge of the model parameters $\bm{m}$ given data $\bm{d}$ and represents
the solution to the inverse problem (Tarantola & Valette, 1982; Mosegaard &
Tarantola, 1995). The likelihood function $p(\bm{d|m}{})$, assuming
independent and identically distributed Gaussian noise with known standard
deviation $\sigma$ in the data, is given by

$p(\bm{d|m}{})\propto\mathrm{exp}\left(-\frac{\|\bm{d}-\bm{G}(\bm{m})\|_{2}^{2}}{2\sigma^{2}}\right).$
(3)

Inverse problems are typically ill-posed and often require some form of
regularisation in the form of prior knowledge about the model parameters to
further constrain the inverse problem. In the Bayesian framework, this is
represented by the prior pdf $p(\bm{m}{})$, which can take many forms. Common
examples include the Uniform prior (Jaynes, 2003),
$p(\bm{m}{})\propto\mathrm{constant}$; the Gaussian prior (Tarantola &
Valette, 1982; Golub et al., 1999),
$p(\bm{m}{})\propto\mathrm{exp}(-\mu\|\bm{m}\|_{2}^{2})$; and the Laplace
prior (Donoho, 2006; Candès et al., 2011),
$p(\bm{m}{})\propto\mathrm{exp}(-\mu\|\bm{m}\|_{1})$, to promote sparsity. In
these priors, $\|\cdot\|_{p}$ is the $l_{p}$-norm and $\mu$ is a parameter
that captures the width of the distribution. Henceforth we use the Laplace
prior as we wish to promote sparsity. The constant of proportionality in
equation 2 is known as the Bayesian evidence. This is useful for comparing the
physical model ($\bm{G}$) to alternative hypotheses, but may be ignored for
inference about model parameters $\bm{m}$.

Solving the inverse problem by probabilistic sampling involves sampling the
posterior pdf $p(\bm{m}{}\bm{|d})$, most commonly using MCMC algorithms. The
simplest is the MH algorithm. MH proposes a new sample from a proposal
distribution, which has a certain probability of being accepted. More
sophisticated and better suited to high-dimensional problems than MH are
gradient-based methods, such as the HMC and Langevin algorithms mentioned
previously, which use the gradients of the target pdf to more efficiently
guide the exploration of the parameter space. While the set of samples from
the posterior pdf represents the full solution of the inverse problem, results
are generally reported in terms of summary statistics, such as the mean of the
samples. The maximum a posteriori (MAP) solution,

$\bm{m}^{\mathrm{(MAP)}}=\operatorname*{argmax}_{\bm{m}}\;p(\bm{m}{}\bm{|d})$
(4)

is also a common choice, however, due to the randomness of MCMC there is no
guarantee of the MAP being found. The highest posterior probability sample of
the MCMC provides the closest estimate of the MAP, but can still be relatively
far from the true MAP, particularly in higher dimensional spaces.
Additionally, some measure of uncertainty on individual model parameters can
be calculated from the posterior samples. We give further details of this in
the next section.

In inverse imaging problems, the model parameters $\bm{m}$ represent a 2D (or
possibly 3D) image. MCMC can be used to sample either the direct pixel values
$\bm{x}$ of the image, or the coefficients $\bm{\alpha}$ representing the
image in a particular basis. The forward operator $\bm{G}$ in equation 1 can
be described as a (possibly non-linear) measurement operator $\bm{\Phi}$
acting on the image $\bm{x}$. In many physical situations, the image $\bm{x}$
has a sparse representation in some basis, which can be exploited to
regularise the inverse problem. For some sparsifying basis $\bm{\Psi}$, we
have $\bm{x}=\bm{\Psi\alpha}$, where $\bm{\alpha}$ is the set of coefficients
representing $\bm{x}$ in the basis encoded in $\bm{\Psi}$. The model
parameters $\bm{m}$ sampled by MCMC are the coefficients $\bm{\alpha}$, and as
such the forward operator is given by

$\bm{G}(\bm{\alpha})=\bm{\Phi}(\bm{\Psi\alpha}).$ (5)

For non-linear inverse problems, i.e. where $\bm{\Phi}$ is a non-linear
measurement operator, a common choice is to linearise via a Taylor expansion
around an initial guess model. While this usually leads to acceptable results
in weakly non-linear problems, linearisation errors and the need for a good
initial guess model can lead to artefacts in the recovered model. If the
forward problem (equation 5) is computationally fast, sampling methods are
typically best for strongly non-linear problems.

###  2.2 Harmonic representations of spherical images

We consider functions $f(\theta,\phi)$ defined on the sphere $\mathbb{S}^{2}$
with colatitude $\theta\in[0,\pi]$ and longitude $\phi\in[0,2\pi)$. The
canonical orthogonal basis functions for $\mathbb{S}^{2}$ are the spherical
harmonics

$Y_{lm}(\theta,\phi)=\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}}P_{\ell}^{m}(\cos\theta)e^{im\phi}$
(6)

for non-negative integers $\ell$ and integers $m\leq|\ell|$, denoting angular
degree and angular order, respectively. As per the Condon-Shortley phase
convention, the associated Legendre functions $P_{\ell}^{m}(x)$ include a
$(-1)^{m}$ phase factor, ensuring the conjugate symmetry relation

$Y_{\ell,-m}(\theta,\phi)=(-1)^{m}Y_{\ell m}^{*}(\theta,\phi),$ (7)

where ∗ denotes complex conjugation. Any square integrable scalar function on
$\mathbb{S}^{2}$ can be represented in the frequency domain by its harmonic
coefficients $f_{lm}\in\mathbb{C}$, obtained by projecting $f$ onto the
spherical harmonic basis functions using the inner product (forward spherical
harmonic transform)

$f_{lm}=\langle f,Y_{\ell m}\rangle=\int_{\mathbb{S}^{2}}f(\theta,\phi)Y_{\ell
m}^{*}(\theta,\phi)\sin\theta d\theta d\phi.$ (8)

The function $f(\theta,\phi)$ can be recovered exactly from its spherical
harmonic coefficients (inverse spherical harmonic transform) by

$f(\theta,\phi)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}f_{\ell m}Y_{\ell
m}(\theta,\phi).$ (9)

This formalism can be generalised to spin-$s$ fields (Newman & Penrose, 1966;
Goldberg et al., 1967), common to astrophysical spherical signals (McEwen et
al., 2015b; Wallis et al., 2021), which are characterised by local rotations
$\chi\in[0,2\pi)$ in the tangent plane at a point $(\theta,\phi)$,

${}_{s}f^{\prime}(\theta,\phi)=e^{-is\chi}\ _{s}f(\theta,\phi)$ (10)

The spin-weight $s$ of a function can be increased or decreased by applying
the spin-raising or spin-lowering operators, $\eth,\overline{\eth}$,
respectively.

$\eth\equiv-\sin^{s}\theta\left(\frac{\partial}{\partial\theta}+\frac{i}{\sin\theta}\frac{\partial}{\partial\phi}\right)\sin^{-s}\theta,$
(11)
$\overline{\eth}\equiv-\sin^{-s}\theta\left(\frac{\partial}{\partial\theta}-\frac{i}{\sin\theta}\frac{\partial}{\partial\phi}\right)\sin^{s}\theta.$
(12)

Note how these operators are defined based on the spin of the function to
which they are applied. Using these operators, one can obtain the spin-
weighted spherical harmonics from the original scalar (spin-0) spherical
harmonics, which form an orthogonal basis for spin functions on the sphere

${}_{s}Y_{\ell m}(\theta,\phi)=\begin{cases}0,&\ell<|s|\\\
\sqrt{\frac{(\ell-s)!}{(\ell+s)!}}\eth^{s}Y_{\ell m}(\theta,\phi),&0\leq
s\leq\ell\\\
\sqrt{\frac{(\ell+s)!}{(\ell-s)!}}(-1)^{s}\overline{\eth}^{-s}Y_{\ell
m}(\theta,\phi),&-\ell\leq s\leq 0\\\ \end{cases}$ (13)

As such, spin-$s$ fields can easily be decomposed and reconstructed in terms
of spin spherical harmonics in a manner similar to equations 8 and 9.

In this work we consider bandlimited signals, that is signals such that
${}_{s}f_{\ell m}=0,\;\forall\;\ell\geq L$ for some bandlimit $L$. We
discretise spherical signals according to the McEwen-Wiaux (MW) sampling
theorem (McEwen & Wiaux, 2011). This equiangular sampling theorem has
theoretically exact and efficient spherical harmonic transforms and requires
fewer discrete points on the sphere than other sampling theorems (e.g.
Driscoll & Healy, 1994; Górski et al., 2005). It also exploits a relationship
between spin spherical harmonics and Wigner functions (Goldberg et al., 1967;
McEwen & Wiaux, 2011) to avoid repeated applications of spin-raising/lowering
operators for fast spin spherical harmonic transforms.

###  2.3 Wavelet representations of spherical images

For our sparsifying basis, we consider scale-discretised axisymmetric wavelets
(Wiaux et al., 2008; Leistedt et al., 2013; McEwen et al., 2018), i.e., the
wavelets are azimuthally symmetric when placed at the poles. There exist
directional wavelets which vary azimuthally (e.g. Wiaux et al., 2008; McEwen
et al., 2015b, 2018) and may be desirable in certain applications. However,
such wavelets require transforms on the rotation group $\mathrm{SO}(3)$ which
are more computationally expensive than transforms on $\mathbb{S}^{2}$ (McEwen
et al., 2015a). The spherical wavelet transform is defined as the convolution
of $f$ and the wavelets $\Psi^{j}(\theta,\phi)$. Convolution on the sphere is
defined as the inner product of $f$ with wavelets that have been rotated over
the surface of the sphere by some operator $\mathcal{R}_{(\theta,\phi)}$,
analogous to standard convolution on the plane where the wavelet is laterally
translated over the image. As such, the wavelet coefficients
$W^{\Psi^{j}}\in\mathbb{S}^{2}$ are given by

$W^{\Psi^{j}}(\theta,\phi)=\langle
f,\mathcal{R}_{(\theta,\phi)}\Psi^{j}\rangle.$ (14)

The wavelets $\Psi^{j}$ cover a range of scales $J_{0}\leq j\leq J$, which
extract the highest frequency (high $\ell$ information of $f$). The lowest
frequency information (low $\ell$) is extracted by a scaling function
$\Upsilon(\theta,\phi)$ in a similar manner:

$W^{\Upsilon}(\theta,\phi)=\langle
f,\mathcal{R}_{(\theta,\phi)}\Upsilon\rangle.$ (15)

The wavelets and scaling function are defined on the harmonic line such that
the wavelets can be seen as narrow pass-band filters and the scaling function
a low-pass filter that includes the zero frequency component. For more details
on choosing the width of the wavelets on the harmonic line we refer the
interested reader to Leistedt et al. (2013). The wavelets and scaling function
are defined such that they satisfy an admissibility condition which allows $f$
to be decomposed and recovered exactly from its wavelet coefficients. The
reconstruction is given by

$f(\theta,\phi)=\sum_{\Gamma}\int_{\mathbb{S}^{2}}W^{\Gamma}(\theta^{\prime},\phi^{\prime})(\mathcal{R}_{(\theta^{\prime},\phi^{\prime})}\Gamma)(\theta,\phi)\sin{\theta^{\prime}}d\theta^{\prime}d\phi^{\prime}$
(16)

where $\Gamma\in\\{\Upsilon,\Psi^{j}\\}$. Note that in practice the spherical
wavelet transform is computed in harmonic space (Leistedt et al., 2013),
meaning the forward and inverse wavelet transforms (equations 14-16)
implicitly involve both a forward and an inverse spherical harmonic transform
(Wallis et al., 2017), which dominate the computational complexity. The
spherical harmonic and axisymmetric wavelet transforms used in this work are
implemented in the packages PYSSHT111https://github.com/astro-informatics/ssht
and PYS2LET222https://github.com/astro-informatics/s2let, respectively.

##  3 Framework for posterior sampling on the sphere

In this section we outline our framework for sampling the posterior of
spherical inverse problems. We give the details of the high-dimensional
sampling algorithm we use and the modifications necessary for the sphere.
Extra consideration is given to the forward operator, and finally we define a
measure of uncertainty that can be calculated from the posterior samples.

###  3.1 Proximal MCMC

####  3.1.1 Moreau-Yosida envelope and proximal mapping

The $\lambda$-Moreau-Yosida envelope of some concave function $h(\bm{x})$ is
given by (Moreau, 1962; Bauschke & Combettes, 2017)

$h^{\lambda}(\bm{x})=\min_{\bm{u}}\\{h(\bm{u})+\|\bm{u}-\bm{x}\|^{2}_{2}/(2\lambda)\\},$
(17)

where $\lambda>0$. This is a smooth approximation of $h$ with many desirable
properties. Firstly, $h^{\lambda}$ can be made arbitrarily close to $h$ by
making $\lambda$ small. Secondly, the minimisers of $h^{\lambda}$ are the same
as the minimisers of $h$. Thirdly, $h^{\lambda}$ is continuously
differentiable even if $h$ is not. The gradient of $h^{\lambda}$ is given by

$\nabla
h^{\lambda}(\bm{x})=[(\bm{x}-\mathrm{prox}^{\lambda}_{h}(\bm{x}))]/\lambda,$
(18)

where $\mathrm{prox}^{\lambda}_{h}(\bm{x})$ is the proximal mapping of $h$
(Moreau, 1962)

$\mathrm{prox}^{\lambda}_{h}(\bm{x})=\operatorname*{argmin}_{\bm{u}}\\{h(\bm{u})+\|\bm{u}-\bm{x}\|^{2}_{2}/(2\lambda)\\}.$
(19)

Rewriting equation 18 as

$\mathrm{prox}^{\lambda}_{h}(\bm{x})=\bm{x}-\lambda\nabla h^{\lambda}(\bm{x})$
(20)

is reminiscent of the standard first-order finite differences approximation
$h(x_{0}+a)\approx h(x_{0})+ah^{\prime}(x_{0})$, which shows that
$\mathrm{prox}^{\lambda}_{h}$ may be interpreted as a gradient step in
$h^{\lambda}$ with step size $\lambda$. The key property here is that the
proximal map can be used to minimise $h^{\lambda}$ and, by extension, the non-
differentiable $h$ since they share the same minimisers (Parikh & Boyd, 2014).
For an excellent introduction to proximal mappings and algorithms, see Parikh
& Boyd (2014) and Combettes & Pesquet (2011).

####  3.1.2 Proximal Langevin algorithm

The proximal MCMC method developed by Pereyra (2016) is based on Langevin MCMC
(Roberts & Tweedie, 1996), a gradient-based sampling method, which we describe
here. Consider a Langevin diffusion process $Y(t)$ for $0\leq t\leq T$
associated with a stationary pdf $\pi$. This process is given by the
stochastic differential equation

$dY(t)=\frac{1}{2}\nabla\log{\pi[Y(t)]}dt+dW(t)$ (21)

for some Brownian motion $W$. We use the forward Euler discrete time
approximation with step size $\delta$

$\bm{m}^{(i+1)}=\bm{m}^{(i)}+\frac{\delta}{2}\nabla\log{\pi[\bm{m}^{(i)}]}+\sqrt{\delta}\bm{w}^{(i)},$
(22)

where $\bm{m}$ is the discretised Langevin diffusion and
$\bm{w}^{(i)}\sim\mathcal{N}(0,\mathbb{I})$. This is the Unadjusted Langevin
Algorithm (ULA) (Roberts & Tweedie, 1996). Under certain regularity
conditions, ULA produces samples that converge to an ergodic measure close to
$\pi$. A MH acceptance step can be added to remove the approximation error,
giving the Metropolis-Adjusted Langevin Algorithm (MALA) (Roberts & Tweedie,
1996).

A well known limitation of the Langevin and other gradient-based sampling
algorithms is that they require differentiable pdfs, which is not the case for
some popular priors such as the Laplace prior. A proposed solution to this is
to apply a Moreau-Yosida approximation to the non-differentiable terms in the
posterior pdfs (Pereyra, 2016; Cai et al., 2018a; Pereyra et al., 2020). The
posterior pdf for our inverse problem (equation 2) is of the form
$\pi(\bm{m})\propto\mathrm{exp}(-g(\bm{m})-f(\bm{m}))$, where
$g(\bm{m})=\frac{1}{2\sigma^{2}}\|\bm{d}-\bm{G}(\bm{m})\|_{2}^{2}$ is our
Gaussian data fidelity and $f(\bm{m})=\mu\|\bm{m}\|_{1}$ is our non-
differentiable Laplacian prior. Applying a $\lambda$-Moreau-Yosida
approximation to the prior, it then follows that the chain step for the
Moreau-Yosida ULA (MYULA) is given by

$\displaystyle\bm{m}^{(i+1)}$ $\displaystyle{}={}$
$\displaystyle\left(1-\frac{\delta}{\lambda}\right)\bm{m}^{(i)}+\frac{\delta}{\lambda}\mathrm{prox}_{f}^{\lambda}[\bm{m}^{(i)}]$
(23) $\displaystyle{-}\>\delta\nabla
g[\bm{m}^{(i)}]+\sqrt{\delta}\bm{w}^{(i)}.$

The tuning parameter $\delta$ must be small for the forward Euler
approximation (equation 22) to converge, and Pereyra (2016) argued that the
optimal value for $\lambda$ is $\lambda=\delta/2$. We describe how to
calculate the proximal mapping of our prior on the sphere (second term
equation 23) in the following subsection. If $\bm{G}$ is linear, the gradient
of the data fidelity is straightforward to calculate using the adjoint

$\nabla g=\bm{G}^{\dagger}(\bm{Gm}-\bm{d})/\sigma^{2}.$ (24)

It is straightforward to modify equation 24 if the data errors are
characterised by a covariance matrix $\bm{C}$ rather than a single standard
deviation $\sigma$. It is important to note here that for spherical inverse
problems, the forward operator $\bm{G}$ may include spherical harmonic
transforms, and thus the adjoint transforms are also needed for equation 24.
Further discussion about $\bm{G}$ in the spherical setting is given later in
this section. Algorithm 1 outlines the use of the MYULA chain in practice,
highlighting the key steps and equations. Again, a MH acceptance step can be
added (Pereyra, 2016; Cai et al., 2018a), however in our experiments there was
little improvement for the additional computational cost.

Algorithm 1 MYULA on the sphere

INPUTS: observed data $\bm{d}$, data errors $\sigma$, initial sample
$\bm{m}^{(0)}$, $i=0$, $N$, $N_{\mathrm{thin}}$, $N_{\mathrm{burn}}$,
quadrature weights $\bm{q}$, $\delta$, $\lambda$, $\mu$

OUTPUTS: chain $\\{\bm{m}^{(i)}:i=1,...,N\\}$

while $i<N\times N_{\mathrm{thin}}+N_{\mathrm{burn}}$ do

Calculate gradient of data fidelity (eq 24)

Calculate proximity map of prior (eq 26)

Calculate $\bm{m}^{(i+1)}$ (eq 23)

if $i>N_{\mathrm{burn}}$ then

if $\mathrm{mod}(i,N_{\mathrm{thin}})=0$ then

Save $\bm{m}^{(i+1)}$ to chain

end if

end if

$i\mathrel{+}=1$

end while

###  3.2 Modification for the sphere

We sample a set of spherical wavelet coefficients that are defined at each
point on the sphere (equation 14). The MW sampling theorem is equiangular, and
as such we need to account for an overdensity of sampling points near the
poles when we calculate the proximal mapping of our sparsity-promoting prior.
Proximal operators are generally calculated by a small convex optimisation
problem (Parikh & Boyd, 2014). Fortunately, there exist closed-form
representations for the proximal mapping of many common functions (Combettes &
Pesquet, 2011), including the $\ell_{1}$-norm we use in our prior. This is
crucial for MCMC methods where the proximal mapping needs to be repeatedly
computed. It can be shown that for $f(\bm{m})=\mu\|\bm{m}\|_{1}$,

$\mathrm{prox}^{\lambda}_{f}(\bm{m})=\mathrm{soft}_{\lambda\mu}(\bm{m})$ (25)

where $\mathrm{soft}_{\lambda\mu}$ is the soft thresholding operator with
threshold ${\lambda\mu}$ defined as (Combettes & Pesquet, 2011)

$\mathrm{soft}_{\lambda\mu}(m_{i})=\begin{cases}0&\mathrm{if}\;\;m_{i}\leq{\lambda\mu},\\\
m_{i}-{\lambda\mu}\mathrm{sgn}(m_{i})&\mathrm{if}\;\;m_{i}>{\lambda\mu}\end{cases}$
(26)

where $\mathrm{sgn}(x)$ is the sign of $x$, and can be applied component-wise
to the components $m_{i}$ of $\bm{m}$. To account for the spherical sampling
overdensity, a weighting must be applied to parameters $\bm{m}$, which takes
the form of a diagonal matrix $\bm{W}$ of quadrature weights (McEwen & Wiaux,
2011) that vary with colatitude and bandlimit. This can be seen as applying a
component-wise prior on each of the elements of $\bm{m}$. The proximal mapping
changes only in that the threshold in equation 26 becomes $\lambda\mu w_{i}$,
where $w_{i}$ is the $i^{\mathrm{th}}$ diagonal element of $\bm{W}$.

###  3.3 The forward operator

A key consideration for any MCMC sampling algorithm is the speed of the
forward modelling step, i.e. making predictions of the data $\bm{d}$ given a
MCMC sample $\bm{m}$. This is particularly important for spherical problems
where the forward operator $\bm{G}$ may contain spherical harmonic transforms,
which are, unfortunately, slow and typically scale in complexity as
$\mathcal{O}(L^{3})$ (McEwen & Wiaux, 2011). As discussed in Section 2, the
spherical harmonic transforms are implicit in the spherical wavelet transforms
which form our sparsifying basis. Expanding equation 5 to see this, to sample
the spherical wavelet coefficients in pixel space $\bm{\alpha}$ we have in
practice

$\bm{G}(\bm{\alpha})=\bm{\Phi\mathrm{S}}^{-1}\bm{\mathrm{WS}\alpha},$ (27)

where $\bm{\mathrm{S}}$ and $\bm{\mathrm{S}}^{-1}$ are the forward and inverse
spherical harmonic transforms (equations 8 and 9), respectively, and
$\bm{\mathrm{W}}$ is the spherical wavelet transform in harmonic space
(Leistedt et al., 2013; Wallis et al., 2017). This is further complicated by
the need for the adjoints of these transforms (equation 24). These four
spherical harmonic transforms at each MCMC step immediately limits the
bandlimits $L$ for which sampling the posterior may be feasible. For
reference, in our experiments we performed dummy inversions at $L=64$ and
$L=128$ with an identity measurement operator ($\bm{\Phi}=\bm{I}$) which took
around 3 and 17 days, respectively, for $10^{6}$ samples on a 2.5 GHz Intel
Xeon Platinum 8180M processor. While this may not be an issue for applications
where information exists at relatively low degrees, for example in seismic
tomography where maximum bandlimits are typically around $L=40$ (e.g. Ritsema
et al., 2011; Chang et al., 2015), it will be impractical for applications
where the bandlimits of interest are higher (e.g. Price et al., 2020b). To
avoid spherical harmonic transforms in equation 27, one could instead sample
the spherical harmonic coefficients of $\bm{\alpha}$, provided one can
reformulate the measurement operator $\bm{\Phi}$ appropriately. We note
however that it is possible that, as shown in the following section, that the
harmonic measurement operator may be slower than the pixel space operator even
with spherical wavelet transforms. Further, it is conceivable that, depending
on how the measurement operator scales, the choice of harmonic or pixel space
will depend on the bandlimit $L$. In either case, at least one spherical
harmonic transform will be needed for the prior, as there is no reason to
expect the harmonic representation of the wavelet coefficients to also be
sparse. As such, spherical harmonic transforms are unavoidable, and it is
crucial that these are computed as efficiently as possible, for example via
repeated exploitation of fast Fourier transforms as in McEwen & Wiaux (2011).
This all highlights the special consideration that must be given to the
forward operator for spherical inverse problems.

###  3.4 Spherical wavelet parameterisation

As discussed in Section 2, our inverse problem is parameterised using
axisymmetric spherical wavelets (Wiaux et al., 2008; Leistedt et al., 2013;
McEwen et al., 2015b) to promote sparsity. We can exploit the multiresolution
character of these wavelets for further computational savings. By
construction, the wavelets at each scale $j$ each have different bandlimits
$k_{j}=B^{j+1}\leq L$, where $B$ is a wavelet scale parameter (Leistedt et
al., 2013). By using a sampling theorem, the transforms at each spherical
wavelet scale can be performed up to their own bandlimit $k_{j}$, and only the
minimum number of samples on the sphere at that bandlimit are needed. This
multiresolution transform gives a four to five times speed up of the spherical
wavelet transforms (Leistedt et al., 2013), and also dramatically reduces the
dimensionality of our parameter space, compared to a full resolution transform
where each wavelet scale is sampled at the overall bandlimit $L$. As an
example, for parameters $L=32,\;B=2,\;J_{0}=2$, the full resolution algorithm
has $10\,080$ wavelet coefficients, compared to only $4676$ for the
multiresolution algorithm. While a full assessment of how this affects the
convergence speed of the MCMC chain in terms of number of required steps is
beyond the scope of this paper, we found real-time speed-ups and significant
memory savings in our experiments.

###  3.5 Uncertainty quantification

By collecting samples from the posterior we can calculate any measure of
uncertainty. For example, a common choice in Bayesian statistics is the
credible intervals $[\xi_{i}^{-},\xi_{i}^{+}]$ of the model parameters. These
intervals contain the values that can be taken by parameter $m_{i}$ with
probability $(1-\alpha)$, for some chosen small $\alpha$

$p(m_{i}\in[\xi_{i}^{-},\xi_{i}^{+}]|\bm{d})=1-\alpha$ (28)

The lower and upper interval limits are calculated as the $\frac{\alpha}{2}$
and $1-\frac{\alpha}{2}$ quantiles, respectively, of the posterior. We note
that having sampled the spherical wavelet coefficients, inverse spherical
wavelet transforms will be necessary to obtain an uncertainty map in real
space as opposed to wavelet space. This can be expensive for the same reasons
as previously discussed with respect to the forward operator, although to a
much lesser extent after burn-in and thinned samples have been discarded. Of
course if the desired summary statistic is linear in the model parameters
(e.g. mean) then this can be calculated in wavelet space and only requires a
single spherical wavelet transform. Importantly, in this way we get
uncertainties at the pixel level when viewed in real space. As previously
discussed, current uncertainty quantification methods in similar contexts only
work on superpixels (Price et al., 2020a).

##  4 Seismological surface wave phase velocity maps

The deep structures of the Earth’s interior are best illuminated from seismic
waves. Waves generated from earthquakes travel through the Earth, interacting
with the materials along the way. Inverting the measurements of these waves
recorded at seismic receivers, typically located at the Earth’s surface,
reveals the layers and structures in the Earth (Rawlinson et al., 2010).
Depending on the type and frequency of the waves, this can be done to image a
range of scales from local applications (e.g. imaging oil, gas, geothermal
fields) to global problems (e.g. imaging the whole Earth’s mantle, down to
${\sim}$$2800\text{\,}\mathrm{km}$ depth). The seismic imaging the Earth’s
interior in three dimensions is known as seismic tomography.

In this section we introduce the common problem in global seismic tomography
of building surface wave phase velocity maps, which we use as an illustrative
example. Being a 2-D problem that can be described by relatively simple
theory, this is a natural example of application that is well suited to
illustrate our framework for sampling the posterior for spherical inverse
problems. The resolution requirements in global seismic tomography are also
much lower than for full-sky mass-mapping, as discussed in the next section,
typically only requiring $L\leq 40$. This makes our method directly applicable
to this problem.

###  4.1 Surface wave phase velocity maps

Mapping the phase velocity of surface waves is a common problem in seismic
tomography (Trampert & Woodhouse, 1995; Ekström et al., 1997; Ekström, 2011).
Seismic surface waves generated from earthquakes travel along the Earth’s
surface. These waves are dispersive in that their velocity depends on the wave
period, with different wave periods being sensitive to the structures at
different depths in the Earth’s interior (Dahlen & Tromp, 1998). Phase
velocity maps show how the velocity of surface waves at a given period varies
due to lateral heterogeneities in the Earth’s composition and temperature.
Hence, creating these maps for waves of different periods is often a first
step towards building 3D models of the Earth’s mantle (e.g. Durand et al.,
2015).

For a particular wave period, one can measure from a seismogram the average
phase velocity along the path that the wave has travelled between a seismic
source and receiver. Inverting measurements from many crossing paths, that
would ideally cover the Earth uniformly, produces a phase velocity map. In
practice, the distribution of paths is determined by the locations of
earthquakes, typically along tectonic plate boundaries, and the locations of
seismic stations predominantly in the continents of the northern hemisphere,
making this an ill-posed inverse problem. On global scales, long-period
($T>25$ s) phase velocity measurements are typically modelled using linearised
ray theory (an infinite frequency approximation analogous to geometrical
optics), also known as the great-circle approximation (e.g. Woodhouse &
Dziewonski, 1984; Parisi & Ferreira, 2016) and inverted using least-squares
algorithms (e.g. Tarantola & Valette, 1982; Trampert & Woodhouse, 1995;
Ekström, 2011; Durand et al., 2015). In this framework, the path travelled by
the seismic wave is assumed to correspond to the great circle between the
source and the receiver. The observed mean phase velocity anomaly
$\langle{\delta c/c_{0}}\rangle$ for a given source-receiver pair is given by
the average of the phase velocity field along the minor arc great circle $S$,

$\left\langle{\frac{\delta
c}{c_{0}}}\right\rangle=\frac{1}{\Delta}\int_{S}\frac{\delta
c}{c_{0}}(\theta,\phi)dS\,,$ (29)

where $c_{0}$ is the phase velocity value computed for a reference Earth
model, $\delta c=c-c_{0}$, with $c$ being the actual phase velocity, and
$\Delta$ is the length of path $S$.

Although there exist more complete theories describing phase velocity
anomalies, notably including non-linear and finite frequency effects, recent
studies have shown that the great-circle approximation accurately predicts the
phase of long-period fundamental mode surface waves for current global
tomography models (e.g. Parisi & Ferreira, 2016; Godfrey et al., 2019). Hence,
in this study we focus on the great-circle approximation. For simplicity, we
also only invert for isotropic lateral variations in phase velocity. We note,
however, that our proximal MCMC approach can in principle handle these
additional physical parameters (e.g. anisotropy, the directional dependence of
wave velocity) with some slight modification.

###  4.2 Pixel space path integration

Figure 1: Time taken to perform the forward modelling in pixel space (red), in
pixel space with a spherical wavelet transform (green), and in harmonic space
(blue). The point for the harmonic approach at bandlimit $L=64$ is well beyond
the vertical scale at around 0.2 s. Timings performed on a 2020 MacBook Pro
with an Apple M1 processor.

The common way to compute path integrals on the sphere (equation 29) is to
rotate the coordinate system such that the path lies along the equator, which
is easily expressed in terms of Wigner-D matrices and the spherical harmonic
coefficients of the spherical signal (Dahlen & Tromp, 1998). For a dataset
with $N_{\mathrm{paths}}$ paths, the measurement operator can be represented
by a dense matrix $\bm{\Phi}_{h}\in\mathbb{C}^{N_{\mathrm{paths}}\times
L^{2}}$ acting on the spherical harmonic coefficients
$\hat{\bm{x}}\in\mathbb{C}^{L^{2}}$. For seismic datasets typically consisting
of $\mathcal{O}(10^{5})$ paths, this dense matrix multiplication can be quite
slow. Our approach is to instead measure the path integral directly on the
pixelised sphere using a sparse matrix
$\bm{\Phi}_{p}\in\mathbb{R}^{N_{\mathrm{paths}}\times N_{\mathrm{pixels}}}$,
where each element of $\bm{\Phi}_{p}$ is a weight representing the normalised
distance each path travels in a pixel. This effectively approximates the
integral as a weighted Riemann sum over the pixelised function
$\bm{x}\in\mathbb{R}^{N_{\mathrm{pixels}}}$. The adjoint operator is trivially
the transpose of $\bm{\Phi}_{p}$. The first step for our numerical path
integration is to find discrete geographical points along the great circle
minor arc between a source and a receiver using spherical trigonometry. For
the discretisation, a sampling rate of about 200 points per radian (3.5 points
per degree) was generally sufficient for this work. Each of the geographical
points along the path is then mapped to its nearest MW sampling point. This
mapping assigns a weight to each MW sampling point, which is given by

$w_{tp}=\frac{n_{tp}}{s\Delta},$

where $n_{tp}$ is the number of geographical points on the path that are
closest to MW sampling point indexed in the $\theta$ and $\phi$ directions by
$t$ and $p$, respectively, $s$ is the path sampling rate and $\Delta$ is the
full path length. This can easily be done for each path of the dataset in
parallel to build the full measurement operator.

Table 1: Accuracy of pixel space path integration Bandlimit $L$ |  Mean Diff. (%) | R2E
---|---|---
20 | -0.01 | $2.14\text{\times}{10}^{-4}$
28 | -0.02 | $1.52\text{\times}{10}^{-4}$
32 | -0.02 | $1.34\text{\times}{10}^{-4}$
64 | -0.01 | $5.64\text{\times}{10}^{-5}$

Figure 1 compares the forward modelling time for both the harmonic and pixel
space path integrations, $\bm{\Phi}_{h}$ and $\bm{\Phi}_{p}$, respectively,
for bandlimits $L\in\\{20,28,32,64\\}$ and a realistic set of ray paths (see
Figure 3 bottom). Also shown are the times for the pixel space integration
with an additional spherical wavelet transform $\bm{\Phi}_{p}\bm{\Psi}$, as
required when sampling wavelet coefficients (equation 5) instead of sampling
the image directly. Clearly integration in pixel space is much faster than in
harmonic space, even with the computational overhead of the spherical wavelet
transforms. Pixel space integration also scales better to higher bandlimits.
This is due to the extreme sparsity (less than 2% nonzero elements) of
$\bm{\Phi}_{p}$, whereas $\bm{\Phi}_{h}$ is generally dense. Table 1 shows the
mean percentage difference and the relative squared error
$\mathrm{R2E}=\|\bm{d}_{\mathrm{harm}}-\bm{d}_{\mathrm{pix}}\|_{2}^{2}/\|\bm{d}_{\mathrm{harm}}\|_{2}^{2}$
between predictions made in harmonic space, $\bm{d}_{\mathrm{harm}}$, and in
pixel space, $\bm{d}_{\mathrm{pix}}$, when performed on the ground truth map
we use in our synthetic experiment (see Figure 2 top). As can be expected, as
the bandlimit increases the error in the pixel space integration decreases.
Crucially, even at relatively low bandlimits the pixel space integration is
sufficiently accurate.

###  4.3 Results and Discussion

In this section we present the results of a synthetic test and real data
inversions. We use MYULA to sample the axisymmetric wavelet coefficients
$\bm{\alpha}$ of the spherical image $\bm{x}$, using the sparse measurement
operator described in the previous section.

####  4.3.1 Synthetic experiment

Figure 2: Synthetic GDM52 recovery experiment. (top) Ground truth from Ekström
(2011). (Middle) Mean solution from proximal MCMC. (Bottom) difference between
the truth and our solution. All the maps show perturbations in phase velocity
($\delta c/c_{0}$) with respect to the reference model PREM (Dziewonski &
Anderson, 1981). Green lines show the tectonic plate boundaries (Bird, 2003).
The bottom panel shows the unnormalised posterior probability throughout the
MCMC chain, indicating that the sampling has converged.

As a synthetic example, we use the global phase velocity model GDM52 (Ekström,
2011) at a wave period of $T=40$ s as a ground truth $\bm{x}$. Surface waves
of this period are mainly sensitive to Earth structure at depths of around 100
km (Dahlen & Tromp, 1998; Durand et al., 2015), so the image (Figure 2) shows
well-known tectonic features such as slow anomalies along spreading ridges
(e.g. Ekström, 2011). The model is bandlimited to $L=28$, corresponding to
$3724$ wavelet coefficients using the multiresolution wavelet transform. We
create a synthetic dataset $\bm{d}$

$\bm{d}=\bm{\Phi x}+\bm{n},$ (30)

where $\bm{\Phi}$ is our pixel-space forward operator and
$\bm{n}\sim\mathcal{N}(0,\sigma)$, where $\sigma$ is the standard deviation of
the predictions $\bm{\Phi x}$ to simulate noise in observed data. $\bm{\Phi}$
is constructed using the same paths as those used to originally build GDM52
(Ekström et al., 1997; Ekström, 2011), thereby ensuring a realistic and non-
uniform spatial distribution of the data. In this case we have $179\,657$
paths. We use the signal-to-noise ratio

$\mathrm{SNR}(\bm{x}_{0})=20\log_{10}\left(\frac{\|\bm{x}\|_{2}}{\|\bm{x}-\bm{x}_{0}\|_{2}}\right)$
(31)

and relative squared error

$\mathrm{R2E}(\bm{x}_{0})=\frac{\|\bm{d}-\bm{\Phi\bm{x}}_{0}\|_{2}^{2}}{\|\bm{d}\|_{2}^{2}}$
(32)

to assess the reconstruction accuracy and predictive accuracy, respectively,
of our chosen point solution $\bm{x}_{0}$, which in this case we choose to be
the post-burn mean mean of our MCMC samples. We perform $10^{6}$ MCMC steps,
saving every 500th sample. The first half of the saved samples are discarded
as a burn-in when calculating our mean solution and uncertainty. This takes
the sampling well beyond the point of convergence, which we take to be the
point where there is no longer a significant change in posterior probability
(see Figure 2). This inversion takes about 20 hours on a 2.5 GHz Intel Xeon
Platinum 8180M processor. The tuning parameters are set to $\mu=500$ and
$\delta=10^{-6}$. Figure 2 shows the ground truth, the post-burn mean of our
proximal MCMC samples and the difference between the two. Our solution has an
excellent data fit ($\mathrm{R2E}=$ $9.96\text{\times}{10}^{-3}$) and model
recovery ($\mathrm{SNR}=8.81$ dB). Differences between the ground truth and
our mean solution are small on average (0.5%) with some small-scale blobs of
larger differences. The majority of these blobs occur in the southern
hemisphere where data coverage is poorer (see Figure 3).

We show the map of 95% credible interval ranges as well as a map showing the
density of ray paths of our data set in Figure 3. Here the uncertainty is
calculated on the image space representation of our solution (i.e. $\bm{x}$).
As can be expected, on the whole we find lower uncertainties where we have a
higher density of data in areas such as eastern Asia and western US. We also
see much smaller scale regions of higher uncertainty. Looking at the
uncertainties in wavelet space (i.e. $\bm{\alpha}$, which is sampled by the
proximal MCMC) in Figure 4 it is clear that the smaller scale wavelet
coefficients have higher uncertainty. Thus, the patterns of differences
between our chosen point solution and the truth (Figure 2) are captured by the
uncertainty of our sampled parameters.

Figure 3: Synthetic GDM52 recovery experiment. (top) 95% credible interval
range (equation 28) calculated on our MCMC samples in image space ($\bm{x}$)
in units of percentage deviation from the reference value, as in Figure 2.
(bottom) Row sum of path matrix $\bm{\Phi}$ (unitless) representing the
density of rays in the dataset. Figure 4: Synthetic GDM52 recovery experiment.
95% credible interval ranges of spherical wavelet coefficients at different
scales. The range of the colourbar is the same for all maps in units of
percentage deviation from the reference value, as in Figure 2 and 3.

####  4.3.2 Real Data Inversions

To demonstrate our method on real data, we invert the same data that was used
to build GDM52 (Ekström et al., 1997; Ekström, 2011) at wave periods
$T=25,\;40,\;75$ s, which have strong sensitivity to about 100 km, 160 km and
300 km depths, respectively (Dahlen & Tromp, 1998; Durand et al., 2015). The
results of these inversions are shown in Figure 5. The datasets at these wave
periods consist of about $103\,633$, $179\,657$ and $286\,302$ paths,
respectively. Again we use a bandlimit of $L=28$. For these inversions, we
perform $750\,000$ chain steps, saving every $500^{\mathrm{th}}$ sample and
discarding the first 500 saved samples as burn-in. This takes between 8 and 33
hours, depending on the number of ray paths, on a 2.5 GHz Intel Xeon Platinum
8180M processor. The tuning parameter $\delta$ is chosen on a case-by-case
basis. Our mean solutions show all the expected velocity anomalies, being very
similar to the GDM52 phase velocity maps (Figure 5, left). For example, the
$T=25$ s map (Figure 5, top) shows a clear distinction between the slow
continents and the fast oceans. On the other hand, the $T=40$ s map (Figure 5,
middle) depicts a good correlation between slow anomalies and plate
boundaries, while the $T=75$ s map (Figure 5, bottom) shows deeper mantle
signals, such as high velocities associated with cratons. A notable difference
with GDM52 (Figure 5, left) is a north-south streak of fast velocities off the
coast of the western US seen at all wave periods. This is a well-known
artefact resulting from not modelling azimuthal anisotropy (Forsyth, 1975;
Ekström, 2011). This streak corresponds to a region of high uncertainty in our
solutions (Figure 5, right). Again, the uncertainty maps correlate with ray
density as expected.

Figure 5: Real data inversions. Phase velocity maps from GDM52 (Ekström, 2011)
(left), mean proximal MCMC solutions from this study (middle) and credible
interval ranges (right) for wave periods 25, 40 and 75 s (top to bottom). All
the maps show perturbations in phase velocity ($\delta c/c_{0}$) with respect
to the reference model PREM (Dziewonski & Anderson, 1981).

We emphasise that the main purpose of this study is not to build improved
phase velocity maps, which can be constructed quickly using, e.g. least
squares approaches (Tarantola & Valette, 1982; Trampert & Woodhouse, 1995),
but rather to illustrate and validate our framework for sampling the posterior
of spherical inverse problems with a useful, well-known first application.
Hence, we do not consider more sophisticated theoretical formulations, such
as, e.g. full ray theory (Ferreira & Woodhouse, 2007), finite frequency theory
(Zhou et al., 2005), including anisotropic effects (Ekström, 2011), etc.
Future work will focus on the application of the method to more sophisticated
problems, such as e.g. depth inversions using non-linear theory.

##  5 Low-resolution Cosmological Mass-Mapping

One of the predictions of Einstein’s theory of general relativity was that the
gravitational influence of massive objects will cause light to bend around
them (Einstein, 1905). As a result, distant objects often look distorted when
observed by astronomers. This is the phenomenon of gravitational lensing.
Light from distant stars and galaxies travels through the universe, bending
around all the mass along the way to Earth. Inverting measurements of the
distortions in the images of the light source reveals the hidden masses that
the light passed by (e.g. Kaiser & Squires, 1993), including dark matter
(Heavens, 2009), producing so-called cosmological mass maps.

In this section we demonstrate our framework on the retrieval of cosmological
mass-maps from simulated data. This is a 2-D problem described by a simple
linear forward model. We treat this as a simple demonstration and perform the
inversion at a much lower bandlimit than would normally be required for
cosmological mass-mapping. Due to the poor scaling of the spherical wavelet
transforms with $L$, performing the millions of transforms required at
$L\sim\mathcal{O}(10^{3})$ (Wallis et al., 2021; Price et al., 2020b) would
simply be too slow. However, with future computational improvements to speed
up the transforms (e.g. by exploiting GPUs), our framework will be well-suited
for this particular problem.

###  5.1 Mass-mapping on the celestial sphere

Images of galaxies are typically distorted, as the light they emit gets lensed
by the mass between us and the source galaxies. Gravitational lensing occurs
regardless of the nature of the intervening mass, and as such lensing is an
excellent probe for dark matter (Heavens, 2009). Mass-mapping maps the total
density perturbation along a line of sight between a source galaxy and the
observer based on measurements of the distortion of galaxy images. Up until
recently, lensing surveys only covered relatively small sky-fractions, so
planar approximations were made. As the area of coverage has increased with
newer surveys, planar approximations are no longer valid (Wallis et al.,
2021), resulting in mass-maps now being constructed on the sphere (e.g.
Jeffrey et al., 2021; Price et al., 2020b; Wallis et al., 2021).

Gravitational lensing studies deal with two main fields: the convergence field
${}_{0}\kappa(\theta,\phi)$, which causes magnification of the galaxy image;
and the shear field ${}_{2}\gamma(\theta,\phi)$, which causes rotation and
stretching of the galaxy image (Dodelson, 2017). Mass maps reveal the
convergence field ${}_{0}\kappa(\theta,\phi)$, which can be shown to be the
integrated mass density along the line of sight (Bartelmann & Schneider,
2001). This is linearly related in spherical harmonic space to the spin-2
shear field ${}_{2}\gamma(\theta,\phi)$, measured from observations of galaxy
shapes, by a linear kernel $\mathcal{K}_{\ell}$ given by (Kaiser & Squires,
1993; Wallis et al., 2021)

$\mathcal{K}_{\ell}=\frac{-1}{\ell(\ell+1)}\sqrt{\frac{(\ell+2)!}{{(\ell-2)!}}}.$
(33)

Thus, our forward model in this case is given by
$\bm{\Phi}=\bm{M}\;_{2}\bm{S}^{-1}\bm{K}\;_{0}\bm{S}$, where $\bm{K}$ encodes
the linear kernel above and $\bm{M}$ is a masking matrix to account for areas
on the sky without reliable data (e.g. the galactic plane and ecliptic).
${}_{2}\bm{S}$ and ${}_{0}\bm{S}$ denote the spin-2 and spin-0 spherical
harmonic transforms needed for to account for the spin symmetries of the shear
and convergence fields, respectively. From hereon in we drop the spin
subscripts for clarity. We note that in reality there exists a degeneracy
between $\gamma$ and $\kappa$, and as a result the true observable is not the
shear but the reduced shear $g=\gamma/(1-\kappa)$. In the weak lensing regime,
this non-linear effect is very small, so we ignore it here. Accounting for
non-linearities in our MCMC is possible, provided the relevant gradients of
the non-linear forward model can be computed efficiently.

The choice of using sparsity promoting priors in this case is motivated by the
need to recover non-Gaussian structures, particularly at high $\ell$. These
are created by non-linear structure growth of the density field throughout the
evolution of the Universe. Promoting sparsity has been used previously (e.g.
Leonard et al., 2014; Lanusse et al., 2016; Price et al., 2020b; Starck et
al., 2021) to incorporate non-Gaussian structures, as have Wiener filters on a
Gaussian prior (e.g. Jeffrey et al., 2018),log-normal priors (e.g. Böhm et
al., 2017; Fiedorowicz et al., 2022) and physically informed power spectrum
priors (e.g. Porqueres et al., 2021).

A Bayesian sampling method for mass-mapping was recently implemented by
Fiedorowicz et al. (2022). They used HMC to efficiently navigate the large
parameter space with a log-normal prior. In comparison with our framework, the
log-normal prior is differentiable, which allows the gradient-based HMC to be
used, with the exploitation of efficient autodifferentiation making this
computationally tractable. However, the log-normal prior makes certain
assumptions about the cosmological parameters, and thus the mass maps produced
by their method cannot be used for inference about cosmology, without
additional further work (Fiedorowicz et al., 2022). Our framework makes no
such assumptions.

###  5.2 Results and Discussion

As a synthetic example, we attempt to reconstruct a simulated mass map from
the Takahashi N-body simulation (Takahashi et al., 2017)333Data available at
http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allskyraytracing/.. Slices are
provided at a range of redshifts, and we select redshift slice 16
corresponding to $z\sim 1$. This slice is bandlimited at $L=64$, giving our
ground truth mass map, to which we apply a basic galactic plane and ecliptic
mask.

Synthetic shear data is generated by applying the measurement operator
$\bm{\Phi}$ to the ground truth map. The noise in weak lensing surveys depends
on the number of galaxy counts per unit area and the variance of the intrinsic
ellipticity distribution $\sigma_{e}\sim 0.37$. We choose an overall number
density of galaxies per arcmin2, $n_{\mathrm{gal}}$, and add zero-mean
Gaussian noise to the synthetic data, with the variance of the noise given by

$\sigma^{2}_{t}=\frac{\sigma_{e}^{2}}{\sqrt{2n_{t}}},$ (34)

where $t$ is the colatitude index, and $n_{t}$ is the expected number of
galaxies in a pixel at colatitude $t$ given the overall number density
$n_{\mathrm{gal}}$. This dependence on colatitude comes from the equiangular
nature of the MW sampling theorem (McEwen & Wiaux, 2011).

When reporting the SNR, R2E summary statistics and our solution maps, we use a
slightly larger mask than what was applied to the synthetic shear data. This
is to remove leakage artefacts that occur around the edge of the mask due to
the sudden lack of data and also wavelets that have support both inside and
outside the mask. At the high resolutions typically of interest for mass-
mapping, this leakage will be minimal.

Maps of the ground truth, mean of our MCMC samples, the difference between
them and our measured uncertainty are shown in Figure 6. For this inversion,
we perform $12.5\times 10^{6}$ chain steps, saving every 500th sample and
using the last 3000 samples for our results, taking about 4 days on a 2.5 GHz
Intel Xeon Platinum 8180M processor. The step-size parameter $\delta$ is
chosen to be $10^{-9}$, the largest value that allows a stable Euler
approximation of the Langevin diffusion. The regularisation parameter $\mu$ is
chosen to be $5\times 10^{5}$, to constrain the wavelet coefficients to the
appropriate order of magnitude. At this resolution ($L=64$), there are
$18\,916$ wavelet coefficients to be sampled. The mean solution has an SNR of
7.83 dB and a R2E of 0.1. The uncertainty map, measured as the range of the
95% credible intervals (equation 28) show, as expected, high uncertainty
within the masked regions where the solution is not constrained by any data.

Figure 6: Results of a simple mass-mapping example where the ground truth is
known. (top left) Ground truth simulation from Takahashi et al. (2017). (top
right) Mean of the proximal MCMC samples. (bottom left) Difference between the
top two maps. (bottom right) 95% credible interval range obtained from
proximal MCMC. A galactic plane and ecliptic mask is show in all maps except
the uncertainty map, which, as expected, shows high uncertainty in these
regions.

While we obtain encouraging results here with a good reconstruction and
physically reasonable pixel-level uncertainties, for this method to be adopted
in full-sky mass mapping, computational advances are needed in the
implementation of the spherical harmonic and wavelet transforms such that the
scaling with $L$ is not as severe. Currently, the complexity scaling of the
spherical harmonic transforms is $\mathcal{O}(L^{3})$ (McEwen & Wiaux, 2011),
and dominates over the efficient harmonic space wavelet transform (Leistedt et
al., 2013) discussed in Section 2. With at least two spherical harmonic
transforms per iteration of MYULA, current implementations will not allow for
the desired $L\sim\mathcal{O}(10^{3})$ bandlimits. New efforts implementing
these transforms on GPUs should go a significant way towards pushing our
method to higher resolutions. Further theoretical advances for the
representation of spherical signals could also lead to computational savings
(e.g. Ocampo et al., 2022).

##  6 Conclusions

In this paper we have presented a general framework for posterior sampling of
inverse problems on the sphere with sparsity promoting priors that allows for
flexible uncertainty quantification and extends naturally to non-linear
problems. We demonstrated the potential applicability of this method to both
astrophysical and geophysical problems. As with all MCMC methods, the
suitability of this method depends on the time taken to take the next chain
sample, particularly in relation to the forward modelling step. The
computational complexity of transforms on the sphere mean that this framework
is generally feasible for problems of low to moderate resolution (roughly
$L\leq 64$), such as those commonly considered in global seismic tomography.
At higher resolutions, as needed for full-sky mass-mapping, posterior sampling
quickly becomes intractable largely due to the poor scaling of spherical
harmonic transforms present in the forward operator. In either case, special
consideration must be given to the forward operator and whether it should be
formulated in harmonic or pixel space, and also if its adjoint is known. A
harmonic formulation would avoid repeated spherical harmonic transforms, but
these savings could be lost on the measurement operator. Making use of a more
efficient proximal algorithm based on a stochastic Runge-Kutta approximation
of the Langevin equation (Pereyra et al., 2020) could be a promising way
forward for higher resolution spherical inverse problems. This algorithm is
more complex but converges to a solution much faster than the Euler
approximation algorithm used in this work, thereby potentially requiring fewer
spherical harmonic transforms. Additionally, faster implementations of the
spherical harmonic transforms leveraging GPUs would immediately increase the
potential of our method.

The examples shown in this work are largely illustrative, as simplifications
have been made. In the mass-mapping example, the bandlimit is much lower than
the angular orders at which the convergence spectrum has the most power due to
the computational restrictions imposed by the spherical harmonic transforms.
We also have not fully considered the effect of reduced shear. Nonetheless,
ongoing work to implement the spherical harmonic and wavelet transforms on
accelerators (GPUs) provide a route to scale to higher resolutions, making our
framework is a promising addition to other recent methods for obtaining mass-
maps with uncertainties (e.g. Price et al., 2020b; Fiedorowicz et al., 2022).
For the surface wave tomography example, we have used the great circle
approximation and not accounted for anisotropic effects, in which case a
least-squares approach is fast and efficient. Nonetheless, our results
demonstrate the feasibility of our framework methods on global scale
tomographic inverse problems. In this case, the more commonly used harmonic
formulation of the forward problem proved to be too slow, and we were able to
reformulate it in pixel space such that it was much faster even with the
computational overhead of spherical wavelet transforms. Further, the
uncertainties calculated from our posterior samples make physical sense, being
correlated with the distribution of data. Bayesian methods in seismic
tomography on large-to-global scales have largely been used for independent 1D
inversions (e.g. Shapiro & Ritzwoller, 2002; Khan et al., 2011; Ravenna &
Lebedev, 2017), although new advances in gradient-based or variational methods
(e.g. Fichtner et al., 2018; Gebraad et al., 2020; Zhang & Curtis, 2020; Zhao
et al., 2020) are promising for 2D and 3D probabilistic tomography. Our
framework is a further contribution to this advance in methodology, with the
novelty of being able to use a non-differentiable prior.

## Acknowledgements

A.M. is supported by the STFC UCL Centre for Doctoral Training in Data
Intensive Science (grant number ST/P006736/1). A.M.G.F. is grateful to funding
from the European Research Council (ERC) to the UPFLOW consolidator grant
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 101001601).

## Data Availability

We provide a new Python package,
pxmcmc444https://github.com/auggiemarignier/pxmcmc that contains
implementations of the proximal MCMC algorithms used in this work as well as
the measurement operators, wavelet transforms and priors with their proximal
mappings. The code is designed to be flexible, with base classes that will
allow users to implement their own forward (measurement and transform)
operators and priors. The proximal MCMC algorithms implemented are themselves
not restricted to spherical problems, as the spherical aspects of the
inversions appear in the likelihood and prior proximal calculations. We also
provide scripts and data to reproduce the synthetic experiment described in
this paper. MCMC chains for the real data inversions are available from the
authors upon request. The discretisation of great circle paths is implemented
in our Python code
greatcirclepaths555https://github.com/auggiemarignier/GreatCirclePaths, which
is publicly available.

## References

  * Bartelmann & Schneider (2001) Bartelmann, M. & Schneider, P., 2001.  Weak gravitational lensing, Physics Reports, 340(4), 291–472.
  * Bauschke & Combettes (2017) Bauschke, H. H. & Combettes, P. L., 2017.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer New York, New York, NY.
  * Bird (2003) Bird, P., 2003.  An updated digital model of plate boundaries, Geochem. Geophy. Geosy., 4(3).
  * Bodin et al. (2012) Bodin, T., Sambridge, M., Gallagher, K., & Rawlinson, N., 2012.  Transdimensional inversion of receiver functions and surface wave dispersion, J. Geophys. Res.-Sol. Ea., 117(B2).
  * Böhm et al. (2017) Böhm, V., Hilbert, S., Greiner, M., & Enßlin, T. A., 2017.  Bayesian weak lensing tomography: Reconstructing the 3d large-scale distribution of matter with a lognormal prior, Phys. Rev. D, 96, 123510.
  * Cai et al. (2018a) Cai, X., Pereyra, M., & McEwen, J. D., 2018a.  Uncertainty quantification for radio interferometric imaging - I. Proximal MCMC methods, Mon. Not. R. Astron. Soc., 480(3), 4154–4169.
  * Cai et al. (2018b) Cai, X., Pereyra, M., & McEwen, J. D., 2018b.  Uncertainty quantification for radio interferometric imaging: Ii. map estimation, Mon. Not. R. Astron. Soc., 480(3), 4170–4182.
  * Cai et al. (2022) Cai, X., McEwen, J. D., & Pereyra, M., 2022.  Proximal nested sampling for high-dimensional bayesian model selection, Statistics and Computing, 32(5), 87.
  * Candès et al. (2011) Candès, E. J., Eldar, Y. C., Needell, D., & Randall, P., 2011.  Compressed sensing with coherent and redundant dictionaries, Appl. Comput. Harmon. A., 31(1), 59–73.
  * Carrillo et al. (2014) Carrillo, R. E., McEwen, J. D., & Wiaux, Y., 2014.  PURIFY: a new approach to radio-interferometric imaging, Mon. Not. R. Astron. Soc., 439(4), 3591–3604.
  * Chang et al. (2015) Chang, S.-J., Ferreira, A. M. G., Ritsema, J., van Heijst, H. J., & Woodhouse, J. H., 2015.  Joint inversion for global isotropic and radially anisotropic mantle structure including crustal thickness perturbations, J. Geophys. Res.-Sol. Ea., 120(6), 4278–4300.
  * Charléty et al. (2013) Charléty, J., Voronin, S., Nolet, G., Loris, I., Simons, F. J., Sigloch, K., & Daubechies, I. C., 2013.  Global seismic tomography with sparsity constraints: Comparison with smoothing and damping regularization, J. Geophys. Res.-Sol. Ea., 118(9), 4887–4899.
  * Combettes & Pesquet (2011) Combettes, P. L. & Pesquet, J.-C., 2011.  Proximal Splitting Methods in Signal Processing, pp. 185–212, Springer New York, New York, NY.
  * Corless et al. (2009) Corless, V. L., King, L. J., & Clowe, D., 2009.  A new look at massive clusters: weak lensing constraints on the triaxial dark matter haloes of A1689, A1835 and A2204, Mon. Not. R. Astron. Soc., 393(4), 1235–1254.
  * Dahlen & Tromp (1998) Dahlen, F. A. & Tromp, J., 1998.  Theoretical Global Seismology, Princeton University Press.
  * Dodelson (2017) Dodelson, S., 2017.  Gravitational Lensing, Cambridge University Press, Cambridge.
  * Donoho (2006) Donoho, D. L., 2006.  Compressed sensing, IEEE Trans. Inf. Theory, 52(4), 1289–1306.
  * Driscoll & Healy (1994) Driscoll, J. R. & Healy, D. M., 1994.  Computing Fourier Transforms and Convolutions on the 2-Sphere, Adv. Appl. Math., 15(2), 202–250.
  * Duane et al. (1987) Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D., 1987.  Hybrid monte carlo, Physics Letters B, 195(2), 216–222.
  * Durand et al. (2015) Durand, S., Debayle, E., & Ricard, Y., 2015.  Rayleigh wave phase velocity and error maps up to the fifth overtone, Geophys. Res. Lett., 42(9), 3266–3272.
  * Dziewonski & Anderson (1981) Dziewonski, A. M. & Anderson, D. L., 1981.  Preliminary reference earth model, Phys. Earth Planet. In., 25(4), 297–356.
  * Einstein (1905) Einstein, A., 1905.  Zur elektrodynamik bewegter körper, Annalen der Physik, 322(10), 891–921.
  * Ekström (2011) Ekström, G., 2011.  A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25-250 s, Geophys. J. Int, 187(3), 1668–1686.
  * Ekström et al. (1997) Ekström, G., Tromp, J., & Larson, E. W. F., 1997.  Measurements and global models of surface wave propagation, J. Geophys. Res.-Sol. Ea., 102(B4), 8137–8157.
  * Ferreira & Woodhouse (2007) Ferreira, A. M. G. & Woodhouse, J. H., 2007.  Source, path and receiver effects on seismic surface waves, Geophys. J. Int, 168(1), 109–132.
  * Ferreira et al. (2020) Ferreira, A. M. G., Marignier, A., Attanayake, J., Frietsch, M., & Berbellini, A., 2020.  Crustal structure of the Azores Archipelago from Rayleigh wave ellipticity data, Geophys. J. Int, 221(2), 1232–1247.
  * Fichtner et al. (2018) Fichtner, A., Zunino, A., & Gebraad, L., 2018.  Hamiltonian Monte Carlo solution of tomographic inverse problems, Geophys. J. Int, 216(2), 1344–1363.
  * Fiedorowicz et al. (2022) Fiedorowicz, P., Rozo, E., Boruah, S. S., Chang, C., & Gatti, M., 2022.  Karmma - kappa reconstruction for mass mapping, Monthly Notices of the Royal Astronomical Society, 512(1), 73–85.
  * Forsyth (1975) Forsyth, D. W., 1975.  The early structural evolution and anisotropy of the oceanic upper mantle, Geophys. J. Int, 43(1), 103–162.
  * Gebraad et al. (2020) Gebraad, L., Boehm, C., & Fichtner, A., 2020.  Bayesian Elastic Full-Waveform Inversion Using Hamiltonian Monte Carlo, J. Geophys. Res.-Sol. Ea., 125(3).
  * Godfrey et al. (2019) Godfrey, K. E., Dalton, C. A., Ma, Z., Hjörleifsdóttir, V., & Ekström, G., 2019.  A comparison of approaches for the prediction and inversion of surface wave phase delays, Geophys. J. Int, 217(3), 1496–1514.
  * Goldberg et al. (1967) Goldberg, J. N., Macfarlane, A. J., Newman, E. T., Rohrlich, F., & Sudarshan, E. C. G., 1967.  Spin-s spherical harmonics and ð, Journal of Mathematical Physics, 8(11), 2155–2161.
  * Golub et al. (1999) Golub, G. H., Hansen, P. C., & O’Leary, D. P., 1999.  Tikhonov regularization and total least squares, SIAM J. Matrix Anal. A., 21(1), 185–194.
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., Wandelt, B. D., Hansen, F. K., Reinecke, M., & Bartelmann, M., 2005.  HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, Astrophys. J., 622(2), 759–771.
  * Heavens (2009) Heavens, A., 2009.  Weak lensing: Dark matter, dark energy and dark gravity, Nuclear Physics B - Proceedings Supplements, 194, 76–81.
  * Jaynes (2003) Jaynes, E. T., 2003.  Probability Theory: The Logic of Science, Cambridge University Press, Cambridge.
  * Jeffrey et al. (2018) Jeffrey, N., Abdalla, F. B., Lahav, O., Lanusse, F., Starck, J.-L., Leonard, A., Kirk, D., Chang, C., Baxter, E., Kacprzak, T., Seitz, S., Vikram, V., Whiteway, L., Abbott, T. M. C., Allam, S., Avila, S., Bertin, E., Brooks, D., Carnero Rosell, A., Carrasco Kind, M., Carretero, J., Castander, F. J., Crocce, M., Cunha, C. E., D’Andrea, C. B., da Costa, L. N., Davis, C., De Vicente, J., Desai, S., Doel, P., Eifler, T. F., Evrard, A. E., Flaugher, B., Fosalba, P., Frieman, J., García-Bellido, J., Gerdes, D. W., Gruen, D., Gruendl, R. A., Gschwend, J., Gutierrez, G., Hartley, W. G., Honscheid, K., Hoyle, B., James, D. J., Jarvis, M., Kuehn, K., Lima, M., Lin, H., March, M., Melchior, P., Menanteau, F., Miquel, R., Plazas, A. A., Reil, K., Roodman, A., Sanchez, E., Scarpine, V., Schubnell, M., Sevilla-Noarbe, I., Smith, M., Soares-Santos, M., Sobreira, F., Suchyta, E., Swanson, M. E. C., Tarle, G., Thomas, D., Walker, A. R., & Collaboration, D., 2018.  Improving weak lensing mass map reconstructions using gaussian and sparsity priors: application to des sv, Monthly Notices of the Royal Astronomical Society, 479(3), 2871–2888.
  * Jeffrey et al. (2021) Jeffrey, N., Gatti, M., Chang, C., Whiteway, L., Demirbozan, U., Kovacs, A., Pollina, G., Bacon, D., Hamaus, N., Kacprzak, T., Lahav, O., Lanusse, F., et al., 2021.  Dark energy survey year 3 results: Curved-sky weak lensing mass map reconstruction, Mon. Not. R. Astron. Soc..
  * Kaiser & Squires (1993) Kaiser, N. & Squires, G., 1993.  Mapping the dark matter with weak gravitational lensing, Astrophys. J., 404, 441.
  * Khan et al. (2011) Khan, A., Boschi, L., & Connolly, J. A. D., 2011.  Mapping the Earth’s thermochemical and anisotropic structure using global surface wave data, J. Geophys. Res.-Sol. Ea., 116(B1).
  * Lanusse et al. (2016) Lanusse, F., Starck, J. L., Leonard, A., & Pires, S., 2016.  High resolution weak lensing mass mapping combining shear and flexion, Astronomy & Astrophysics, 591.
  * Leistedt et al. (2013) Leistedt, B., McEwen, J. D., Vandergheynst, P., & Wiaux, Y., 2013.  S2LET: A code to perform fast wavelet analysis on the sphere, Astron. Astrophys., 558, A128.
  * Leonard et al. (2014) Leonard, A., Lanusse, F., & Starck, J.-L., 2014.  Glimpse: accurate 3d weak lensing reconstructions using sparsity, Monthly Notices of the Royal Astronomical Society, 440(2), 1281–1294.
  * Lewis & Bridle (2002) Lewis, A. & Bridle, S., 2002.  Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys. Rev. D, 66, 103511.
  * Loris et al. (2007) Loris, I., Nolet, G., Daubechies, I., & Dahlen, F. A., 2007.  Tomographic inversion using $\ell$1-norm regularization of wavelet coefficients, Geophys. J. Int, 170(1), 359–370.
  * Malinverno (2002) Malinverno, A., 2002.  Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem, Geophys. J. Int, 151(3), 675–688.
  * McEwen & Wiaux (2011) McEwen, J. D. & Wiaux, Y., 2011.  A novel sampling theorem on the sphere, IEEE Trans. Signal Process., 59(12), 5876–5887.
  * McEwen et al. (2015a) McEwen, J. D., Büttner, M., Leistedt, B., Peiris, H. V., & Wiaux, Y., 2015a.  A novel sampling theorem on the rotation group, IEEE Signal Processing Letters, 22(12), 2425–2429.
  * McEwen et al. (2015b) McEwen, J. D., Leistedt, B., Büttner, M., Peiris, H. V., & Wiaux, Y., 2015b.  Directional spin wavelets on the sphere, IEEE Trans. Sig. Proc., submitted.
  * McEwen et al. (2018) McEwen, J. D., Durastanti, C., & Wiaux, Y., 2018.  Localisation of directional scale-discretised wavelets on the sphere, Applied and Computational Harmonic Analysis, 44(1), 59–88.
  * Moreau (1962) Moreau, J. J., 1962.  Fonctions convexes duales et points proximaux dans un espace hilbertien, CR Hebd. Acad. Sci., 255, 2897–2899.
  * Mosegaard & Tarantola (1995) Mosegaard, K. & Tarantola, A., 1995.  Monte Carlo sampling of solutions to inverse problems, J. Geophys. Res.-Sol. Ea., 100(B7), 12431–12447.
  * Neal (2012) Neal, R. M., 2012.  MCMC using Hamiltonian dynamics, in Handbook of Markov Chain Monte Carlo, pp. 113–162, CRC Press.
  * Newman & Penrose (1966) Newman, E. T. & Penrose, R., 1966.  Note on the bondi-metzner-sachs group, Journal of Mathematical Physics, 7(5), 863–870.
  * Ocampo et al. (2022) Ocampo, J., Price, M. A., & McEwen, J. D., 2022.  Scalable and Equivariant Spherical CNNs by Discrete-Continuous (DISCO) Convolutions.
  * Parikh & Boyd (2014) Parikh, N. & Boyd, S., 2014.  Proximal algorithms, Found. Trends Opt., 1(3), 127–239.
  * Parisi & Ferreira (2016) Parisi, L. & Ferreira, A. M. G., 2016.  Empirical assessment of the validity limits of the surface wave full ray theory using realistic 3-D Earth models, Geophys. J. Int, 205(1), 146–159.
  * Pereyra (2016) Pereyra, M., 2016.  Proximal Markov chain Monte Carlo algorithms, Stat. Comput., 26(4), 745–760.
  * Pereyra (2017) Pereyra, M., 2017.  Maximum-a-posteriori estimation with bayesian confidence regions, SIAM J. Imaging. Sci., 10(1), 285–302.
  * Pereyra et al. (2020) Pereyra, M., Mieles, L. V., & Zygalakis, K. C., 2020.  Accelerating Proximal Markov Chain Monte Carlo by Using an Explicit Stabilized Method, SIAM J. Imaging. Sci., 13(2), 905–935.
  * Planck Collaboration (2016) Planck Collaboration, 2016.  Planck 2015 results - I. Overview of products and scientific results, Astron. Astrophys., 594, A1.
  * Porqueres et al. (2021) Porqueres, N., Heavens, A., Mortlock, D., & Lavaux, G., 2021.  Bayesian forward modelling of cosmic shear data, Monthly Notices of the Royal Astronomical Society, 502(2), 3035–3044.
  * Pratley et al. (2017) Pratley, L., McEwen, J. D., d’Avezac, M., Carrillo, R. E., Onose, A., & Wiaux, Y., 2017.  Robust sparse image reconstruction of radio interferometric observations with purify, Mon. Not. R. Astron. Soc., 473(1), 1038–1058.
  * Price et al. (2020a) Price, M. A., Cai, X., McEwen, J. D., Pereyra, M., & Kitching, T. D., 2020a.  Sparse bayesian mass mapping with uncertainties: local credible intervals, Mon. Not. R. Astron. Soc., 492(1), 394–404.
  * Price et al. (2020b) Price, M. A., McEwen, J. D., Pratley, L., & Kitching, T. D., 2020b.  Sparse Bayesian mass-mapping with uncertainties: Full sky observations on the celestial sphere, Mon. Not. R. Astron. Soc., 500(4), 5436–5452.
  * Ravenna & Lebedev (2017) Ravenna, M. & Lebedev, S., 2017.  Bayesian inversion of surface-wave data for radial and azimuthal shear-wave anisotropy, with applications to central Mongolia and west-central Italy, Geophys. J. Int, 213(1), 278–300.
  * Rawlinson et al. (2010) Rawlinson, N., Pozgay, S., & Fishwick, S., 2010.  Seismic tomography: A window into deep earth, Physics of the Earth and Planetary Interiors, 178(3-4), 101–135.
  * Ritsema et al. (2011) Ritsema, J., Deuss, A., van Heijst, H. J., & Woodhouse, J. H., 2011.  S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophys. J. Int, 184(3), 1223–1236.
  * Roberts & Rosenthal (1998) Roberts, G. O. & Rosenthal, J. S., 1998.  Optimal scaling of discrete approximations to Langevin diffusions, J. Roy. Stat. Soc. B, 60(1), 255–268.
  * Roberts & Tweedie (1996) Roberts, G. O. & Tweedie, R. L., 1996.  Exponential Convergence of Langevin Distributions and Their Discrete Approximations, Bernoulli, 2(4), 341–363.
  * Shapiro & Ritzwoller (2002) Shapiro, N. M. & Ritzwoller, M. H., 2002.  Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle, Geophys. J. Int, 151(1), 88–105.
  * Skilling (2006) Skilling, J., 2006.  Nested sampling for general bayesian computation, Bayesian Analysis, 1(4), 833–859.
  * Starck et al. (2021) Starck, J. L., Themelis, K. E., Jeffrey, N., Peel, A., & Lanusse, F., 2021.  Weak-lensing mass reconstruction using sparsity and a gaussian random field, Astronomy & Astrophysics, 649.
  * Takahashi et al. (2017) Takahashi, R., Hamana, T., Shirasaki, M., Namikawa, T., Nishimichi, T., Osato, K., & Shiroyama, K., 2017.  Full-sky Gravitational Lensing Simulation for Large-area Galaxy Surveys and Cosmic Microwave Background Experiments, Astrophys. J., 850(1), 24.
  * Tarantola & Valette (1982) Tarantola, A. & Valette, B., 1982.  Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. Space Ge., 20(1), 219–232.
  * Trampert & Woodhouse (1995) Trampert, J. & Woodhouse, J. H., 1995.  Global phase velocity maps of Love and Rayleigh waves between 40 and 150 seconds, Geophys. J. Int, 122(2), 675–690.
  * Wallis et al. (2017) Wallis, C. G., Wiaux, Y., & Mcewen, J. D., 2017.  Sparse image reconstruction on the sphere: Analysis and synthesis, IEEE Trans. Image Process., 26(11), 5176–5187.
  * Wallis et al. (2021) Wallis, C. G. R., Price, M. A., McEwen, J. D., Kitching, T. D., Leistedt, B., & Plouviez, A., 2021.  Mapping dark matter on the celestial sphere with weak gravitational lensing, Mon. Not. R. Astron. Soc..
  * Wiaux et al. (2008) Wiaux, Y., McEwen, J. D., Vandergheynst, P., & Blanc, O., 2008.  Exact reconstruction with directional wavelets on the sphere, Mon. Not. R. Astron. Soc., 388(2), 770–788.
  * Wiaux et al. (2009) Wiaux, Y., Jacques, L., Puy, G., Scaife, A. M. M., & Vandergheynst, P., 2009.  Compressed sensing imaging techniques for radio interferometry, Mon. Not. R. Astron. Soc., 395(3), 1733–1742.
  * Woodhouse & Dziewonski (1984) Woodhouse, J. H. & Dziewonski, A. M., 1984.  Mapping the upper mantle: Three-dimensional modeling of earth structure by inversion of seismic waveforms, J. Geophys. Res.-Sol. Ea., 89(B7), 5953–5986.
  * Zhang & Curtis (2020) Zhang, X. & Curtis, A., 2020.  Variational full-waveform inversion, Geophys. J. Int, 222(1), 406–411.
  * Zhao et al. (2020) Zhao, X., Curtis, A., & Zhang, X., 2020.  Bayesian seismic tomography using normalizing flows, eartharxiv.
  * Zhou et al. (2005) Zhou, Y., Dahlen, F. A., Nolet, G., & Laske, G., 2005.  Finite-frequency effects in global surface-wave tomography, Geophys. J. Int, 163(3), 1087–1111.