An Optimal-Dimensionality Sampling Scheme on the Sphere with Fast Spherical Harmonic Transforms
We develop a sampling scheme on the sphere that permits accurate computation of the spherical harmonic transform and its inverse for signals band-limited at using only samples. We obtain the optimal number of samples given by the degrees of freedom of the signal in harmonic space. The number of samples required in our scheme is a factor of two or four fewer than existing techniques, which require either or samples. We note, however, that we do not recover a sampling theorem on the sphere, where spherical harmonic transforms are theoretically exact. Nevertheless, we achieve high accuracy even for very large band-limits. For our optimal-dimensionality sampling scheme, we develop a fast and accurate algorithm to compute the spherical harmonic transform (and inverse), with computational complexity comparable with existing schemes in practice. We conduct numerical experiments to study in detail the stability, accuracy and computational complexity of the proposed transforms. We also highlight the advantages of the proposed sampling scheme and associated transforms in the context of potential applications.
2-sphere (unit sphere), spherical harmonic transform, sampling, harmonic analysis, spherical harmonics.
Signals are inherently defined on the sphere in a variety of fields of science and engineering. These include geodesy , cosmology , computer graphics , medical imaging , astrophysics , quantum chemistry , wireless communication , acoustics  and planetary science , to name a few. In signal processing analysis on the sphere (e.g.,[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 9, 26, 27, 28]) the signal is often analysed in both the spherical (spatial) domain and harmonic (spectral) domain. The transformation from spatial to spectral is through the spherical harmonic transform (SHT) (see, e.g., [29, 30, 10, 11, 28]), which is the well-known counterpart of the Fourier transform. For example, analysis of signals in the spectral domain through the SHT has been instrumental in refining the standard cosmological model and in the study of the anisotropies in the cosmic microwave background (CMB) . Consequently, the ability to compute the SHT of a signal is of significant importance. Furthermore, since data-sets on the sphere can be of considerable size, and the cost of acquiring samples on the sphere can be large[8, 32], the computation of the SHT of the signal should require the minimum possible number of samples, and be computationally accurate and efficient.
The development of sampling schemes on the sphere and computationally efficient methods to compute the spherical harmonic transform from samples has been investigated extensively in the literature[10, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 11]. Sampling schemes and their associated SHT computational methods can be evaluated by three key criteria: (1) the number of samples, defined as the spatial dimensionality; (2) the computational complexity; and (3) the numerical accuracy. In this work, we propose a sampling scheme for band-limited signals on the sphere which requires the same number of samples on the sphere as the number of the degrees of freedom of the signal in harmonic space. Furthermore, we develop an accurate method to compute the SHT with complexity scaling, in practice, comparable with the existing schemes. We first review the developments made in the literature followed by a summary of the contributions of this paper.
1.1 Relation to Prior Work
Among existing sampling schemes in the literature, iso-latitude sampling schemes (e.g.,[10, 36, 37, 35, 38, 40, 42, 11, 33, 39]), where the samples along longitude are taken over iso-latitude rings (annuli), enable a separation of variables in the computation of the SHT, which results in a reduction in computational complexity. For the computation of spherical harmonic transforms, sampling theorems have been constructed [10, 37, 42, 11], which lead to theoretically exact SHTs, in addition to other numerical approaches, such as approximate quadrature [36, 35]), least squares [33, 39] or spherical designs [43, 44], which nevertheless often lead to accurate transforms. We focus our attention on an iso-latitude sampling scheme that facilitates accurate computation of the SHT of a signal that is band-limited at (formally defined in Section 2.2). For the accurate computation of the SHT of a signal band-limited at , the optimal spatial dimensionality, denoted by , attainable by any sampling scheme on the sphere is given by , which is the number of degrees of freedom in harmonic space. Existing schemes either require or samples and therefore do not achieve the optimal spatial dimensionality. We refer the reader to  for a more comprehensive review of existing sampling schemes.
An exact method to compute the SHT, based on a sampling theorem on the sphere, was developed by Driscoll and Healy in  exploiting an equiangular sampling comprised of iso-latitude rings of points, where the number of points along longitude, in each ring, is the same and equal to . Thus the spatial dimensionality of Driscoll and Healy sampling scheme is . The well-known Gauss-Legendre quadrature on the sphere [45, 46] may also be used to construct a sampling theorem and exact SHT from samples on the sphere. The placement of iso-latitude rings is given by the roots of the Legendre polynomials of order , as dictated by Gauss-Legendre quadrature, and the number of points in each ring remains . More recently, a new sampling theorem based on an equiangular sampling scheme has been proposed by McEwen and Wiaux, which achieves spatial dimensionality , and requries fewer samples than the Gauss-Legendre approach. For all of these sampling schemes the associated spherical harmonic transforms (that are stable and do not require precomputation) have computational complexity of . We note that an algorithm was developed for the Driscoll and Healy sampling theorem  with complexity , but it requires precomputation and storage. The precomputation is practicable for applications in acoustics , quantum chemistry , medical imaging , where the band-limit is of the order . However, the precomputation becomes infeasible for applications in astrophysics  and cosmology , where the band-limit is of the order , as the precomputation requires GB of storage for , which scales to approximately GB for the band-limit [42, 11].
SHTs using least squares approaches have also been developed for equiangular sampling schemes [33, 39], which require samples and achieve good accuracy. However, a naive application of least squares is computationally inefficient since the computational complexity to compute SHTs scales as . However, a separation of variables can be employed to reduce the complexity to . Moreover, a least squares method can also be developed to compute SHTs using the optimal number () of spatial samples (either regularly or irregularly) distributed over the sphere, but the complexity of such method scales with and therefore least squares becomes computationally infeasible even for small band-limits. We note that complexity of the least squares approach using samples cannot be reduced by a separation of variables since aliasing errors are introduced if the number of samples along longitude in each iso-latitude ring is smaller than . Furthermore, the accuracy and stability of the least squares approach cannot be guaranteed. For example, the reconstruction error in the least squares system using samples is (obtained through numerical experiments on the regular grid) for band-limit and the error grows with the band-limit.
We also note the work on spherical designs for computing integrals over the sphere using quadrature based on uniform weighting (see [43, 44] for a comprehensive review). Is it numerically, but rigorously, proved that the computation of the SHT of a band-limited signal using spherical designs can be performed with samples . It is also conjectured that in fact samples may be used , although this is not proved. Spherical designs with samples have been constructed successfully for band-limit up to only .
A summary of the contributions of this paper are as follows.
We develop a sampling scheme on the sphere that permits accurate computation of the SHT for band-limited signals, attaining the optimal spatial dimensionality of .
We develop a computationally efficient method to compute the SHT and its inverse (called the forward and inverse SHT in the sequel) using our optimal spatial dimensionality sampling scheme, which has complexity with scaling, in practice, comparable to the existing methods, which do not achieve optimal spatial dimensionality. Furthermore, once the sample positions are determined no additional precomputation is required.
We characterize the numerical accuracy and computational complexity of our proposed SHT as a function of the band-limit parameter demonstrating its feasibility on large real-world data-sets.
1.3 Paper Organization
We review the mathematical background and harmonic analysis on the sphere in Section 2. In Section 3, we propose the sampling scheme on the sphere, which achieves optimal spatial dimensionality, and present the associated novel SHT. We evaluate the numerical accuracy of the proposed SHT, present its computational complexity analysis and outline potential applications in Section 4. The concluding remarks are made in Section 5.
In this section we review the mathematical background for signals and harmonic analysis on the sphere.
2.1 Signals on the Sphere
In this work, we consider square integrable complex functions of the form , defined on unit sphere , where denotes the Euclidian norm, denotes the co-latitude and denotes the longitude. The inner product of two functions and defined on is defined as 
where denotes the complex conjugate, denotes the differential area element on the sphere and the integration is carried out over the sphere, that is, . With the inner product in (1), the space of square integrable complex valued functions on the sphere forms a complete Hilbert space . Also, the inner product in (1) induces a norm . We refer to the functions with finite induced norm as signals on the sphere.
2.2 Harmonic Analysis on the Sphere
The Hilbert space is separable and the spherical harmonic functions (or spherical harmonics for short) [28, 30, 29] of all degrees and orders form the archetype complete orthonormal set of basis functions. By completeness, any signal can be expanded as
is the spherical harmonic coefficient of degree and order . Some background properties of spherical harmonics used in this work are given in Appendix 6.
The signal is defined to be band-limited at degree if for . The set of bandlimited signals forms an dimensional subspace of , which is denoted by .
2.3 Spin Functions on the Sphere
The spin functions on the sphere, denoted by and parameterized by integer spin , are special functions which are defined by their behaviour under local rotations. The local rotation by , rotates the spin function by in the tangent plane formed at a point on the sphere characterized by and . Under such rotation, the rotated spin function is related to the original spin function through
The spin spherical harmonics (sometimes also referred to as spin weighted spherical harmonics), denoted by and defined for degree , order and spin form a complete set of basis functions for spin functions on the sphere (see Appendix 6 for the definition of ). Spin spherical harmonics also satisfy the property in (4) and therefore serve as a more suitable choice of basis functions for the following expansion of spin functions
The spin function is said to be band-limited at if for all . The set of such band-limited spin functions for each form a subspace of and is denoted by . Furthermore, we note that , that is, the spin spherical harmonic become the standard spherical harmonic for . In the sequel, any reference to a function (or signal) and spherical harmonic means finite energy scalar function () on the sphere and scalar spherical harmonic () respectively, unless otherwise explicitly stated that the signal or spherical harmonic under consideration is spin weighted.
3 Optimal Sampling Scheme and Novel Spherical Harmonic Transform
We first consider the spherical harmonic transform (SHT) of band-limited scalar functions , for which the summation over degree in (2) is truncated to , that is,
We first present our sampling scheme for the discretization of a band-limited signal on the sphere. Then we develop the proposed forward SHT which determines the spherical harmonic coefficients for and using the discretized signal. We also present the inverse SHT to determine the signal efficiently over the proposed sampling scheme from its spherical harmonic coefficients. We later show that the proposed sampling scheme and SHTs (both forward and inverse) are also applicable for band-limited spin functions .
3.1 Proposed Sampling Scheme
We propose an iso-latitude sampling of the sphere. Define the indexed vector as
which consists of points along . Also define as a vector of arbitrary points along for . We shortly present the location of these sample points. For discretization along , we consider equally spaced sampling points along for each . Define be a vector of equally spaced sampling points along in the ring placed at , given by
In this way, we will have iso-latitude rings of sampling points along for each sampling point along , where the number of points, , in each ring depends on the location of the ring along latitude. We note that the number of samples in the proposed sampling scheme attains the optimal spatial dimensionality, that is,
which also represents the number of degrees of freedom in harmonic space for a signal band-limited at . We first develop the forward and inverse SHTs and later provide details about the location of the samples in the vector .
3.2 Forward Spherical Harmonic Transform – Formulation
We develop the forward SHT to compute the spherical harmonic coefficients of a band-limited signal sampled over the samples of our sampling scheme. First, we first develop the necessary mathematical formulation and later we present the philosophy of our approach and develop the forward SHT. For order , define a vector
for each , where denotes scaled associated Legendre functions (see Appendix 6 for the definition of associated Legendre functions and spherical harmonic )
By defining a matrix as
and a vector containing spherical harmonic coefficients of order given by
we can write as follows
where the vector contains the values for . It is possible to recover by inverting this system, which is elaborated shortly.
Using the formulation in (14), number of spherical harmonic coefficients of order (or ) contained in the vector (or ) can be determined by first computing (or ) over number of samples or equivalently by evaluating (or ) for all and the matrix in (12) and then solving (14). We note that (using (30)), is not dependent on the signal and must be chosen such that can be inverted accurately, which enables the accurate computation of (or ).
3.3 Forward Spherical Harmonic Transform – Philosophy
Following Remark 1, we need to compute formulated in (3.2) for all and for given . Using an FFT, for each can be computed exactly by evaluating the integral as a summation, provided the following two conditions are satisfied:
, for all and for all , which ensures that the univariate signal along for each is band-limited at , given the complex exponentials as basis functions along , and
there are at least sample points in each -ring placed for all .
These conditions and the following remark form the foundation of our proposed method to compute the forward SHT.
For a band-limited signal on the sphere given by (6), we have contributions from the complex exponentials for ; thus we require samples in each of the rings placed at in order to compute correctly regardless of the choice of . However, if the spherical harmonic coefficients of all degrees (and orders) greater than are known for some , their contributions can be removed. Such a removal makes the signal band-limited along with respect to the contributions of the complex exponentials for . Consequently, can be computed correctly by taking an FFT over only samples (instead of samples) in a ring placed at some . We elaborate this philosophy below.
Fig. 1 shows the graphical representation of the spectral domain of a signal band-limited at . There is one order spherical harmonic coefficient and one order spherical harmonic coefficient , which can be determined using (14) by first computing and using an FFT over only one ring of samples along placed at . It should be noted that an FFT over samples, in fact, computes correctly for all orders . Once is computed, the signal at the other sample positions (for all and associated samples along in each ring) can be updated as
denotes the part of the signal composed of contribution of spherical harmonics of order and and all degrees . Once the signal is updated at other sample positions as given in (15), there is no contribution of spherical harmonics of orders and or equivalently there is no contribution of complex exponentials and in the signal. Thus, following the conditions stated earlier in this subsection, we need (instead of ) samples along in the rings placed at all . We only require samples along the -ring placed at to determine the spherical harmonic coefficients of order and , which once computed can be used to update the signal as
at other sample positions for all and samples in the associated samples along in each ring. Proceeding in this similar manner, we note that the forward SHT can be computed using number of samples along each ring placed at . Using the philosophy mentioned above, we summarize the proposed forward SHT in the form of the procedure below.
3.4 Inverse Spherical Harmonic Transform
The inverse SHT computes the signal from its spherical harmonic coefficients. Using the separation of variables technique (also adopted in [10, 40, 38, 11]), changing the order of summation in (2) and using (3.2) and (3.3), we write the inverse SHT as follows
where and are the sample points belonging to the proposed sampling scheme. The inverse SHT can be computed by the following procedure.
3.5 Placement of Samples along Co-latitude
For the representation of band-limited signals with band-limit we yet need to choose the location of samples, , along co-latitude for the sampling scheme presented in Section 3.1, where a ring of equiangular sample points along is placed at each . The simplest choice is to use the equiangular set of samples given by
for the placement of the rings, where rings are placed such that the rings with greater numbers of samples are placed nearer to the equator . In this setting, the ring of sample is placed at , the ring of samples is placed at and so on. For such a placement of samples, the vector for band-limit is given by
As an example, the samples on the sphere for this scheme are shown in Fig. 2 for .
We note that the proposed forward SHT requires the samples in the vector to be chosen such that the matrix becomes invertible so that the system in (14) can be inverted accurately. The sampling along the co-latitude as given in (20), although an attractive choice, may not be appropriate as the matrix may become ill-conditioned. For example, the condition number (ratio of the largest eigenvalue to the smallest eigenvalue), denoted by , of the matrix , constructed with the sample positions given in (20), for and for different values of is plotted in Fig. 3(a) and the samples in a vector are shown in Fig. 3(b). Furthermore, the maximum of the condition number, denoted by over for different values of band-limit is plotted in Fig. 3(c), where it can be observed that at least one matrix for becomes more ill-conditioned for larger band-limit .
In order to address this issue, we propose the following recipe to determine the optimal ordering of samples in a vector , which we refer to as the condition number minimization method. For a given and samples in given by (19), the vector is constructed as follows:
Choose farthest from the poles, which is a natural choice for the ring of samples along .
For each , choose from the set , given in (19), which minimizes the condition number of the matrix .
Such a placement of samples along co-latitude ensures the robust inversion of the system given by (14), thus resulting in an accurate computation of spherical harmonic coefficients by the proposed forward SHT. As an illustration, we again plot the condition number of the matrix obtained using the optimal sample positions for in Fig. 4(a), the optimal sample positions in Fig. 4(b) and the maximum condition number ) for different band-limits in Fig. 4(c). In comparison to the plots in Fig. 3, the condition number is significantly smaller for the case of optimal sampling, which leads to an accurate implementation of the proposed forward SHT. Furthermore, when computing multiple harmonic transforms for different , the optimal position of samples in a vector needs to be computed once only for each . The vector can be stored in a double precision with the storage requirement of only KB for , for example, and approximately MB for all band-limits . Moreover, we highlight that this storage determines the placement of iso-latitude rings on the sphere and is required to be known for either forward or inverse transform. Once the sample positions are known, we do not require any further precomputation for the computation of SHTs. We discuss the computational complexity of the proposed transforms later in the paper.
3.6 Alternative Placement of Samples along Co-latitude
The equiangular samples in the set are placed along co-latitude according to a uniform measure . Alternatively, the samples along co-latitude can be placed according to different measures. For example, in the context of compressive sensing, it has been proved in  that a sparse (in spectral domain) band-limited signal can be recovered from fewer measurements if samples are drawn from the measure , compared to sampling with respect to the uniform measure , which in turn has been shown by  to require fewer samples than sampling with respect to the measure .
We compare the equiangular placement of rings with the placement of rings according to the measures and . Define and for as sets of samples along co-latitude, where samples are placed according to the measures and , respectively. The sets and are constructed by choosing (South Pole) and using the following relation between the consecutive samples
for . Using the sets and , we construct the indexed vectors and , respectively, similar to in (20) defined for the set in (19), to determine the sample positions of iso-latitude rings such that the ring with samples is placed nearer to the equator.
For , we show the sample positions in Fig. 5(a) and the condition number of the matrix , constructed with the sample positions for different values of in Fig. 5(c), where it can be observed that the placement of rings according to the measure results in greater ill-conditioning of the matrices as compared to the placement of rings according to the uniform measure . We also optimize the sample positions by applying the proposed condition number minimization method. The optimal sample positions are shown in Fig. 5(b) and the condition number of the matrix , obtained using the optimal sample positions for different values of , is also plotted in Fig. 5(c). Since the matrices for the original are highly ill-conditioned, the proposed condition number minimization method, that performs the re-ordering of the sample positions along co-latitude, does not find an ordering that significantly improves the ill-conditioning of the matrices.
We also carry out a similar analysis for the sample positions shown in Fig. 6(a). The optimal sample positions obtained by applying the condition number minimization method is shown in Fig. 6(b) and the condition number of the matrix using the sample positions or optimal sample positions is plotted in Fig. 6(c) for different values of , which illustrates that the placement of rings according to the measure also results in greater ill-conditioning of the matrices as compared to the placement of rings according to the uniform measure . Thus, we conclude that the use of equiangular placement (with uniform measure ) of samples along co-latitude in the proposed sampling scheme performs well compared to the use of sampling methods which place the samples according to the measures and . We note that the equiangular placement of samples is not the only choice to place samples along co-latitude and more sophisticated sampling methods can be developed. For example, the sampling scheme that also takes into account the location of samples in each ring along longitude can be designed. However, the development of such a design is beyond the scope of the current paper and is considered as a direction for future research.
3.7 Computation of Spherical Harmonics
The implementation of both forward and inverse transforms require the computation of scaled associated Legendre functions for all degrees and orders and for all . Different recursion relations can be used for the computation of associated Legendre functions for given . For example, we can use the recursion proposed by Risbo  that computes for all orders for given degree in each step of recursion, or alternatively, we can use the three-term recursion which grows with degree and recursively computes for all for a given .
We note that the forward transform iterates over different values of and uses (which is composed of of different and ) in each iterative step for the computation of spherical harmonic coefficients in a vector (or ) as given in (14). Furthermore, the computation of in the implementation of the inverse transform also requires for all for a given . Therefore the three-term recursion is a natural choice to compute associated Legendre functions in our proposed transforms. For given and , the three term recursion relation is given by
which grows with the degree for given with the following initial condition for
and symmetry relation which follows from (32) given in the appendix
The variant of recursion relation in (21) has been adopted in the literature for the computation of associated Legendre polynomials [37, 40]. As demonstrated in , the three-term recurrence relation in (21) is stable when the recurrence is carried out in the direction of increasing , provided the initial condition in (22) is computed accurately (either by using higher than double precision arithmetic or by an adaptive rescaling).
3.8 Extension to Spin Functions on Sphere
By comparing the expansion of a spin function into spin spherical harmonics in (5) with the expansion of the non-spin standard function in (2), we note that the forward and inverse transforms developed for non-spin functions are also applicable to band-limited spin functions with the following associations
The extension of forward and inverse transforms to the spin functions require the computation of spin weighted spherical harmonics , which can be carried out using the recurrence relation given in Appendix 6, which is the generalized version of the recurrence relation in (21). We do not further investigate the application of our transforms to spin function and limit our explorations for the non-spin standard functions in the rest of the paper.
4 Analysis of Proposed Sampling Scheme and Transforms
In this section we evaluate the proposed sampling scheme and associated SHTs using the following criteria: (1) the number of samples required to accurately represent a band-limited signal; (2) the computational complexity of the associated forward and inverse transforms; and (3) the accuracy of the forward and inverse transforms. We have carried out the implementation of the proposed transforms in double precision arithmetic. The code to compute the scaled associated Legendre function for given , and using the recursion relation in (21) is written in C in order to speed up the computation, and uses an adaptive rescaling so that double precision arithmetic is accurate. The forward and inverse transforms, outlined as procedures in the previous section, are implemented in MATLAB.
4.1 Numerical Accuracy
We analyse the numerical accuracy of our forward and inverse transforms that implement our proposed optimal sampling scheme on the sphere. The accuracy of the proposed transforms means that the inverse (or forward) SHT of any band-limited signal followed by the forward (or inverse) SHT yields the same band-limited signal, with error on the order of the numerical precision. In order to evaluate the numerical accuracy of the proposed SHTs, we carry out two numerical experiments, before further studying the error distribution in harmonic space.
Experiment 1 (Spectral-Spatial-Spectral)
In our first experiment, we generate the test signal spherical harmonic coefficients for with real and imaginary parts uniformly distributed in the interval . The inverse SHT is then used to synthesise the band-limited test signal in the spatial domain over the samples of our sampling scheme, followed by the forward SHT to compute the reconstructed spherical harmonic coefficients, denoted by . The experiment is repeated 10 times and the average values for the maximum error and the mean error , given by
are recorded and plotted in Fig. 7 for band-limits in the range .
Experiment 2 (Spatial-Spectral-Spatial)
In the second experiment to test the numerical accuracy of the proposed transforms, we randomly generate the complex valued band-limited test signal with real and imaginary parts uniformly distributed in the interval over the samples of our sampling scheme. The forward SHT, followed by inverse SHT is applied on the test signal to obtain the reconstructed signal . We repeat the experiment 10 times and the obtain the average values for the maximum error and the mean error between the test signal and reconstructed signal , defined as
where the sum is over all points in the proposed sampling scheme. Both and are plotted in Fig. 8 for the band-limit in the range .
It can be observed that both the maximum error and the mean error grows quadratically with the band-limit , which is due to the computational flow of the proposed transform. The proposed transform sequentially computes the spherical harmonic coefficients , first for order , proceeding to order . The computation of the spherical harmonic coefficients of order requires knowledge of all coefficients with order greater than because the forward transform eliminates the effect of all coefficients greater than from the sample positions in the rings placed at to avoid aliasing, before taking FFTs along these rings. Any error introduced in the computation of the spherical harmonic coefficients of order propagates in the computation of coefficients of order less than . In order to further elaborate, we plot the error for and averaged over 10 realization of experiment in Fig. 9, where it can be observed that error is comparatively smaller for higher order coefficients and increases as order decreases from to .
4.2 Why is Spatial Dimensionality Important?
The fundamental property of any sampling scheme is the number of samples required to accurately represent a band-limited signal. The existing sampling schemes in the literature [45, 46, 10, 37, 11] that support an accurate SHT do not attain the optimal spatial dimensionality on the number of samples, as also highlighted earlier. In comparison, our proposed sampling scheme and associated transforms require samples, which is the optimal spatial dimensionality attainable by any sampling scheme since the band-limited signal belongs to the dimensional subspace .
Now, we briefly discuss the significance of achieving optimal spatial dimensionality. In addition to the practical considerations [8, 32, 54] that desire fewer samples for the representation of band-limited signals, we highlight that experiment 2, composed of forward transform of a signal randomly generated over samples followed by the inverse transform, yields the original signal in the spatial domain. This is not the case for existing sampling schemes since a random signal over samples in the spatial domain, where typically  or , may not belong to the dimensional subspace of band-limited signals with band-limit .
4.3 Computational Complexity Analysis
First, we analyse the computational complexity of the proposed forward SHT. Following the forward SHT procedure, the complexity to compute for each , which only requires one point FFT, is . Using the recursive relation in (21), for all and for all and for each and can be computed in time.
Since the computation of requires solving the system in (14), that can be carried out naively using the least squares approach with complexity . However, the system in (14) can be solved more efficiently in practice by employing fast algorithms. For example, the system of size can be solved in , instead of , using the algorithm of . Once is computed, the effect of higher order spherical harmonics is removed from the signal, which can be carried out in two steps: 1) evaluation of given in (3.3), which can be done asymptotically in for each ; and (2) updating the signal as for all and all associated sampling points along , which is again performed with complexity . Since these operations need to be repeated for each , and is of the order , the complexities mentioned above are scaled by and therefore the overall asymptotic complexity of the forward transform scales naively as . The dominant factor is due to the inverting of system in (14) only, which can be more efficiently implemented in practice and we demonstrate later in this section that the complexity of both the inversion of the system and the forward transform scales close to in practice. If the matrices and the inverse matrices (or pseudo-inverse) for all are pre-computed, the theoretical complexity reduces to . However, the pre-computation requires the storage of the order and quickly becomes infeasible for higher band-limits [42, 11]. We note that the complexity to compute the scaled associated Legendre function