Iterative Residual Fitting for Spherical Harmonic Transform of Band-Limited Signals on the Sphere: Generalization and Analysis
We present the generalized iterative residual fitting (IRF) for the computation of the spherical harmonic transform (SHT) of band-limited signals on the sphere. The proposed method is based on the partitioning of the subspace of band-limited signals into orthogonal subspaces. There exist sampling schemes on the sphere which support accurate computation of SHT. However, there are applications where samples (or measurements) are not taken over the predefined grid due to nature of the signal and/or acquisition set-up. To support such applications, the proposed IRF method enables accurate computation of SHTs of signals with randomly distributed sufficient number of samples. In order to improve the accuracy of the computation of the SHT, we also present the so-called multi-pass IRF which adds multiple iterative passes to the IRF. We analyse the multi-pass IRF for different sampling schemes and for different size partitions. Furthermore, we conduct numerical experiments to illustrate that the multi-pass IRF allows sufficiently accurate computation of SHTs.
Spherical harmonics, basis functions, spherical harmonic transform, residual fitting, band-limited signals, 2-sphere (unit sphere).
Signals are defined on the sphere in a variety of fields including geodesy , computer graphics , cosmology , astrophysics , medical imaging , acoustics  and wireless communication  to name a few. Spherical harmonic (SH) functions  are a natural choice of basis functions for representing the signal on the sphere in all these applications. Analysis on the sphere is done in both spatial (spherical) and spectral (spherical harmonic) domains. The transformation between the two domains is enabled by the well known spherical harmonic transform (SHT) [8, 9]. For harmonic analysis and signal representation (reconstruction), the ability to accurately compute the SHT of a signal from its samples taken over the sphere is of great importance.
Sampling schemes have been devised in the literature for the accurate and efficient computation of SHTs [10, 11]. However, the samples may not be available, in practice (e.g., [12, 5]), over the grid defined by these sampling schemes. To support the computation of SHTs in applications where samples or data-sets are not available on the pre-defined grid, least squares fitting (LSF) methods have been investigated for efficient computation of the SHTs [12, 13, 14, 15, 16, 17, 18]. LSF methods formulate a large linear system of basis functions and then attempt to solve it efficiently. However, due to memory overflow, it is not suitable for systems with large band-limits, [19, 20, 21]. To solve this problem, an iterative residual fitting (IRF) method has been proposed in  as an extension of LSF and incorporates a divide and conquer technique for the computation of SHTs. The basic idea of IRF is to divide the subspace spanned by all spherical harmonics into smaller partitions and then perform least squares on each partition iteratively. Although IRF is fast, it creates a less accurate reconstruction  as the size of the harmonic basis increases for large band-limits. To improve the reconstruction accuracy, a multi-pass IRF approach is used which includes multiple passes for fitting. This is same as IRF but it involves multiple IRF operations rather than one. A variant of this scheme is presented in , where reconstruction for 3D surfaces is carried out by taking large number of samples.
In this paper, we present an IRF method for the computation of the SHT of a band-limited signal in a general setting that partitions the subspace of band-limited signals into orthogonal subspaces, where each orthogonal subspace can be spanned by different numbers of basis functions. We also formulate multi-pass IRF to improve the accuracy of computation of the SHT. We analyze multipass IRF for different choices of partitioning of the subspace and sampling schemes [10, 11, 22, 19] and show that the computation of the SHT converges in all cases. We also show that the convergence is fast for the partition choice considered in this work.
The remainder of this paper is organized as follows. Section 2 provides the necessary mathematical background and notation required to understand the work. IRF and multi-pass IRF methods are formulated in section 3. In section 4, we carry out accuracy analysis of the proposed IRF method for different partition choices and sampling schemes. Finally, concluding remarks are presented in section 5.
2 Mathematical Background
2.1 Signals on the Sphere
A point on the unit sphere , is parameterized by , where represents the transpose, represents the co-latitude and denotes the longitude. The space of square integrable complex functions of the form , defined on the unit sphere, form a complex separable Hilbert space, denoted by , with inner product defined as by 
where represents the complex conjugate operation. The functions with finite induced norm are referred to as signals on the sphere.
2.2 Spherical Harmonics
Spherical harmonic (SH) functions, denoted by for integer degree and integer order , serve as complete basis for . Due to the completeness of the SH functions, any function on the sphere can be expanded as
where are the SH coefficients of degree and order and form the spectral domain representation of the signal , given by the spherical harmonic transform (SHT) defined as
The signal is band-limited at degree if for all , . A set of band-limited signals forms an dimensional subspace of , denoted by .
3 Generalized Iterative Residual Fitting
3.1 Iterative Residual Fitting (IRF) – Formulation
The IRF method is based on the idea to partition the subspace into smaller subspaces and carry out least-squares estimation on these partitions iteratively. In this way, a large linear problem is divided into manageable small subsets of linear problems. The subspace has graphical representation of the form shown in Fig. 1, which also represents the SH (spectral) domain formed by the SH coefficients of the band-limited signal in . We partition into orthogonal subspaces , each of dimension . We analyse different choices for partitioning later in the paper. We index the SH functions that span the subspace as . We also define .
Given samples (measurements) of the band-limited signal , we wish to compute SH coefficients. By defining a vector
of measurements (samples) of the signal on the sphere and the matrix , with entries , of size containing SH functions that span the subspace evaluated at sampling points, the vector of SH coefficients can be iteratively computed (estimated) in the least-squares sense as
where represents the Hermetian of a matrix and
is the residual between the samples of the signal and the signal obtained by using the coefficients for and the estimation of coefficients is carried out iteratively for . We note that the computational complexity for (5) for each would be of the order of . The computational complexity to compute (6) is . We later analyse the estimation accuracy of the IRF method for different sampling schemes on the sphere and different partitions of the subspace of band-limited signals. For a special case of partitioning the subspace into subspaces based on the degree of spherical harmonics , it has been shown that the iterative residual fitting allows sufficiently accurate estimation of SH coefficients .
The proposed IRF method enables accurate computation of the SHT of signals with a sufficient number of randomly distributed samples. The IRF algorithm finds significance use in applications where samples on the sphere are not taken over a predefined grid. For example, the samples are taken over the cortical surface in medical imaging , where IRF allows sufficient accurate parametric modeling of cortical surfaces.
3.2 Multi-Pass IRF and Residual Formulation
To improve the estimation accuracy, we employ the so-called multi-pass IRF  which is based on the use of IRF method in an iterative manner. In multi-pass IRF, the IRF algorithm is run for a number of iterations, denoted by . To clarify the concept, we incorporate the iteration index in the formulation in (5) and (6) as
After -th iteration, can be computed for each as
the residual after the -th iteration is given by
In general, the residual in (11) depends on the distribution of sampling points and nature of partitioning of . In the next section, we show that the residual converges to zero for a variety of sampling schemes and different partitions.
4 Analysis of Multi-Pass IRF
4.1 Partition Choices
In order to understand the partitioning of , we refer to the graphical representation of shown in Fig. 1, which describes the position of spherical harmonic coefficients with respect to degree and order . We give numbers to the spectral harmonic coefficients (basis functions) shown in Fig. 1 from 1 to in a way that we start the domain from and then traverse the whole domain by to for increasing values of . In a similar way, we can also traverse the whole domain by fixing for all values of . We analyse four different type of partitions, whose sizes vary with the increasing or decreasing values of degrees and orders . The size of each partition is denoted by . In all the partitions, the generalized IRF is run for all values of and for a fixed value of .
Partition Choice 1
We first consider the partitioning of based on the spherical harmonic degree . We take partitions each for degree such that the subspace is spanned by spherical harmonics of degree . Consequently, the dimension of each subspace is . As mentioned earlier, the IRF has been applied already for this choice of partition . We show through numerical experiments that alternative choices for partitioning result in faster convergence and more accurate computation of the SHT.
Partition Choice 2
For partition choice 2, we combine the -th partition choice 1 and -th partition choice 1, to obtain or partitions for even or odd band-limit respectively. For even , each partition 2 is of size for . For odd , we have partitions with for and one partition of size .
Partition Choice 3
Here, we consider partitioning with respect to each order (see Fig. 1). Consequently, we have partitions, one for each order and spanned by SH functions of order .
Partition Choice 4
Partition choice is obtained by combining the partitions in partition choice 3. We obtain partitions by combining partition choice 3 for and for . With such combining, we have partitions of each of size .
Here we analyse the accuracy of the computation of the SHT, that is, the computation of SH coefficients, of the band-limited signal sampled over different sampling schemes. For the distribution of samples on sphere, we consider equiangular sampling  and optimal-dimensionality sampling  in our analysis as these schemes support the accurate computation of the SHT for band-limited signals. Among the sampling schemes on the sphere, which do not support the highly accurate computation of the SHT, we consider the HEALPix sampling scheme  and random samples with uniform distribution with respect to the differential measure .
In order to analyse accuracy, we take a test signal by first generating the spherical harmonic coefficients with real and imaginary part uniformly distributed in and using (2) to obtain the signal over the samples for each sampling scheme. For a meaningful comparison, we take approximately the same number of points for each sampling scheme. We apply the proposed multi-pass IRF for each choice of partition and each sampling scheme to compute the estimate of SH coefficients and record the maximum error between reconstructed and original SH coefficients given by
which is plotted in logarithmic scale in Fig. 2 for band-limit . Different partition choices and different sampling schemes (see caption for number of samples for each sampling scheme) against the number of iterations of the proposed multi-pass IRF are plotted, where it can be observed that 1) the error converges to zero (, double precision) for all partition choices and sampling schemes, and 2) the error converges quickly for partition choice . We also validate the formulation of the residual in (11) by computing after each iteration of the multi-pass IRF. To illustrate the effect of the number of samples on the accuracy of the proposed multi-pass IRF, we have taken , and samples of optimal dimensionality sampling  and plot the error in Fig. 2(d)-(f), where it is evident that the error converges quickly for a greater number of samples. The convergence of the error is in agreement with the formulation of the residual in (11), however, convergence changes with the sampling scheme and nature of the partition of the subspace of band-limited signals. This requires further study and is the subject of future work.
We have presented the generalized iterative residual fitting (IRF) method for the computation of the spherical harmonic transform (SHT) of band-limited signals on the sphere. Proposed IRF is based on partitioning the subspace of band-limited signals into orthogonal spaces. In order to improve the accuracy of the transform, we have also presented a multi-pass IRF scheme and analysed it for different sampling schemes and for four different size partitions. We have performed numerical experiments to show that accurate computation of the SHT is achieved by multi-pass IRF. For different partitions and different sampling distributions, we have analysed the residual (error) and demonstrated the convergence of the residual to zero. Furthermore, it has been demonstrated that the rate of convergence of error depends on the sampling scheme and choice of partition. Rigorous analysis relating the nature of partitions and convergence of the proposed method and application of proposed method in medical imaging, computer graphics and beyond are subjects of future work.
- A. Amirbekyan, V. Michel, and F. J. Simons, “Parametrizing surface wave tomographic models with harmonic spherical splines,” Geophys. J. Int.l, vol. 174, no. 2, pp. 617–628, Jan. 2008.
- R. Ng, R. Ramamoorthi, and P. Hanrahan, “Triple product wavelet integrals for all-frequency relighting,” ACM Trans. Graph., vol. 23, no. 3, pp. 477–487, Aug. 2004.
- Y. Fantaye, C. Baccigalupi, S. Leach, and A. Yadav, “Cmb lensing reconstruction in the presence of diffuse polarized foregrounds,” J. Cosmol. Astropart. Phys., vol. 2012, no. 12, p. 017, July 2012.
- 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, “Seven-year wilkinson microwave anisotropy probe (wmap) observations: Sky maps, systematic errors, and basic results,” ApJS, vol. 192, no. 2, p. 14, Jan. 2011.
- M. Chung, K. Dalton, L. Shen, A. Evans, and R. Davidson, “Weighted fourier series representation and its application to quantifying the amount of gray matter,” IEEE Trans. Med. Imag., vol. 26, no. 4, pp. 566–581, Apr. 2007.
- W. Zhang, M. Zhang, R. A. Kennedy, and T. D. Abhayapala, “On high resolution head-related transfer function measurements: An efficient sampling scheme,” IEEE/ACM Trans. Audio, Speech, Language Process., vol. 20, no. 2, pp. 575–584, Dec. 2012.
- Y. F. Alem, Z. Khalid, and R. A. Kennedy, “3D spatial fading correlation for uniform angle of arrival distribution,” IEEE Commun. Lett., vol. 19, pp. 1073–1076, June 2015.
- R. A. Kennedy and P. Sadeghi, Hilbert Space Methods in Signal Processing. Cambridge, UK: Cambridge University Press, March 2013.
- J. J. Sakurai, Modern Quantum Mechanics, 2nd ed. Reading, MA: Addison Wesley Publishing Company, Inc., 1994.
- J. D. McEwen and Y. Wiaux, “A novel sampling theorem on the sphere,” IEEE Trans. Signal Process., vol. 59, no. 12, pp. 5876–5887, Dec 2011.
- Z. Khalid, R. A. Kennedy, and J. D. McEwen, “An optimal-dimensionality sampling scheme on the sphere with fast spherical harmonic transforms,” IEEE Trans. Signal Process., vol. 62, no. 17, pp. 4597–4610, Sep. 2014.
- N. Sneeuw, “Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective,” Geophys. J. Int., vol. 118, no. 3, pp. 707–716, Sep. 1994.
- D. M. Healy, Jr., D. Rockmore, P. J. Kostelec, and S. S. B. Moore, “FFTs for the 2-sphere - improvements and variations,” J. Fourier Anal. and Appl., vol. 9, pp. 341–385, 2003.
- J. Blais and M. Soofi, “Spherical harmonic transforms using quadratures and least squares,” in Computational Science â ICCS 2006, ser. Lecture Notes in Computer Science. Springer Berlin Heidelberg, 2006, vol. 3993, pp. 48–55.
- K. Ivanov and P. Petrushev, “Irregular sampling of band-limited functions on the sphere,” Appl Comput Harmon Anal., vol. 37, no. 3, pp. 545 – 562, Nov. 2014.
- J. Keiner, S. Kunis, and D. Potts, “Efficient reconstruction of functions on the sphere from scattered data,” J. Fourier Anal. Appl., vol. 13, no. 4, pp. 435–458, May 2007.
- S. Kunis and D. Potts, “Fast spherical fourier algorithms,” J. Comput. Appl. Math., vol. 161, no. 1, pp. 75 – 98, Dec. 2003.
- ——, “Stability results for scattered data interpolation by trigonometric polynomials,” SIAM J. Sci. Comput., vol. 29, no. 4, pp. 1403–1419, Feb. 2007.
- L. Shen and M. K. Chung, “Large-scale modeling of parametric surfaces using spherical harmonics,” in 3D Data Processing, Visualization, and Transmission, Third International Symposium on, June 2006, pp. 294–301.
- C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations. Society for Industrial and Applied Mathematics, 1995.
- R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics, 1994.
- 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, no. 2, p. 759, Apr. 2005.