{'fraction_non_alphanumeric': 0.03614155021614557, 'fraction_numerical': 0.012843068261787473, 'mean_word_length': 4.672797262617622, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse image reconstruction on the sphere: analysis and synthesis

We develop techniques to solve ill-posed inverse problems on the sphere by sparse regularisation, exploiting sparsity in both axisymmetric and directional scale-discretised wavelet space. Denoising, inpainting, and deconvolution problems, and combinations thereof, are considered as examples. Inverse problems are solved in both the analysis and synthesis settings, with a number of different sampling schemes. The most effective approach is that with the most restricted solution-space, which depends on the interplay between the adopted sampling scheme, the selection of the analysis/synthesis problem, and any weighting of the l1 norm appearing in the regularisation problem. More efficient sampling schemes on the sphere improve reconstruction fidelity by restricting the solution-space and also by improving sparsity in wavelet space. We apply the technique to denoise Planck 353 GHz observations, improving the ability to extract the structure of Galactic dust emission, which is important for studying Galactic magnetism.


I. INTRODUCTION
S PHERICAL images arise in many fields, from cosmology (e.g. [1]) to biomedical imaging (e.g. [2]), where inverse problems are often encountered. Sparse priors have proved highly effective in regularising Euclidean inverse problems, where sparsity is imposed in a wavelet space or sparsifying dictionary. In the spherical setting, wavelet theory is only recently starting to approach maturity, while a mature, general, and robust framework for sparse regularisation is lacking.
Over the last couple of decades there have been many developments regarding wavelet theory in spherical settings. Many initial attempts to extend wavelet transforms to the sphere differed primarily in the manner in which dilations are defined on the sphere [3]- [13]. These constructions were essentially based on continuous methodologies, which, although insightful, limited practical application to problems where the exact synthesis of a function from its wavelet coefficients is not required. A number of early discrete constructions followed [14]- [18]; however, many of these constructions do not necessarily lead to stable bases [19]. More recently, a number of discrete wavelet frameworks have emerged that have found considerable application, particularly in cosmology (e.g. [20]- [23]), including needlets [24]- [26]; scale-discretised wavelets [27]- [33]; and the isotropic undecimated and pyramidal wavelet transforms [34]. All three of these approaches have also been extended to analyse signals defined on the threedimensional ball formed by augmenting the sphere with the radial line [35]- [39], such as the galaxy distribution.
Solving Euclidean inverse problems by imposing sparse regularising priors has become increasingly popular in recent years. This trend has been driven by improving theoretical foundations for the recovery of sparse signals, facilitated by the theory of compressive sensing [40]- [42], and empirical results that have demonstrated the effectiveness of sparse priors for wide classes of natural images. Sparse reconstruction problems can be posed in either the synthesis or analysis settings [43]. In the synthesis setting, the sparse (e.g. wavelet) coefficients of the signal are recovered, from which the signal is synthesised. In the analysis setting, although sparsity is imposed in some sparsifying (e.g. wavelet) dictionary, the signal is recovered directly. When the dictionary considered is not an orthonormal basis but a redundant dictionary, the synthesis and analysis approaches exhibit quite different properties since the solution-space of the analysis problem is more restrictive than the synthesis problem [43]- [45]. Empirical studies have shown promising results for the analysis setting (e.g. [43], [46], [47]), which is hypothesised to be due to its more restrictive solution-space. Some progress has been made towards solving sparse regularisation problems on the sphere (e.g. [48]- [50]). Compressive sensing for signals sparse in spherical harmonic space is considered in [49], while inpainting problems are considered in [48], [50], imposing sparsity in a redundant dictionary [48] and the signal gradient [50].
In this work we consider general linear inverse problems on the sphere, including denoising, inpainting, and deconvolution problems, and combinations thereof, and apply sparse regularising priors in scale-discretised wavelet space, using both axisymmetric and directional wavelets [27]- [30]. Moreover, for the first time we study in detail the properties and empirical performance of the analysis and synthesis problems on the sphere. Furthermore, we investigate the impact of sampling theorems and schemes on the sphere [51]- [53] for sparse image recovery, which have already been shown to play a significant role [50]. In particular, we study the impact of the efficiency of sampling schemes in both the synthesis and analysis settings.
While sparse regularisation is often very effective, we close this introduction by cautioning against the blind application of sparse priors. For example, for the cosmic microwave background (CMB), which is to very good approximation a realisation of a Gaussian random field on the sphere, we recall that inpainting by imposing sparsity in spherical harmonic space (via the 1 norm) has the undesirable property of breaking statistical isotropy [54]. One must therefore take care in applying priors appropriate for the problem at hand.
The remainder of the article is structured as follows. In Sec. II we review harmonic analysis on the sphere S 2 and rotation group SO(3) and associated sampling theorems and schemes. In Sec. III we present the general framework to solve inverse problems on the sphere using sparse reconstruction. In Sec. IV we study sparse image reconstruction on the sphere through numerical experiments, comparing the analysis and synthesis settings and evaluating the impact of the sampling scheme used. In Sec. V we apply the methods proposed to denoise the Planck 353 GHz data. In Sec. VI we conclude.

II. SAMPLING ON THE SPHERE AND ROTATION GROUP
In this section we review the representation of signals on the sphere and the rotation group, in both the spatial and harmonic domains. We consider discretised signals, sampled according to different sampling schemes and sampling theories, that differ in the number of samples required to capture all of the information content of signals.

A. Signals on the Sphere
We consider the space of square integrable functions defined on the sphere S 2 . The canonical basis for the space of square integrable functions on the sphere is given by the spherical harmonics Y m ∈ L 2 (S 2 ), with natural ∈ N, integer m ∈ Z and |m| ≤ . Due to the orthogonality and completeness of the spherical harmonics, any square integrable function on the sphere x ∈ L 2 (S 2 ) may be represented by its spherical harmonic expansion where the spherical harmonic coefficients are given by the usual projection onto each basis function: where dΩ(ω) = sin θ dθ dϕ is the usual invariant measure on the sphere and ω = (θ, ϕ) denote spherical coordinates with colatitude θ ∈ [0, π] and longitude ϕ ∈ [0, 2π). Complex conjugation is denoted by the superscript * . Throughout, we consider signals on the sphere band-limited at L, that is signals such that x m = 0, ∀ ≥ L, in which case the summation over in Eq. (1) may be truncated to the first L terms.
In the discrete setting we can write the forward and inverse spherical harmonic transforms as linear operators, respectively: whereỸ ∈ C L 2 ×N S 2 and Y ∈ C N S 2 ×L 2 , with N S 2 denoting the number of samples on the sphere required to capture all the information content of a signal band limited at L. We denote the concatenated vector of N S 2 spatial measurements by x ∈ C N S 2 and the concatenated vector of L 2 harmonic coefficients byx ∈ C L 2 . Here and throughout we denote the forward harmonic transform with a tilde and the inverse transform without. Since sampling theorems on the sphere do not reach optimal dimensionality, as discussed in more detail below, the operatorsỸ and Y are not necessarily inverses of one another, e.g. YỸ = I (although we noteỸY = I, where I is the identity).
When considering images on the sphere the sampling theorem adopted can be of great significance. A sampling theorem allows one to transform from real space to harmonic space and back, without loss of information, from a finite number of samples N S 2 . Sampling theorems on the sphere differ in the number of samples N S 2 required. No existing sampling theorem on the sphere achieves the optimal number of samples of L 2 suggested by the harmonic dimensionality of a bandlimited signal. The canonical Driscoll & Healy [51] sampling theorem on the sphere (hereafter DH) requires ∼ 4L 2 samples to capture the information content of a signal band-limited at L. Recently, McEwen and Wiaux [52] (hereafter MW) developed a novel sampling theorem requiring ∼ 2L 2 samples only, thereby reducing the spherical Nyquist sampling rate by a factor of two. More recently, Khalid, Kennedy and McEwen [53] developed a new sampling scheme (hereafter KKM) that achieves the optimal number of L 2 samples. However, this scheme does not lead to a sampling theorem with theoretically exact spherical harmonic transforms; nevertheless, good numerical accuracy is achieved in practice.
Fast algorithms to compute spherical harmonic transforms, which avoid any pre-computation 1 , have been developed for the DH and MW sampling theorems, which scale as O(L 3 ) [51], [52], [55]. The complexity of the fast algorithm for the KKM sampling schemes scales as O(L 4 ), which can be reduced by appealing to algorithms to perform fast matrixvector multiplications, thereby reaching close to O(L 3 ) in practice [53].
Alternative sampling schemes also exist (e.g. HEALPix [56], IGLOO [57], GLESP [58]), although these are typically oversampled and the accuracy of numerical quadrature can in some cases be limited (e.g. HEALPix). Finite rate of innovation schemes to recover signals on the sphere comprised of a number of Dirac delta functions have also been developed [59]- [61]. However, such super-resolution approaches cannot be applied to recover general signals on the sphere. In this article we focus on efficient sampling schemes for general band-limited signals, which are also highly accurate (with accuracy close to machine precision): namely, the KKM [53], MW [52] and DH [51] schemes.

B. Signals on the Rotation Group
When considering directional wavelets it is necessary to be able to decompose and reconstruct square integrable signals defined on the rotation group SO(3), the space of threedimensional rotations, where rotations are parameterised by the Euler angles ρ = (φ, θ, ψ), with φ ∈ [0, 2π), θ ∈ [0, π] and ψ ∈ [0, 2π). We adopt the zyz Euler convention corresponding to the rotation of a physical body in a fixed coordinate system about the z, y and z axes by ψ, θ and φ, respectively.
The Wigner D-functions D mn ∈ L 2 (SO(3)), with natural ∈ N and integer m, n ∈ Z, |m|, |n| ≤ , are the matrix elements of the irreducible unitary representation of the rotation group SO(3) [62]. Consequently, the D * mn also form an orthogonal basis in L 2 (SO(3)). 2 Due to the orthogonality and completeness of the Wigner D-functions, any square integrable function on the rotation group x ∈ L 2 (SO(3)) may be represented by its Wigner expansion where the Wigner coefficients are given by the projection onto each basis function: where d (ρ) = sin θ dθ dφ dψ is the usual invariant measure on the rotation group. Note that ·, · is used to denote inner products over both the sphere and the rotation group (the case adopted can be inferred from the context). Throughout, we consider signals on the rotation group band-limited at L, that is signals such that x mn = 0, ∀ ≥ L, in which case the summation over in Eq. (5) may be truncated to the first L terms.
In the discrete setting we can write the forward and inverse Wigner transforms as linear operators, respectively: whereD denoting the number of samples on the rotation group required to capture all the information content of a signal band-limited at L. The harmonic dimensionality of a band-limited signal on the rotation group reads L(4L 2 − 1)/3. We denote the concatenated vector of N SO(3) spatial measurements by x ∈ C N SO(3) and the concatenated vector of harmonic coefficients byx ∈ C L(4L 2 −1)/3 . Again, we denote the forward harmonic transform with a tilde, and the inverse transform without, and note that the operatorsD and D are not necessarily inverses of one another.
For signals with limited directional sensitivity, it is convenient to consider a directional band-limit N , such that x mn = 0, ∀n ≥ min(N, ). In settings where the adopted band-limit (i.e. resolution) is not clear from the context, we adornD and D with superscripts denoting the spherical and directional band-limits adopted, e.g.D L,N . Sampling theorems on the rotation group can be constructed from a straightforward extension of sampling theorems defined on the sphere. Kostelec et al. [63] extend the DH sampling theorem, leading to a sampling theorem on the rotation group requiring ∼ 8L 3 samples. McEwen et al. [64] extend the MW sampling theorem, leading to sampling theorem requiring ∼ 4L 3 samples. The KKM sampling scheme has not yet been extended to the rotation group. 3 No sampling theorem on the rotation group reaches the optimal harmonic dimensionality of ∼ 4L 3 /3. We adopt the sampling theorem on the rotation group of McEwen et al. [64] in this work, where fast algorithms are developed to compute forward and inverse Wigner transforms that scale as O(N L 3 ).

III. SPARSE IMAGE RECONSTRUCTION ON THE SPHERE
We develop the proposed framework to solve inverse problems on the sphere in this section. We begin by reviewing wavelet transforms on the sphere, before presenting a discrete, operator formulation that illuminates the adjoint operators of the wavelet transform. Sparse regularisation problems on the sphere are then posed, in both analysis and synthesis settings, before the properties of these problems are discussed, along with algorithmic details for solving the problems, which require fast adjoint operators.

A. Wavelet Analysis and Synthesis
We adopt the scale-discretised wavelet transform on the sphere [27], [28], [30], [33], which supports directional wavelets. The wavelet transform is given by the directional convolution of each wavelet, Ψ j ∈ L 2 (S 2 ), with the signal of interest, x ∈ L 2 (S 2 ): where wavelet coefficients w j ∈ L 2 (SO (3)) incorporate directional information and so are defined on the rotation group. R ρ is the rotation operator that rotates by the Euler angles ρ. Wavelets are considered for a range of scales j, which runs from J min to J max . For further details on the wavelet construction and transform see, e.g., [30], [33]. The lowest frequency content of the signal of interest is extracted by the axisymmetric convolution of a scaling function, Υ ∈ L 2 (S 2 ), with the function of interest: where scaling coefficients s ∈ L 2 (S 2 ) are defined on the sphere since low-frequency directional structure is not typically of interest. In harmonic space these directional and axisymmetric convolutions read, respectively, (see, e.g., [30] The original image on the sphere can be reconstructed from its wavelet and scaling coefficients by (16) provided the wavelets and scaling function satisfy an admissibility criterion [27], [28], [30], [33]. In harmonic space, reconstruction reads In practice, for directional wavelet transforms we consider wavelets with an azimuthal band-limit N , i.e. Ψ j n = 0, ∀n ≥ min(N, ), which implies the directional wavelet coefficients w j also exhibit a directional band-limit N . Furthermore, directional wavelets with even or odd azimuthal symmetry are typically considered, in which case only N (rather than 2N −1) directions are required [27], [30], [33].

B. Discrete Wavelet Analysis and Synthesis
We formulate discrete, operator representations of the forward and inverse wavelet transforms that permit a clear construction of the adjoint wavelet operators. We consider the harmonic representation of the wavelet transform, which is inherently discretised, where we concatenate the harmonic coefficients into a single vector. The wavelet transform can be represented by its action on harmonic coefficients, followed by inverse harmonic transforms. A similar represententation is formulated for the inverse wavelet transform. By formulating wavelet transforms as a concatenation of operators, it is straightforward to construct operators representing adjoint wavelet transforms, which are required for solving sparse regularisation problems.
The harmonic expressions for the wavelet transform given by Eq. (14) and Eq. (15) may be written in terms of linear operators:ŵ whereŵ j denotes Wigner coefficients of the wavelet coefficients w j andŝ denotes spherical harmonic coefficients of the scaling coefficients s. The operators W j ∈ C N j (L j ) 2 ×L 2 and S ∈ C L 2 s ×L 2 implement harmonic space multiplication by the wavelet Ψ j * n and scaling function Υ * 0 , respectively, as described by Eq. (14) and Eq. (15), where the normalisation factor is not included in the former but is included in the latter (L j and N j are defined in detail below). The normalisation for the wavelets, given by 8π 2 /(2 +1), is applied by the operator . We separate out the normalisation in this case as it does not apply in the reconstruction of the signal seen in Eq. (17).
Harmonic space representations of wavelet and scaling coefficients are represented at the minimum resolution required to capture all signal content. Thus, the band-limit for each wavelet scale j is limited to L j and for the scaling function to L s . Wavelet Ψ j has support in the range [λ j−1 , λ j+1 ], where λ ∈ R is a scaling parameter that defines the scale dependance of each wavelet (for a standard dyadic scaling λ = 2), while the scaling function Υ has support in the range < λ Jmin (see [27], [28], [30], [33] for further details). Consequently, L j = λ j+1 and L s = λ Jmin . The azimuthal band limit of a wavelet scale is limited by the overall azimuthal band limit or the band limit of that scale, therefore N j = min(N, L j ).
We collect the harmonic representation of all wavelet and scaling coefficients in a single vector: where · † denotes the Hermitian transpose or adjoint, N = diag(I Ls , N Jmin , N Jmin+1 , . . . , N Jmax ), and W = diag(S, W Jmin , W Jmin+1 , . . . , W Jmax ). The collection of scaling and wavelet coefficients can be calculated from their harmonic representations by a series of inverse spherical harmonic and Wigner transforms by where H = diag(Y, D Jmin , D Jmin+1 , . . . , D Jmax ).
The forward wavelet transform, denoted by the operator Ψ, can then be expressed by the concatenation of operators defined above, yielding In other words, the wavelet transformΨ is composed of a spherical harmonic transformỸ, harmonic wavelet multiplication W, harmonic normalisation N, and inverse spherical harmonic and Wigner transforms H. The inverse wavelet transform, denoted by the operator Ψ can be represented in a similar manner. Firstly, we note the spherical harmonic and Wigner coefficients of the wavelet and scaling coefficients can be calculated by a series of forward harmonic transforms byα =Hα, where the final equality follows by noting (NW) −1 = W † , which can be inferred from Eq. (17), which in turn follows by the wavelet admissibility criterion. In other words, the inverse wavelet transform Ψ is composed of forward spherical harmonic and Wigner transformsH, harmonic wavelet multiplication and summation W † , and an inverse spherical harmonic transform Y.

C. Adjoints of Discrete Wavelet Analysis and Synthesis
The operator descriptions of the forward and inverse wavelet transforms formulated in the previous subsection permit a straightforward construction of the adjoint operators, which are required to solve inverse problems on the sphere when imposing sparsity in wavelet space. From Eq. (25), the wavelet transform operator readsΨ = HNWỸ, with corresponding adjoint: where we have noted that N is self-adjoint, i.e. N † = N. From Eq. (27), the inverse wavelet transform operator reads Ψ = YW †H , with corresponding adjoint: In the discrete setting, recall that the adjoint and inverse op- Consequently, for practical application, fast algorithms must be developed to compute not only forward transforms and their inverses but also the adjoints of both the forward and inverse transforms. This has been performed already [50] for the spherical harmonic transforms of the MW sampling theorem [52], i.e. to applyỸ † and Y † . To apply directional wavelets, we also require fast algorithms to compute adjoint Wigner transforms. In Appendix A, we develop fast algorithms to computeD † and D † for the Wigner transforms corresponding to the sampling theorem on the rotation group of [64].

D. Sparse Regularisation
We consider linear, ill-posed inverse problems defined on the sphere, including, for example, denoising, inpainting, and deconvolution problems. Consider M measurements y ∈ R M of the signal on the sphere x ∈ R N S 2 , acquired according to the measurement equation where Φ ∈ R M ×N S 2 is the measurement operator and n ∈ R N S 2 is measurement noise, assumed to be Gaussian, i.e. n ∼ N (0, σ), where σ = x 2 × 10 (−SNR/20) , where SNR is in dB. For example, the measurement operator Φ may model the beam or point-spread function of a sensor (in a deconvolution problem) or a masking of the signal (in an inpainting problem). We regularise the spherical ill-posed inverse problem of Eq. (30) by promoting sparsity in wavelet space by posing synthesis and analysis problems on the sphere. The synthesis problem reads where the signal is then recovered from its wavelet coefficients by x = Ψα . The analysis problems reads where we recover the signal x directly. 4 The 1 norm appearing in the sparsity constraint must be defined appropriately for the spherical setting [50], as discussed in more detail below. The square of the residual noise follows a scaled χ 2 distribution with M degrees of freedom, i.e. y − Ψx 2 2 ∼ σ 2 χ 2 (M ). Consequently, we choose to correspond to a given percentile of the χ 2 distribution [50].
When solving the synthesis and analysis problems of Eq. (31) and Eq. (32) we are free to choose different sampling schemes (e.g. KKM, MW, DH sampling). In Euclidean space, the analysis problem has shown promising results in empirical studies, which we recall is hypothesised to be due to the more restrictive solution-space of the analysis setting. This relationship between the size of the solution-space and the analysis and synthesis settings does not in general carry over to the spherical setting since on the sphere sampling is not typically optimal. Consequently, recovering the signal directly in the analysis setting does not necessarily lead to the most restrictive solution-space. The most restrictive solutionspace depends on the interplay between the adopted sampling scheme, the selection of the analysis/synthesis problem, and any weighting of the 1 norm, which is made explicit in the following subsection. 4 Our framework can be applied to spin s ∈ Z signals on the sphere (see e.g. [52]) in a straightforward manner, noting that the spin wavelet transform of [30] can be represented by the operators sΨ = sHNWsỸ and Ψ = sYW † sH, where the forward and inverse scalar spherical harmonic transforms are replaced by spin versions, e.g. replacing Y by sY.

E. Algorithmic Details
The 1 norm appearing in Eq. (31) and Eq. (32) must be defined appropriately for the spherical setting, taking into account the sampling scheme adopted. In [50], where the total variation (TV) norm is considered, the associated discrete TV norm is weighted by the exact quadrature weights of the sampling theorem adopted in order to approximate the continuous norm. Through numerical experiments we have found the 1 norm, and the solution of the sparse reconstruction problems, to be relatively insensitive to the exact form of weights: provided weights capture the area of each pixel the underlying continuous norm is well approximated and it is not necessary to use exact quadrature weights. 5 Consequently, for all sampling schemes we adopt the following weights for the wavelet and scaling coefficients, respectively, corresponding to scale j and pixel p: where n j φ , n j θ and n j ψ are the number of samples in the φ, θ and ψ directions, θ p is the θ coordinate of the sample p, and η ∈ R + is a decay parameter.
The weights approximate the area of each pixel, normalised by the energy of the wavelet and scaling function for the given scale j. Furthermore, for the weighting of wavelet coefficients the additional factor (λ j ) η is introduced. The term λ j corresponds to the middle harmonic multipole to which the wavelet Ψ j is sensitive, while η is introduced as a free parameter to incorporate prior knowledge of natural signals, i.e. to control the wavelet decay imposed as a prior when solving the sparse regularisation problems. Increasing η promotes large scale features by increasing the weight applied to small wavelet scales, thereby increasing their penalty to the 1 norm. Moreover, for the synthesis setting, increasing η reduces the effective size of the solution-space.
We use the Douglas-Rachford (DR) [67] splitting algorithm to solve the sparse regularisation problems posed in Sec. III-D. The Douglas-Rachford algorithm [67] is based on a splitting approach that requires the computation of two proximity operators. In our case, one proximity operator is based on the 1 norm and the other on the data fidelity constraint. The adaptation of the DR algorithm to the sphere is discussed further in [50]. The DR algorithm requires the adjoint of the operators that appear in the problem specification, e.g. the adjoint sparsifying operators shown in Eq. (28) and Eq. (29). In numerical experiments, if inverse operators are used in place of the adjoints, we have seen convergence failures. In [50], fast adjoints for the spherical harmonic transform corresponding to the MW sampling scheme [52] were dervied. In an analogous manner, we derive in Appendix A fast adjoints for the Wigner transforms of [64]. These efficient adjoints have a numeric complexity of order N L 3 compared to the naive adjoint operations which have complexity of order N 2 L 4 . The power method is used to calculate the norms of the operators (required when solving the optimisation problems).

IV. NUMERICAL EXPERIMENTS
We perform numerical experiments to both assess the effectiveness of imposing sparsity in wavelet space and to test the impact of the sampling scheme used and whether or not the problem is solved in the analysis or synthesis setting.
We compare the analysis and synthesis settings for the KKM, DH and MW sampling schemes. These experiments are performed at low resolution as fast adjoint transforms for the DH and KKM sampling theorems are lacking. We chose to solve a noisy inpainting problem for these tests (as an example of a common inverse problem). At high resolution we demonstrate image reconstruction with both axisymmetric and directional wavelet sparsity priors using MW sampling (for which we have constructed fast adjoint operators). We test the method at high resolution on inpainting and deconvolution problems, and a combined inpainting and deconvolution problem, all in the presence of noise.
We generate low and high resolution test images from Earth topography data. The original Earth topography data are taken from the Earth Gravitational Model (EGM2008) publicly released by the U.S. National Geospatial-Intelligence Agency (NGA) EGM Development Team. 6 The ground truth for the low and high resolution experiments performed in the remainder of this section are shown in Fig. 1.
Much of the work in this section takes advantage of publicly available codes: we use SSHT 7 [52] and NSHT 8 [53] to compute spherical harmonic transforms; SO3 9 [64] to compute harmonic transforms on the rotation group; S2LET 10 [28], [30] to compute wavelet transforms on the sphere; and SOPT 11 [47] to solve inverse problems.
The main code used to perform these experiments was not optimised and significant parts are run in serial in the high level language MATLAB. The runtimes for this naive implementation were as follows: the low resolution experiments described in Section IV-A take approximately 30 seconds for the MW sampling, 1 minute for the DH sampling, and 1 hour for the KKM sampling. The high resolution experiments performed with the MW sampling described in Section IV-B take around 10 minutes. All experiments were performed on a MacBook Pro (early 2015), with a 2.9 GHz Intel Core i5 processor and 16 GB of RAM.

A. Low Resolution Axisymmetric Experiments
We first solve a simple inpainting problem at low resolution (L = 32) using axisymmetric wavelets (N = 1). In this case and represents a uniformly random masking of the spherical image, with one non-zero, unit value on each row specifying the location of the measured datum. The adjoint of the operator can be calculated trivially. Measurements are taken according to Eq. (30) with noise included corresponding to a signal-tonoise-ratio (SNR) of 46 dB, where SNR = 20 log( x 2 / x * − x 2 ), defined in harmonic space to avoid differences due to the number of samples of each sampling scheme. We vary the number of measurements taken as M = N m L 2 , where N m = [0.3, 0.5, 1.0, 1.5, 1.9]. We run these experiments for each of the sampling schemes we consider, specifically the KKM, MW and DH sampling schemes. Results can be seen in Fig. 2. The improvement given by the lower number of samples of the KKM or MW sampling schemes can be seen clearly. There is also typically an improvement when solving Eq. (31), the synthesis problem, as opposed to Eq. (32), the analysis problem. We set η = 2.5 (as we do for the remainder of the article unless otherwise stated), since this was shown to be the value resulting in the reconstructions of the highest SNR, although it should be noted that the resulting SNR was very similar for 2.5 < η < 4.5 when an emperical investigation was conducted. In Fig. 3 we show the average SNR for 10 reconstructions, which supports the findings inferred from Fig. 2. These experiments were repeated on another data set from a different domain (natural light probe images) with similar results.

B. High Resolution Experiments
We run three example high resolution (L = 128) experiments on the high resolution Earth map shown in Fig. 1. We solve all the problems with the MW sampling in the synthesis setting with η = 2.5. We consider the MW sampling as there are currently no fast adjoint algorithms for the other sampling theorems and for the synthesis setting as it was shown to be superior in Sec. IV-A We consider both axisymmetric wavelets (N = 1) and directional wavelets with N = 4 which leads to wavelets with odd azimuthal symmetry. Firstly, we solve a simple inpainting problem with a measurement operator given by Eq. (35), with N m = 1.0. Secondly, we consider a deconvolution problem, where the measurement operator is where G ∈ R L 2 ×L 2 is a diagonal matrix whose elements are, where σ = π/L. Thirdly we solve the combined problem, We consider measurement noise with SNR of 46 dB. We show the results of these experiments and a band limited version of the measured data in Fig. 4. The SNR of the recovered images are encouragingly high, showing a good similarity with the ground truth. The deconvolution problem shows minor visual artefacts in the axisymmetric wavelet case. The visual artefacts are reduced in the directional setting for the two problems involving deconvolution, however SNRs are slightly lower. The inpainting problem visually shows a marked improvement in the directional case over the axisymmetric case, as also illustrated by the improved SNR.

V. DENOISING Planck 353 GHZ DATA
The Planck satellite observed the entire sky at a range of microwave frequencies [1], yielding high resolution maps of the polarised Galactic dust emission from its high frequency polarised channel centred on 353 GHz and total intensity maps at even higher frequencies from other channels [68]. One of the many uses of these maps is the study of the Galactic magnetic field, where it is important to have high SNR maps of the clumps of dust in the Galaxy. It is common practise to smooth the data with a Gaussian kernel in order to suppress the high frequency noise [69]. This smoothing has the undesirable effect of not denoising large scales and, perhaps more damaging, removing important structure on small scales. Here we examine the use of sparsity in wavelet space as a prior to denoise the Planck 353GHz total intensity map.
The Planck 353 GHz data is available to download 12 in HEALPix 13 format [56]. We use the HEALPix software to compute the spherical harmonic coefficients of this spherical image. We then band limit these to L = 2048 and use SSHT [52] to obtain a MW sampled image of the sphere. This is taken as our input data to then be denoised and can be seen in Fig. 5. An estimate of the noise level is made by downloading the noise maps from the same archive, performing the same operation, and averaging the noise over all of the sky. This leads to an initial estimate of , which is subsequently optimised through experimentation.
We solve the denoising problem in the synthesis setting with measurement operator set to the identity and using the axisymmetric wavelets (N = 1). We set η = 3.0 but have found the specific value have little effect on the reconstruction. Fig. 6 shows the original map, the result from denoising and a smoothed map. The smoothed map is the original map convolved with a 5 arcmin kernel to replicate the current denoising techniques adopted. It is clear the noise is reduced by our sparse denoising approach, while preserving small scale structure.

VI. CONCLUSIONS
We develop a general framework to solve image reconstruction problems on the sphere by sparse regularisation, minimising the 1 norm of wavelet coefficient representations of spherical images. By developing fast adjoint operators, we recover convergence guarantees for the resulting convex optimisation problems. As examples, we have demonstrated that using our framework one can solve denoising, inpainting, and deconvolution problems effectively, and combinations thereof.
We study and compare the analysis and synthesis settings for solving inverse problems on the sphere for the first time. The analysis problem has shown promising results in Euclidean space, hypothesised to be due to its more restrictive nature. However, the more restrictive nature of the analysis framework in Euclidean space does not carry over to the spherical setting. The most restrictive solution-space on the sphere depends on the interplay between the adopted sampling scheme, the selection of the analysis/synthesis problem, and any weighting of the 1 norm. We examine a variety of sampling schemes on the sphere, including the DH [51] and

Acquired Data
Recovered Image Axisymmetric Recovered Image Directional  MW [52] sampling theorems (leading to theoretically exact spherical harmonic transforms) and the KKM [53] sampling scheme (leading to approximate but highly accurate spherical harmonic transforms). DH, MW and KKM sampling requires ∼ 4L 2 , ∼ 2L 2 , and L 2 samples, respectively.
To examine the analysis and synthesis problems and the impact of the various sampling schemes considered, we study results from a simple inpainting problem at low resolution. In the numerical results shown in Fig. 2 and Fig. 3 we find that the synthesis setting typically out-performs the analysis setting for suboptimal sampling schemes. Moreover, reconstruction fidelity is enhanced further by adopting more efficient sampling schemes that require fewer samples to capture the information content of signals on the sphere. These findings are robust to the choice of test signal.
We hypothesise that differences between the analysis and synthesis settings on the sphere are due to restrictions of the solution-space. Elad et al. [43] showed theoretically and with simulations that there are fundamental differences between the analysis and synthesis settings due to the solution-spaces. As in Euclidean space, we find settings with a more restricted solution-space yield superior performance. In contrast to Euclidean space, it is the synthesis setting rather than the analysis setting that typically results in a more restrictive solution-space on the sphere. This is due to the inefficiency of spherical sampling schemes and the weighting introduced in the 1 norm. Since the KKM sampling exhibits the optimal number of samples on the sphere there is no appreciable different between the size of the solution-spaces of the analysis and synthesis problems and the fidelity of spherical images recovered by the analysis and synthesis approaches are similar. In addition, the sparsity of band-limited signals is further promoted by more efficient sampling schemes when considering a sparse representation that captures spatial localisation [50], such as wavelets.
We also demonstrate solving inverse problems in a number of high resolution settings, facilitated by our fast adjoint operators. We solve inpainting, deconvolution, and combined inpainting and deconvolution problems, all in the presence of noise, using both axisymmetric and directional wavelets.

Acquired Data
Recovered Image Smoothed Image Fig. 6. Results from denoising the Planck 353 GHz total intensity map. The left column shows the acquired data, the middle column shows the denoised data and the right column shows the data denoised by smoothing with a 5 arcmin Gaussian kernel (the standard approach [69]). Each plot shows a zoomed region of the sphere. The acquired data shown here is a zoom in of that shown in Fig. 5, the top figure is the blue region and the bottom figure is the green region.
For all inverse problems considered our sparse regularisation techniques yield excellent reconstruction fidelity. Our framework for solving inverse problems on the sphere can be applied to many real-world problems. We have shown that our framework can be used to effectively denoise 353 GHz channel observations from the Planck satellite, which will be useful for studying Galactic magnetic fields.

APPENDIX
Standard convex optimisation methods require not only the application of the operators that appear in the optimisation problem but also their adjoints. Moreover, these methods are typically iterative, necessitating repeated application of each operator and its adjoint. Thus, to solve optimisation problems that incorporate Wigner transform operators fast algorithms to apply both the operator and its adjoint are required to render high-resolution problems computationally feasible.
Here we develop fast algorithms to perform adjoint forward and adjoint inverse Wigner transforms for the extension of the MW sampling scheme to the rotation group [64]. The fast adjoint follows by taking the adjoint of each stage of the fast standard transforms [64] and applying these in reverse order. The forward and inverse transforms can be found in [64,Sec. 3]. Here we use notation consistent with that work, where ∆ m m = D mn (0, π/2, 0). The fast adjoint of the forward transform is as follows: Similarly, the fast adjoint of the inverse transform is as follows: