A novel sampling theorem on the sphere
Abstract
We develop a novel sampling theorem on the sphere and corresponding fast algorithms by associating the sphere with the torus through a periodic extension. The fundamental property of any sampling theorem is the number of samples required to represent a bandlimited signal. To represent exactly a signal on the sphere bandlimited at , all sampling theorems on the sphere require samples. However, our sampling theorem requires less than half the number of samples of other equiangular sampling theorems on the sphere and an asymptotically identical, but smaller, number of samples than the GaussLegendre sampling theorem. The complexity of our algorithms scale as , however, the continual use of fast Fourier transforms reduces the constant prefactor associated with the asymptotic scaling considerably, resulting in algorithms that are fast. Furthermore, we do not require any precomputation and our algorithms apply to both scalar and spin functions on the sphere without any change in computational complexity or computation time. We make our implementation of these algorithms available publicly and perform numerical experiments demonstrating their speed and accuracy up to very high bandlimits. Finally, we highlight the advantages of our sampling theorem in the context of potential applications, notably in the field of compressive sampling.
armonic analysis, sampling methods, spheres.
1 Introduction
\PARstartIn many fields of science and engineering data are measured on a spherical manifold. Applications where data are defined inherently on the sphere are found in computer graphics (e.g. [1]), planetary science (e.g. [2, 3, 4, 5]), geophysics (e.g. [6, 7, 8]), quantum chemistry (e.g. [9, 10]) and astrophysics (e.g. [11, 12]), to quote only a few. In many of these applications a harmonic analysis of the data is insightful. For example, spherical harmonic analyses have been remarkably successful in cosmology over the past decade, leading to the emergence of a standard cosmological model. Observations of the anisotropies of the cosmic microwave background (CMB), which are made on the celestial sphere, contain a wealth of information about the early Universe. Cosmologists extract this information from the angular power spectrum of observations of the CMB, computed through a harmonic transform on the sphere (e.g. [13]). Recent and upcoming fullsky observations of the CMB are of considerable size, containing approximately three [12] and fifty [14] million samples respectively. Furthermore, observations are made of both the temperature and polarisation of the CMB, which give rise to scalar and spin functions on the sphere respectively. Consequently, the ability to perform fast scalar and spin spherical harmonic transforms on large data sets is of considerable importance in cosmology and beyond.
Algorithms to perform spherical harmonic transforms have received considerable attention already. Some correspond to sampling theorems on the sphere, where the forward and inverse transform are theoretically exact for a bandlimited signal on the sphere. Others adopt approximate quadrature rules on the sphere, resulting in approximate harmonic transforms that do not correspond to sampling theorems on the sphere. However, these approximate algorithms typically arise from particular pixelisations of the sphere which meet desirable practical criteria, such as pixels of equal area, and are no less important. We focus here on approaches that lead to sampling theorems on the sphere with theoretically exact transforms for signals on the sphere bandlimited at (where the harmonic bandlimit is defined formally in Sec. 3.2). Note that the current [12] and forthcoming [14] CMB observations discussed previously support bandlimits of and respectively. All sampling theorems require samples on the sphere of order , however the exact number of samples required varies for each sampling theorem. For many applications reducing the number of samples required to represent a bandlimited signal on the sphere is of fundamental importance.
Sampling theorems on the sphere and their associated numerical algorithms are evaluated by four criteria: (i) the number of samples required to represent a bandlimited signal exactly; (ii) their computational complexity; (iii) their speed; and (iv) issues surrounding any precomputation. From an information theoretic viewpoint, the fundamental property of any sampling theorem is the number of samples required to represent a bandlimited signal exactly. In this article we review algorithms to compute spherical harmonic transforms accurately and efficiently. We also present a novel sampling theorem on the sphere and corresponding fast algorithms. Our approach compares favourably to the stateoftheart as evaluated by the four criteria listed previously. Furthermore, our algorithms apply to both scalar and spin functions on the sphere without any change in asymptotic complexity or computation time.
The remainder of this article is structured as follows. In Sec. 2 we review comprehensively the literature regarding the computation of spherical harmonic transforms, placing our novel sampling theorem and fast algorithms in the context of preceding work. Harmonic analysis on the sphere is reviewed concisely in Sec. 3, to present the mathematical preliminaries required subsequently. Our sampling theorem on the sphere and the corresponding fast algorithms to compute spherical harmonic transforms are derived in Sec. 4. In Sec. 5 we evaluate our algorithms numerically and discuss the advantages of our sampling theorem in the context of potential applications, notably in the field of compressive sampling. Concluding remarks are made in Sec. 6.
2 Review
The development of sampling theorems on the sphere and fast algorithms to compute spherical harmonic transforms has been driven largely by researchers in the fields of computational harmonic analysis, geophysics and astrophysics. Due to the diverse nature of these fields, the literature on this topic appears to be somewhat disjoint. We attempt to unify these works here and to present a comprehensive review of the historical development of the field.
For isolatitude sampling schemes, where the samples are gathered in isolatitude annuli, a separation of variables may be performed to rewrite the scalar spherical harmonic transform as a Fourier transform in longitude and an associated Legendre transform in colatitude. The computation is then dominated by the associated Legendre transform, reducing the complexity from to . Isolatitude samplings are thus very common and hence we restrict our attention to such schemes. Spin lowering and raising operators can be used to relate the harmonic transform of spin functions to the scalar case, hence scalar transforms have received the majority of attention. Approaches to improve the computational performance of scalar spherical harmonic transforms attempt either to reduce the asymptotic complexity further through fast associated Legendre transforms or to reduce the constant prefactor associated with the asymptotic scaling of the algorithm.
First attempts to compute a scalar spherical harmonic transform through a fast Legendre transform were performed by Orszag [15] and were based on a WentzelKramersBrillouin (WKB) approximation. Alternative approaches using the fast multipole method (FMM) [16] have been considered by Alpert & Rokhin [17] and by Suda & Takani [18]. The complexity of these algorithms scale linearly with the desired accuracy. An alternative approximate algorithm using WKB frequency estimates has been developed by Mohlenkamp [19, 20] for functions with bandlimits that are a power of two, however this approximation can be controlled independently of complexity, which scales as . In any case, these types of approach are necessarily approximate and do not yield exact sampling theorems on the sphere.
Other approximate approaches based solely on the separation of variables have been developed more recently for pixelisations of the sphere that satisfy certain practical requirements, such as HEALPix
Exact transforms with associated sampling theorems have been constructed for particular pixelisation schemes. It is wellknown that GaussLegendre quadrature may be used to construct exact spherical harmonic transforms. To our knowledge, this result was first highlighted in published work by Shukowsky [23], which in turn refers to unpublished (and inaccessible) work by Payne from 1971 [24]. An exact sampling theorem can be constructed from samples on the sphere, where the sample locations in colatitude are chosen as the roots of the Legendre polynomials of order , as dictated by GaussLegendre quadrature. Through a separation of variables, the resulting algorithm is . The GLESP
The first theoretically exact sampling theorem on an equiangular pixelisation was developed by Shukowsky [23], requiring samples on the sphere; however, the exactness of this approach was not studied numerically. Although this algorithm remains , a separation of variables may be used to reduce the computational complexity to . An alternative sampling theorem on the sphere for an equiangular pixelisation was developed by Driscoll and Healy [26].
Moreover, a divideandconquer approach to computing a fast associated Legendre transform in the cosine basis was derived [26]. The resulting algorithm is exact in exact precision arithmetic, and has computational complexity , but is known to suffer from stability problems [27, 28]. Healy et al. [27] readdressed the work of Driscoll & Healy [26], reformulating the sampling theorem on the sphere and developing some variants of the original algorithm,
Other approaches to reduce the cost of computing spherical harmonic transforms focus on reducing the constant prefactor associated with the asymptotic complexity of algorithms. These approaches have typically exploited fast Fourier transforms (FFTs) on equiangular pixelisations to reduce computation time through an association between the sphere and the torus, while their complexity remains . To our knowledge, the first algorithm based on this technique was developed by Dilts [33], where the North and South poles of the sphere were identified to map the sphere to the torus. However, this algorithm is approximate and does not result in a sampling theorem on the sphere. One of the authors of this article developed a sampling theorem on the sphere [34] by making periodic extensions of the sphere in colatitude in order to make an association with the torus. This sampling theorem requires samples on the sphere and, moreover, applies to any spin number without the application of spin lowering and raising operators. However, the forward algorithm associated with this sampling theorem proved unstable (the inverse algorithm did not suffer from stability issues). Recently, Huffenberger & Wandelt [35] adopted the inverse transform of this approach and resolved the instability in the forward algorithm, although in doing so increased the number of points required to sample a bandlimited function on the sphere to . In this article we readdress sampling theorems derived by associating the sphere with the torus through periodic extensions and develop a sampling theorem requiring samples on an equiangular pixelisation, with corresponding fast algorithms that do not suffer from any stability issues.
3 Harmonic Analysis on the Sphere
In this section we review harmonic analysis on the twosphere . We first review the scalar spherical harmonic transform, before generalising to the spin case. Associations are then made between the spin spherical harmonics and the Wigner functions, where the latter provide an orthogonal basis for the decomposition of square integrable functions on the rotation group .
3.1 Scalar spherical harmonics
We consider the space of square integrable functions on the sphere , with the inner product of defined by
where is the usual invariant measure on the sphere and define spherical coordinates with colatitude and longitude . Complex conjugation is denoted by the superscript .
The scalar spherical harmonic functions form the canonical orthogonal basis for the space of scalar functions on the sphere and are defined by
for natural and integer , , where are the associated Legendre functions. We adopt the CondonShortley phase convention, with the phase factor included in the definition of the associated Legendre functions, to ensure that the conjugate symmetry relation holds. The orthogonality and completeness relations for the spherical harmonics read and
respectively, where is the Kronecker delta symbol and is the Dirac delta function.
Due to the orthogonality and completeness of the scalar spherical harmonics, any square integrable scalar function on the sphere may be represented by its spherical harmonic expansion
where the spherical harmonic coefficients are given by the usual projection onto each basis function: . The conjugate symmetry relation of the spherical harmonic coefficients of a real function is given by , which follows directly from the conjugate symmetry of the scalar spherical harmonic functions.
3.2 Spin spherical harmonics
Square integrable spin functions on the sphere , with integer spin , are defined by their behaviour under local rotations. By definition, a spin function transforms as
(1) 
under a local rotation by , where the prime denotes the rotated function. It is important to note that the rotation considered here is not a global rotation on the sphere, such as that represented by an element of the rotation group , but rather a rotation by in the tangent plane at . The sign convention that we adopt here for the argument of the complex exponential in (1) differs to the original definition [36] but is identical to the convention used recently in the context of the polarisation of the CMB [37].
The spin spherical harmonics form an orthogonal basis for spin functions on the sphere for . Spin spherical harmonics were first developed by Newman & Penrose [36] and were soon realised by Goldberg [38] to be closely related to the Wigner functions. We therefore defer the explicit definition of the spin spherical harmonic functions until Sec. 3.3. The conjugate symmetry relation given for the spin spherical harmonics is given by . The spin spherical harmonics satisfy identical orthogonality and completeness relations as the scalar spherical harmonics.
Due to the orthogonality and completeness of the spin spherical harmonics, any square integrable spin function on the sphere may be represented by its spherical harmonic expansion
where the spin spherical harmonic coefficients are given by the usual projection onto each basis function: . Note that the spin spherical harmonics and transforms simply generalise the scalar equivalents to spin signals, reducing to the standard scalar case for . When deriving our novel sampling theorem we consider signals on the sphere bandlimited at , that is signals such that , . The conjugate symmetry relation of the spin spherical harmonic coefficients is given by for a function satisfying (which for a spin function equates to the usual reality condition) and follows directly from the conjugate symmetry of the spin spherical harmonics.
3.3 Wigner functions
The Wigner functions , for natural and integer , form an orthogonal basis for the space of square integrable functions on the rotation group,
and are parameterised by the Euler angles , where , and .
(2) 
where the real polar functions are defined by [39]
(3) 
where are the Jacobi polynomials. Note that recursion formulae are available to compute rapidly the Wigner functions (e.g. [40, 41]). The functions satisfy a number of symmetry relations; in this work we make use of the symmetry relations [39]
(4) 
(5) 
and
(6) 
We are not concerned with decompositions of functions on the rotation group in this article but rather representations of the spherical harmonics by Wigner functions. The spin spherical harmonics may be defined by the Wigner functions through [38]
(7) 
Defining the spherical harmonics in this manner allows us to apply standard Wigner function decompositions to the spherical harmonic functions. We subsequently make considerable use of the Fourier series decomposition of the functions given by [42]:
(8) 
where . This expression follows from a factoring of rotations as highlighted by Risbo [40]. The Fourier series representation of given by (8) allows one to write the spherical harmonic expansion of in terms of a Fourier series expansion of extended appropriately to the twotorus (as discussed in more detail in Sec. 4). Consequently, (8) is fundamental to the derivation of our sampling theorem on the sphere and fast algorithms.
4 Fast Spherical Harmonic Transform
We derive fast algorithms for performing forward and inverse spin spherical harmonic transforms and discuss the corresponding sampling theorem on the sphere. Our approach involves an extension of the sphere to the torus so that FFTs may be exploited to reduce the cost of computation. It is related closely to the algorithms derived by one of the authors in a previous work [34] and to the algorithms derived by Huffenberger & Wandelt [35], however it does not suffer from the instabilities of the former approach and requires only half as many samples on the sphere as the latter approach. We first present the general harmonic formulation of our algorithms, followed by a discussion of periodisation and discretisation, before deriving algorithms to perform exact forward and inverse transforms. Our sampling theorem also leads to a new quadrature rule on the sphere, which we then present, followed by a discussion of symmetries that may be exploited to improve the efficiency of computation for real signals.
4.1 Harmonic formulation
We consider the harmonic transform of spin functions on the sphere , bandlimited at ; consequently, all summations over or up to are truncated to . Furthermore, harmonic coefficients are not defined for , hence we define them to be zero to enforce the contraint when summations are interchanged.
By noting the definition of the spin spherical harmonics in terms of Wigner functions (7), the Wigner decomposition (2) and the Fourier expansion of the Wigner functions (8), the forward transform of may be written
(9) 
where
(10) 
and
(11) 
In Sec. 4.4 we consider implicit quadrature rules to evaluate (10) and (11) exactly. By noting the same substitutions and interchanging the order of summation, the inverse transform of may be written
(12) 
where
(13) 
and
(14) 
Although recasting the forward and inverse spherical harmonic transforms in this manner is no more efficient than the original formulation, expressions (10)–(13) highlight similarities with Fourier series representations. However, the Fourier series expansion is only defined for periodic functions; thus, to recast these expressions in a form amenable to the application of Fourier transforms we must make a periodic extension in colatitude .
4.2 Periodic extension
We make a periodic extension of to the domain so that FFTs may be used to compute the forward and inverse spherical harmonic transform rapidly. When making this periodic extension we must be careful to ensure that the symmetry of our current representation is respected on the new domain; we must apply the symmetries dictated by the inverse transform when imposing periodisation in the forward transform. By substituting (12) into (11) and noting the continuous orthogonality of the complex exponentials, we find the forward and inverse expressions in are related by
(15) 
Consequently, the symmetry we impose in when extended periodically must match the symmetry of . By reflecting , we obtain the following symmetry for :
where we have noted the symmetry
(16) 
following from (14) and (5). Thus, we extend to the domain by constructing
Note that we adopt a different periodic extension to other approaches framed on the torus [34, 35] by applying the extension to the Fourier transform of in , i.e. to . Two periodic extensions, one even and one odd, were required in the approach taken previously by one of the authors [34]. Huffenberger & Wandelt [35] apply the factor as a shift in by , removing the need for two periodic extensions but requiring a even number of samples in , precluding an association with the odd number of points in unless oversampling is performed. We avoid these restrictions by applying the periodic extension to the Fourier transform in of , rather than to directly.
4.3 Discretisation
We adopt an equiangular sampling of the sphere with sample positions given by
(17) 
and
(18) 
In order to extend the domain to we simply extend the domain of the index to include .
An odd number of sample points are required in both and in the extended domain so that a direct association may be made with the harmonic indices and , resulting in samples on the sphere. We also require a symmetric sampling in about the South pole so that samples on the extended domain can be obtained by reflecting samples defined on the original domain. The node positions specified by (17) and (18) eliminate repeated samples at the poles and since these points are excluded from the pixelisation. However, it is not possible to eliminate repeated samples at , since we require a discretisation that is symmetric about but which contains an odd number of sample points.
4.4 Forward transform
The algorithm we derive in this section to compute a forward spin spherical harmonic transform essentially follows the harmonic formulation presented in Sec. 4.1, however we discuss implicit quadrature rules for the exact evaluation of (10) and (11).
Since (11) is simply a Fourier transform we may appeal to the discrete and continuous orthogonality of the complex exponential to express this integral exactly by
for . An FFT may be used to compute for all and with computational complexity . We extend to the domain through the construction
noting .
We now consider an implicit quadrature rule for the exact evaluation of through (10). Firstly, however, we compute from by noting (15) and by appealing to the discrete orthogonality of the complex exponentials to invert (13), giving
An FFT may be used to compute for all and with computational complexity . Now that we have to hand, we substitute (13) into (10), noting (15), to yield
(19) 
where the weights are given by
Note that the definition of these weights is identical to that derived by Huffenberger & Wandelt [35], however we correct some (typographical) errors in their explicit evaluation. Since we are concerned with the values of for , the computation of through (19) explicitly requires weights with argument up to . If the range of is extended to , (19) may be seen as a (reflected) convolution, which may be computed more efficiently following a Fourier transform. However, since only the range is of interest, aliasing may be tolerated provided that it is outside of this range. To ensure that this is the case we zeropad in the domain prior to computing an inverse Fourier transform. We then compute an inverse Fourier transform of the weights on the same extended domain and take the product of these terms, the Fourier transform of which gives . Using FFTs to compute in this manner reduces the computational complexity from for the direct calculation of (19) to .
Once we have computed through the implicit quadrature rule discussed previously, we simply compute the spin spherical harmonic coefficients through (9). The complexity of this computation is , which dominates the overall complexity of the forward algorithm.
The forward algorithm may be summarised conceptually as follows, where we view the upsampling and application of weights in the spatial domain:
Although the complexity of this approach remains identical to a standard separation of variables, the continual use of FFTs reduces the constant prefactor associated with the asymptotic scaling considerably, resulting in an algorithm that may be used to compute harmonic coefficients rapidly. Furthermore, a precomputation is not required, which can otherwise necessitate very large storage requirements for high bandlimits. Finally, note that this algorithm applies to both scalar and spin functions without any change in computational complexity or computation time.
4.5 Inverse transform
The inverse algorithm follows the harmonic formulation presented in Sec. 4.1 closely. Firstly, is computed directly through (14), where we also exploit the symmetry (16) to reduce the number of computations by a factor of two. The complexity of this computation is , which dominates the overall complexity of the inverse algorithm. FFTs are then used to evaluate (13) and (12) rapidly over the extended domain, with complexities for both computations. To facilitate the efficient calculation of (12) through the use of an FFT, we compute function values on the extended domain , however we discard those values computed in the domain . The algorithm presented here to compute the inverse transform is identical to that first proposed by one of the authors [34] and subsequently adopted by Huffenberger & Wandelt [35]. The inverse algorithm may be summarised as follows:
Since we construct algorithms to perform forward and inverse spherical harmonic transforms that are theoretically exact, our construction corresponds to a novel sampling theorem on the sphere. Moreover, to represent a bandlimited signal on the sphere our sampling theorem requires less than half the number of number of samples required by other equiangular sampling theorems on the sphere [23, 26, 35] and an asymptotically identical, but smaller, number of samples that the GaussLegendre sampling theorem.
4.6 Quadrature
The construction of our sampling theorem on the sphere can be used to define an explicit quadrature rule for the integration of a function bandlimited at . This integration, which corresponds to computing the spherical harmonic coefficient , requires approximately half as many samples as needed to compute all spherical harmonic coefficients. We define the explicit quadrature weights to evaluate the following integral exactly by the finite sum:
(20) 
where . As seen from (9), the computation of requires for only. Consequently, aliasing in in all Fourier coefficients except may be tolerated, hence the number of samples required in is reduced from to . From the (reflected) convolution (19), it is apparent that the coefficients are required for all (for only). Consequently, the sampling in remains unchanged, with samples. However, only weights with argument up to are required. The (reflected) convolution thus spans the range . However, since only is of interest aliasing may be tolerated in in all Fourier coefficients except , so that zeropadding is not required before computing the (reflected) convolution as a product in the spatial domain. The absence of upsampling leads to the explicit quadrature (20), with samples on the sphere, where the weights are defined by
and where is the inverse discrete Fourier transform of the reflected weights :
The weights defined on are exactly the samples of the function defined by on and zero on , bandlimited at . The quadrature weights defined on are constructed simply by folding the contributions of on back onto the domain. Both of these weights are plotted and compared to in Fig. 1.
4.7 Real signals
In many practical applications signals observed on the sphere satisfy a reality condition. For spin signals, the reality condition is given by , which implies the conjugate symmetry condition on the spherical harmonic coefficients of the signal. When this reality condition is satisfied (for example, when considering the polarisation of the CMB), we may exploit this symmetry to recover the harmonic coefficients of a spin signal for free from the coefficients of a spin signal. For the spin case, the reality condition reduces to the standard reality condition of a scalar signal, implying . In this case, noting (14) and (4), we obtain the symmetry
We exploit these symmetries to reduce the computational cost of both the forward and inverse algorithms by an additional factor of approximately two for real spin signals.
5 Evaluation
We have implemented our fast algorithms to compute spherical harmonic transforms in double precision arithmetic, exploiting all of the symmetries discussed in Sec. 4 to optimise the implementation.
In this section we evaluate our sampling theorem and fast algorithms in terms of the number of samples required to represent a bandlimited function, numerical precomputation (or lack thereof), the recursions used to compute Wigner functions, and numerical accuracy and computation time. Finally, we evaluate our sampling theorem in the context of potential applications.
5.1 Sampling
The number of samples required to represent a bandlimited function on the sphere exactly is the fundamental property of any sampling theorem, with fewer samples desired. Both the GaussLegendre and our sampling theorems require in general samples, while the Driscoll & Healy sampling theorem requires samples. We therefore provide a reduction in number of samples by a factor of two compared to the canonical equiangular sampling theorem on the sphere. Furthermore, we require fewer samples that the GaussLegendre sampling theorem, which for small bandlimits can be significant (as discussed in Sec. 5.5). The optimal number of samples attainable by a sampling theorem on the sphere is given by the degrees of freedom in harmonic space. No sampling theorem on the sphere reaches this bound in general. However, both the GaussLegendre and our sampling theorem reach the bound in the limiting case of small (we reach the bound for the cases , the GaussLegendre sampling theorem reaches the bound for only, while the Driscoll & Healy sampling theorem never reaches the bound). In Fig. 2 we plot the number of samples against bandlimit for various sampling theorems. We also plot the position of samples for each of these sampling theorems in Fig. 3.


5.2 Precomputation
It is possible to reduce the computational burden of our fast algorithms by precomputing the Wigner functions for argument and for all required harmonic indices. Such a precomputation would require storage, similar to the storage requirements of the precomputation for Driscoll & Healy based algorithms [26, 27, 29]. For these latter approaches precomputation is essential to recover the fast algorithms with complexity below . However, for high bandlimits the storage requirements become impractical at present (recall that the precomputation at is expected to require approximately 77GB of storage). The Wigner functions may be evaluated accurately and rapidly using recursion formulae; therefore to avoid storage problems we do not perform any precomputation and instead compute Wigner functions onthefly using the method of Risbo [40].
5.3 Computing Wigner functions
Since our algorithms require Wigner functions evaluated on the entire plane for argument only, they are flexible with regard to the choice of recursion used to compute the Wigner plane for a given . For example, for the onthefly computation of Wigner functions we may use the recursion of Risbo [40] (which requires the entire plane) or Trapani & Navaza [41] (which is restricted to argument ), which are both of complexity , without altering the overall complexity of our algorithms. However, this is not the case for alternative algorithms.
To compute spherical harmonic transforms for the GaussLegendre and Driscoll & Healy sampling theorems using the algorithms described previously, it is necessary to compute Wigner function values for a single row of the Wigner plane only, but for all and all values of . For the overall algorithms to remain , the onthefly computation of the row of the Wigner plane must be performed in computations. This precludes the use of the recursions devised by Risbo [40] and by Trapani & Navaza [41], for example, both of which would result in an overall algorithm with complexity . Instead, alternative recursions must be used, such as the threeterm recursion in that goes pointwise through the Wigner plane (see e.g. (4.5) of [28]; this recursion is used in our implementation of the GaussLegendre sampling theorem). The inflexibility of these algorithms with regard to the choice of recursion used to compute Wigner functions becomes important when we study the stability of these recursions in the following section. Of course, this issue may be resolved by precomputing Wigner functions using any recursion, but as we have seen this becomes problematic at large bandlimits.
5.4 Numerical accuracy and computation time
We evaluate the numerical accuracy and computation time of our algorithms that implement our new sampling theorem, comparing them to our optimised implementation of the GaussLegendre sampling theorem and to the seminaive algorithm [27] in SpharmonicKit
The maximum absolute error is plotted against bandlimit in Fig. 4 for different sampling theorems. High numerical accuracy is achieved for all sampling theorems at moderate bandlimits, with errors on the order of the machine precision and increasing approximately linearly with bandlimit. For the GaussLegendre sampling theorem we use the gauleg routine of Numerical Recipes [43] to compute GaussLegendre node positions and weights. This method is based on an initial approximation for each node position, followed by an iterative refinement based on Newton’s method, which is likely to explain the slightly inferior error performance of the corresponding algorithms. Furthermore, the algorithms implementing both the GaussLegendre and Driscoll & Healy sampling theorems suffer from their lack of flexibility regarding the recursion used to compute Wigner functions (or associated Legendre functions for the spin case), which necessitates the use of the less accurate pointwise threeterm recursion in , rather than more accurate alternatives. Consequently, the numerical accuracy of our algorithms is superior to the optimised implementations of the alternative sampling theorems. More importantly, however, both the GaussLegendre sampling theorem and the seminaive implementation of the Driscoll & Healy sampling theorem go unstable between and , due to the instability of the pointwise threeterm Wigner recursion. As discussed in Sec. 5.3, these algorithms require a recursion with complexity (), in order to remain for onthefly computation. To resolve the instability of these algorithms it would be necessary to use the alternative recursion of Risbo [40], making them , or to perform a precomputation of Wigner functions. Neither of these solutions are feasible for high bandlimits, hence these algorithms are restricted to moderate bandlimits (unless an alternative pointwise recursion can be found that is stable to high bandlimits). Due to the flexibility of our sampling theorem with regard to Wigner recursion, we are able to use Risbo’s recursion [40], which is stable to at least , without altering the complexity of our algorithms. Finally, note that for the implementations that support spin transforms (the GaussLegendre and our sampling theorems), the error is identical (to statistical noise) for transforms of real and complex signals of different spin.
The computation time for complex signals is plotted against bandlimit in Fig. 5 for different sampling theorems. For all sampling theorems, computation time evolves as as predicted. The seminaive algorithm is slightly faster than our algorithm, which is in turn slightly faster than our optimised implementation of the GaussLegendre sampling theorem. Although we plot performance results for our algorithms using the recursion of Risbo [40], we also implemented the recursion of Trapani & Navaza [41], which we found to be approximately 20% faster than Risbo’s approach but which goes unstable between and . Nevertheless, the Trapani & Navaza [41] recursion can be used at bandlimits at or below to provide a considerable speed enhancement. In this case, at bandlimit we are approximately 25% slower than the seminaive algorithm, but twice as fast as the GaussLegendre algorithm. However, the seminaive algorithm applies for scalar functions only (and requires approximately twice as many samples on the sphere), while the alternative sampling theorems also apply directly for spin functions on the sphere. Finally, note that for the implementations that support spin transforms (the GaussLegendre and our sampling theorems), computation time is identical (to statistical noise) for transforms of signals of different spin, as predicted.
5.5 Applications
We discuss three potential applications of our sampling theorem in the fields of cosmology, neuroscience and compressive sampling (CS). For each case we highlight the enhancements that our sampling theorem will afford.
For the analysis of CMB observations, our sampling theorem and associated algorithms provide the ability to perform fast spherical harmonic transforms that are exact at the very high resolution of current and forthcoming CMB observations. Furthermore, we can compute harmonic transforms of both the temperature and polarisation of the CMB for identical cost.
The number of samples required to represent a bandlimited signal is not only of theoretical interest but also has important practical application. A sampling theorem requiring fewer samples means a bandlimited function can be measured exactly for lower cost. This is particularly important in applications where the cost of acquiring a single sample is large. In neuroscience, for example, diffusion magnetic resonance imaging (MRI) [44] is one such application, where cost is measured in terms of acquisition time. Diffusion MRI has received considerable attention recently as a noninvasive technique to image structural neuronal connectivity in the brain. In this setting, most recent acquisition strategies consider sampling on multiple spherical shells for each voxel of the brain, from which an orientation distribution function (ODF) describing the probability density of neuronal fibre directions is recovered [45, 46]. The ODFs of each voxel are then combined to recover neuronal connectivity. Given the millions of voxels generally considered, at present this imaging modality remains too time consuming for clinical use. Typically, very low bandlimits of order are considered for each spherical shell. Consequently, when adopting an exact sampling theorem even the small reduction in number of samples between our sampling theorem and the GaussLegendre approach of can have a large impact on the total cost of acquisition. For a typical example with acquisitions made on three concentric spherical shells of increasing radius, it has been shown that bandlimits of three, five and nine, respectively, are required to limit aliasing to acceptable levels [47]. For such an acquisition, the total number of samples, and thus total acquisition time, would be reduced by a factor of 13% when replacing GaussLegendre sampling with our sampling theorem. This type of enhancement is of considerable importance in order to make diffusion MRI accessible for clinical use.
The recently developed theory of CS states that it is possible to acquire sparse or compressible signals with fewer samples than standard sampling theorems would suggest [48, 49]. In these settings, the ratio of the number of required measurements to the dimensionality of the signal scales linearly with its sparsity [48]. By reducing the dimensionality of the signal in the spatial domain, our sampling theorem will enhance the performance of CS reconstruction on the sphere when compared to alternative sampling theorems. Furthermore, for sparsity priors defined in the spatial domain, such as signals sparse in the magnitude of their gradient, sparsity is also directly related to the sampling of the signal. For this class of signals, we therefore expect to see an additional enhancement in CS reconstruction performance when adopting our sampling theorem. The use of CS techniques on the sphere is likely to have widespread application for a wide range of problems, including more efficient acquisition, denoising and deconvolution on the sphere. In particular, all of these problems are faced in diffusion MRI and in analysing the CMB, which we are studying currently to evaluate in detail the enhancements provided by our sampling theorem.
6 Conclusions
We have developed a novel sampling theorem on the sphere, with corresponding fast algorithms, by associating the sphere with the torus through a periodic extension. To represent a bandlimited signal on the sphere exactly our sampling theorem requires less than half the number of samples required by other equiangular sampling theorems on the sphere and an asymptotically identical, but smaller, number of samples than the GaussLegendre sampling theorem on the sphere. The complexity of our algorithms to compute both forward and inverse transforms is , with identical scaling to a standard separation of variables. However, the continual use of FFTs reduces the constant prefactor associated with the asymptotic scaling considerably, resulting in algorithms that may be used to compute harmonic transforms rapidly. Numerical experiments have shown our algorithms to be approximately twice as fast as optimised algorithms implementing the GaussLegendre sampling theorem but approximately 25% slower than the seminaive algorithm, the most universally applicable algorithm implementing the equiangular Driscoll & Healy sampling theorem. However, the seminaive algorithm applies for scalar functions only, while our sampling theorem also applies directly for spin functions on the sphere (whereas the computation time for a spin transform using the seminaive algorithm would scale by ). Numerical experiments have also shown our algorithms to be numerically stable to bandlimits of . Conversely, the algorithms that implement the GaussLegendre and Driscoll & Healy sampling theorems on the sphere are restricted in their use of Wigner recursions and, due to the enforced use of the pointwise threeterm Wigner recursion, go unstable between bandlimits and .
Our novel sampling theorem and fast algorithms will be of practical benefit both at very high and low bandlimits. For the analysis of CMB data at very high bandlimits (), our sampling theorem yields exact spherical harmonic transforms with algorithms that are stable and very accurate. For the reconstruction of diffusion MRI images at low bandlimits (), the reduction in the number of samples required by our sampling theorem to represent a bandlimited signal may be exploited to reduce the cost of acquisition significantly. For CS applications in both of these fields and beyond, the reduction in dimensionality and sparsity afforded by our sampling theorem will enhance the performance of CS reconstruction on the sphere.
[]Jason McEwen received a B.E. (Hons) degree in Electrical and Computer Engineering from the University of Canterbury, New Zealand, in 2002 and a Ph.D. degree in Astrophysics from the University of Cambridge in 2006.
He held a Research Fellowship at Clare College, Cambridge, from 2007 to 2008, worked as a Quantitative Analysis in an Investment Bank from 2008 to 2010, and began a position as a Postdoctoral Researcher at Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland, in 2010. His research interests include wavelets on the sphere, compressed sensing, and applications of these theories in cosmology and radio interferometry.
[]Yves Wiaux received the M.S. degree in physics and the Ph.D. degree in theoretical physics from the Université catholique de Louvain (UCL), LouvainlaNeuve, Belgium, in 1999 and 2002, respectively.
He was a Postdoctoral Researcher at the Signal Processing Laboratories of the Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland, from 2003 to 2008. He was also a Postdoctoral Researcher of the Belgian National Science Foundation (F.R.S.FNRS) at the Physics Department of UCL from 2005 to 2009. He is now a Maître Assistant of the University of Geneva (UniGE), Switzerland, with joint affiliation between the Institute of Electrical Engineering and the Institute of Bioengineering of EPFL, and the Department of Radiology and Medical Informatics of UniGE. His research lies at the intersection between complex data processing (including development on wavelets and compressed sensing) and applications in astrophysics (notably in cosmology and radio astronomy) and in biomedical sciences (notably in MRI and fMRI).
Footnotes
 http://healpix.jpl.nasa.gov/
 http://www.mrao.cam.ac.uk/projects/cpac/igloo/
 http://www.glesp.nbi.dk/
 Healy et al. [27] derive a number of variants of the original Driscoll & Healy algorithm [26], including the socalled seminaive, simplesplit and hybrid algorithms. The seminaive algorithm avoids dividing (and conquering) the problem, resulting in complexity. The simplesplit algorithm is a simpler and more stable divideandconquer approach than the original algorithm but with an increased complexity of and is less stable than the seminaive approach. The splitting required by the simplesplit algorithm is costly (in terms of execution time rather than asymptotic complexity), thus for a bandlimit of the seminaive algorithm is greater than two times faster than the simplesplit algorithm [27]. The hybrid algorithm attempts to mitigate the slow execution of the simplespit algorithm and the higher complexity of the seminaive algorithm by splitting the problem between them. The hybrid algorithm appears to achieve a good compromise between stability and efficiency. However, the overall complexity of this algorithm is not clear since it depends on the split between the seminaive and simplesplit algorithms and on user specified parameters.
 http://www.cs.dartmouth.edu/~geelong/sphere/
 We adopt the Euler convention corresponding to the rotation of a physical body in a fixed coordinate system about the , and axes by , and respectively.
 We check that this periodic extension does not impose discontinuities at the poles . To avoid discontinuities we require , which follows trivially for even. From (6) we find , which, due to the continuity of the Wigner functions, implies for odd. Combined with the representation, for odd, this implies for odd; hence our periodic extension does not impose any discontinuity.
 The structure of our algorithms suggest multiple transforms of different spin may in theory be computed simultaneously at lower cost than consecutive computation (since Wigner functions do not need to be recomputed and some computations are independent of spin). However, for simplicity we do not incorporate this optimisation in our current implementation.
 http://www.fftw.org/
 http://www.jasonmcewen.org/
 It is well know that GaussLegendre quadrature may be used to construct an exact sampling theorem on the sphere. GaussLegendre quadrature with points is exact only for a polynomial integrand of order less than or equal to . Since neither the associated Legendre functions nor the Wigner functions are polynomials in , it does not follow immediately that GaussLegendre quadrature results in an exact harmonic transform for scalar and spin signals on the sphere bandlimited at ; nevertheless, this is indeed the case. We prove the result for spin signals on the sphere, thus the scalar case will follow simply by setting . For samples in , we must simply prove that the integrand is polynomial in of maximum degree less than or equal . From inspection of (3.3), we may write the Wigner functions as a polynomial of degree , multiplied by . Similarly, we may write as a polynomial of degree , multiplied by the same factor. The integrand is therefore polynomial with overall degree , which reaches a maximum of . GaussLegendre quadrature with samples in may thus be used to compute exact spherical harmonic transforms of scalar and spin functions bandlimited at .
 http://www.cs.dartmouth.edu/~geelong/sphere/
References
 R. Ramamoorthi and P. Hanrahan, “A signal processing framework for reflection,” ACM Transactions on Graphics, vol. 23, no. 4, pp. 1004–1042, 2004.
 D. L. Turcotte, R. J. Willemann, W. F. Haxby, and J. Norberry, “Role of membrane stresses in the support of planetary topology,” J. Geophys. Res., vol. 86, pp. 3951–3959, 1981.
 M. A. Wieczorek, “Gravity and topology of terrestrial planets,” submitted to Treatise on Geophysics, 2006.
 M. A. Wieczorek and R. J. Phillips, “Potential anomalies on a sphere: applications to the thickness of the lunar crust,” J. Geophys. Res., vol. 103, pp. 383–395, 1998.
 P. Audet, “Directional wavelet analysis on the sphere: application to gravity and topography of the terrestrial planets,” J. Geophysical Research, in press, 2010.
 F. J. Simons, F. A. Dahlen, and M. A. Wieczorek, “Spatiospectral concentration on a sphere,” SIAM Rev., vol. 48, no. 3, pp. 504–536, 2006.
 S. Swenson and J. Wahr, “Methods for inferring regional surfacemass anomalies from GRACE measurements of timevariable gravity,” J. Geophys. Res., vol. 107, pp. 2193–278, 2002.
 K. A. Whaler, “Downward continuation of Magsat lithosphere anomalies to the Earth’s surface,” Geophys. J. R. Astr. Soc., vol. 116, pp. 267–278, 1994.
 C. H. Choi, J. Ivanic, M. S. Gordon, and K. Ruedenberg, “Rapid and stable determination of rotation matrices between spherical harmonics by direct recursion,” J. Chem. Phys., vol. 111, no. 19, pp. 8825–8831, 1999.
 D. W. Ritchie and G. J. L. Kemp, “Fast computation, rotation and comparison of low resolution spherical harmonic molecular surfaces,” J. Comput. Chem., vol. 20, no. 4, pp. 383–395, 1999.
 C. L. Bennett, A. Banday, K. M. Gorski, G. Hinshaw, P. Jackson, P. Keegstra, A. Kogut, G. F. Smoot, D. T. Wilkinson, and E. L. Wright, “4year COBEDMR cosmic microwave background observations: maps and basic results,” Astrophys. J. Lett., vol. 464, no. 1, pp. 1–4, 1996.
 N. Jarosik, C. L. Bennett, J. Dunkley, B. Gold, M. R. Greason, M. Halpern, R. S. Hill, G. Hinshaw, A. Kogut, E. Komatsu, D. Larson, M. Limon, S. S. Meyer, M. R. Nolta, N. Odegard, L. Page, K. M. Smith, D. N. Spergel, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, “Sevenyear Wilkinson Microwave Anisotropy Probe (WMAP) observations: sky maps, systematic errors, and basic results,” Astrophys. J., submitted, 2010.
 E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, “Sevenyear Wilkinson Microwave Anisotropy Probe (WMAP) observations: cosmological interpretation,” Astrophys. J., submitted, 2010.
 Planck collaboration, “ESA Planck blue book,” ESA, Tech. Rep. ESASCI(2005)1, 2005.
 S. A. Orszag, “Fast eigenfunction transforms,” Adv. Math. Supp. Studies, vol. 10, pp. 23–30, 1986.
 R. K. Beatson and L. Greengard, “A short course on fast multipole methods,” in Wavelets, multilevel methods and elliptic PDEs, J. Ainsworth, J. Levesley, W. Light, and M. Marletta, Eds. Oxford University Press, 1997, pp. 1–37.
 B. K. Alpert and V. Rokhlin, “A fast algorithm for the evaluation of Legendre expansions,” SIAM J. Sci. Stat. Comput., vol. 12, no. 1, pp. 158–179, 1991.
 R. Suda and M. Takami, “A fast spherical harmonics transform algorithm,” Math. Comput., vol. 71, no. 238, pp. 703–715, 2002.
 M. J. Mohlenkamp, “A fast transform for spherical harmonics,” Ph.D. dissertation, Yale University, 1997.
 ——, “A fast transform for spherical harmonics,” J. Fourier Anal. and Appl., vol. 5, no. 2/3, pp. 159–184, 1999.
 K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “Healpix – a framework for high resolution discretization and fast analysis of data distributed on the sphere,” Astrophys. J., vol. 622, pp. 759–771, 2005.
 R. G. Crittenden and N. G. Turok, “Exactly azimuthal pixelizations of the sky,” ArXiv, 1998.
 W. Skukowsky, “A quadrature formula over the sphere with application to high resolution spherical harmonic analysis,” J. Geodesy, vol. 60, pp. 1–14, 1986, 10.1007/BF02519350. [Online]. Available: http://dx.doi.org/10.1007/BF02519350
 M. H. Payne, “Truncation effects in geopotential modelling,” Analytical Mechanics Associates, Seabrook, Maryland, Tech. Rep., 1971.
 A. G. Doroshkevich, P. D. Naselsky, O. V. Verkhodanov, D. I. Novikov, V. I. Turchaninov, I. D. Novikov, P. R. Christensen, and L. Y. Chiang, “Gauss–Legendre Sky Pixelization (GLESP) for CMB maps,” Int. J. Mod. Phys. D., vol. 14, no. 2, pp. 275–290, 2005.
 J. R. Driscoll and D. M. J. Healy, “Computing Fourier transforms and convolutions on the sphere,” Advances in Applied Mathematics, vol. 15, pp. 202–250, 1994.
 D. M. J. Healy, D. Rockmore, P. J. Kostelec, and S. S. B. Moore, “FFTs for the 2sphere – improvements and variations,” J. Fourier Anal. and Appl., vol. 9, no. 4, pp. 341–385, 2003.
 P. Kostelec and D. Rockmore, “FFTs on the rotation group,” J. Fourier Anal. and Appl., vol. 14, pp. 145–179, 2008. [Online]. Available: http://dx.doi.org/10.1007/s0004100890135
 Y. Wiaux, L. Jacques, and P. Vandergheynst, “Fast spin 2 spherical harmonics transforms,” J. Comput. Phys., vol. 226, p. 2359, 2005.
 P. J. Kostelec, D. K. Maslen, D. M. Healy, and D. N. Rockmore, “Computational harmonic analysis for tensor fields on the twosphere,” J. Comput. Phys., vol. 162, no. 2, pp. 514–535, 2000.
 Y. Wiaux, L. Jacques, P. Vielva, and P. Vandergheynst, “Fast directional correlation on the sphere with steerable filters,” Astrophys. J., vol. 652, pp. 820–832, 2006.
 D. Potts, G. Steidl, and M. Tasche, “Fast and stable algorithms for discrete spherical Fourier transforms,” Linear Algebra Appl., vol. 275/276, no. 1–3, pp. 433–450, 1998.
 G. A. Dilts, “Computation of spherical harmonic expansion coefficients via FFTs,” J. Comput. Phys., vol. 57, no. 3, pp. 439–453, 1985. [Online]. Available: http://www.sciencedirect.com/science/article/B6WHY4DDR4YFW8/2/f3aa3657495be2d78ddc61d455f67e06
 J. D. McEwen, “Fast, exact (but unstable) spin spherical harmonic transforms,” All Res. J. Phys., vol. 1, no. 1, 2011.
 K. M. Huffenberger and B. D. Wandelt, “Fast and exact spins spherical harmonic transforms,” Astrophys. J. Supp., vol. 189, pp. 255–260, 2010.
 E. T. Newman and R. Penrose, “Note on the BondiMetznerSachs group,” J. Math. Phys., vol. 7, no. 5, pp. 863–870, 1966. [Online]. Available: http://link.aip.org/link/?JMP/7/863/1
 M. Zaldarriaga and U. Seljak, “Allsky analysis of polarization in the microwave background,” Phys. Rev. D., vol. 55, no. 4, pp. 1830–1840, Feb 1997.
 J. N. Goldberg, A. J. Macfarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, “Spin spherical harmonics and ,” J. Math. Phys., vol. 8, no. 11, pp. 2155–2161, 1967. [Online]. Available: http://link.aip.org/link/?JMP/8/2155/1
 D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum theory of angular momentum. Singapore: World Scientific, 1989.
 T. Risbo, “Fourier transform summation of Legendre series and functions,” J. Geodesy, vol. 70, no. 7, pp. 383–396, 1996.
 S. Trapani and J. Navaza, “Calculation of spherical harmonics and Wigner d functions by FFT. Applications to fast rotational matching in molecular replacement and implementation into AMoRe,” Acta Crystallographica Section A, vol. 62, no. 4, pp. 262–269, 2006. [Online]. Available: http://dx.doi.org/10.1107/S0108767306017478
 A. F. Nikiforov and V. B. Uvarov, Classical Orthogonal Polynomials of a Discrete Variable. Berlin: SpringerVerlag, 1991.
 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Fannery, Numerical recipes in Fortran 77, 2nd ed. Cambridge: Cambridge University Press, 1992.
 H. JohansenBerg and T. E. J. Behrens, Diffusion MRI: From quantitative measurement to invivo neuroanatomy. San Diego: Academic Press, 2009.
 D. S. Tuch, “Qball imaging,” Magn. Reson. Med., vol. 52, no. 6, pp. 1358–1372, 2004.
 E. J. CanalesRodríguez, L. MelieGarcía, and Y. IturriaMedina, “Mathematical description of qspace in spherical coordinates: Exact qball imaging,” Magn. Reson. Med., vol. 61, no. 6, pp. 1350–1367, 2009.
 A. Daducci, J. D. McEwen, D. V. D. Ville, J.P. Thiran, and Y. Wiaux, “Harmonic analysis of spherical sampling in diffusion MRI,” in 19th Annual Meeting of International Society for Magnetic Resonance in Medicine, 2011.
 E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory, vol. 52, no. 2, pp. 489–509, feb 2006.
 D. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory, vol. 52, no. 4, pp. 1289–1306, apr 2006.