A Fast and Accurate Algorithm for Spherical Harmonic Analysis on HEALPix Grids with Applications to the Cosmic Microwave Background Radiation

TL;DR

This paper introduces a new fast and accurate spherical harmonic analysis algorithm for HEALPix data, significantly improving computational efficiency and accuracy in analyzing cosmic microwave background radiation.

Contribution

The paper presents a novel spherical harmonic analysis method combining Fourier transforms and Slevinsky's transform, reducing complexity and enhancing accuracy for HEALPix data analysis.

Findings

Computational complexity is reduced to O(N log^2 N).

The method achieves at least twice the convergence rate of existing techniques.

Numerical experiments show improved accuracy in angular power spectrum analysis.

Abstract

The Hierarchical Equal Area isoLatitude Pixelation (HEALPix) scheme is used extensively in astrophysics for data collection and analysis on the sphere. The scheme was originally designed for studying the Cosmic Microwave Background (CMB) radiation, which represents the first light to travel during the early stages of the universe's development and gives the strongest evidence for the Big Bang theory to date. Refined analysis of the CMB angular power spectrum can lead to revolutionary developments in understanding the nature of dark matter and dark energy. In this paper, we present a new method for performing spherical harmonic analysis for HEALPix data, which is a central component to computing and analyzing the angular power spectrum of the massive CMB data sets. The method uses a novel combination of a non-uniform fast Fourier transform, the double Fourier sphere method, and Slevinsky's fast spherical harmonic transform (Slevinsky, 2019). For a HEALPix grid with $N$ pixels (points), the computational complexity of the method is $\mathcal{O}(N\log^2 N)$, with an initial set-up cost of $\mathcal{O}(N^{3/2}\log N)$. This compares favorably with $\mathcal{O}(N^{3/2})$ runtime complexity of the current methods available in the HEALPix software when multiple maps need to be analyzed at the same time. Using numerical experiments, we demonstrate that the new method also appears to provide better accuracy over the entire angular power spectrum of synthetic data when compared to the current methods, with a convergence rate at least two times higher.

Figures

Tables

Table 1. Table 1: Spectral radius of the Richardson iteration matrix from ( 6 ) for different values of N 𝑁 N .

$N_{side}$	$N$	$\rho \left(𝐈-\frac{4\pi }{N}{𝐒}^{\ast }𝐒\right)$
2	48	0.1986
4	192	0.0932
8	768	0.0600
16	3072	0.0475
32	12288	0.0421

$\ell$	176	190	191	230	248	283	292	303	326	366	388	404	421	446	448
$m$	56	81	124	40	155	274	27	145	55	343	200	78	420	284	234

$\ell$	589	633	636	766	829	943	974	1009	1085	1219	1294	1346	1404	1485	1493
$m$	188	269	414	134	517	912	93	483	183	1143	667	259	1400	946	779

Equations

$$C_{\ell }=\frac{1}{2\ell +1}\sum _m|a_{\ell }^m{|}^2.$$

$$\begin{array}{rl}\mathcal{N}& :=\left\{\left(\arccos \left(1-\frac{j^2}{3N_{side}^2}\right),\frac{\pi \left(k+\frac{1}{2}\right)}{2j}\right):j=1,\dots ,N_{side}-1,k=0,\dots ,4j-1\right\}\\ \mathcal{E}& :=\left\{\left(\arccos \left(\frac{2(2N_{side}-j)}{3N_{side}}\right),\frac{\pi \left(k+\frac{(j+1)\text{mod}2}{2}\right)}{2N_{side}}\right):j=N_{side},\dots ,3N_{side},k=0,\dots ,4N_{side}-1\right\}\\ \mathcal{S}& :=\left\{\left(\arccos \left(-\left(1-\frac{j^2}{3N_{side}^2}\right)\right),\frac{\pi \left(k+\frac{1}{2}\right)}{2j}\right):j=1,\dots ,N_{side}-1,k=0,\dots ,4j-1\right\}.\end{array}$$

$$a_{\ell }^m=\frac{4\pi }{N}\sum _{i=1}^N\overline{Y}{}_{\ell }^m({\lambda }_i,{\theta }_i)f({\lambda }_i,{\theta }_i),0\le \ell \le {\ell }_{max},-\ell \le m\le \ell ,$$

$$\hat{f}({\lambda }_i,{\theta }_i)=\sum _{\ell =0}^{{\ell }_{max}}\sum _{m=-\ell }^{\ell }{a}_{\ell }^m{Y}_{\ell }^m({\lambda }_i,{\theta }_i),i=1,\dots ,N$$

$$\text{Analysis:}\mathbf{a}$$

$$\text{Synthesis:}\hat{\mathbf{f}}$$

$$\begin{array}{rl}{\mathbf{r}}^{(k+1)}& =\mathbf{f}-\mathbf{S}{\mathbf{a}}^{(k)},\\ {\mathbf{a}}^{(k+1)}& ={\mathbf{a}}^{(k)}+\mathbf{A}{\mathbf{r}}^{(k+1)},\end{array}$$

$${\mathbf{a}}^{(k+1)}={\mathbf{a}}^{(k)}+\frac{4\pi }{N}{\mathbf{S}}^{\ast }(\mathbf{f}-\mathbf{S}{\mathbf{a}}^{(k)})=\frac{4\pi }{N}{\mathbf{S}}^{\ast }\mathbf{f}+\left(\mathbf{I}-\frac{4\pi }{N}{\mathbf{S}}^{\ast }\mathbf{S}\right){\mathbf{a}}^{(k)}.$$

$$a_{\ell }^m=\sum _{i=1}^N{w}_i\overline{Y}{}_{\ell }^m({\lambda }_i,{\theta }_i)f({\lambda }_i,{\theta }_i),0\le \ell \le {\ell }_{max},-\ell \le m\le \ell ,$$

$${\lambda }_k=\frac{k}{2N_{side}}\pi ,k=0,\dots ,4N_{side}-1.$$

$$\text{f(λ,θj)=:}{f}_j(\lambda )=\sum _{n=-2j}^{2j}{c}_n^{(j)}{e}^{\mathrm{i}n\lambda },$$

$$f\left(\frac{k+\frac{1}{2}}{2j}\pi ,{\theta }_j\right),k=0,\dots ,4j-1.$$

$$f_j\left(\frac{k+\frac{1}{2}}{2j}\pi \right)=\sum _{n=-2j}^{2j-1}{c}_n^{(j)}{e}^{\mathrm{i}n\frac{k+\frac{1}{2}}{2j}\pi }=\sum _{n=-2j}^{2j-1}{c}_n^{(j)}{e}^{\mathrm{i}n\frac{\pi }{4j}}{e}^{\mathrm{i}n\frac{k}{2j}\pi }=\sum _{n=-2j}^{2j-1}{\tilde{c}}_n^{(j)}{e}^{\mathrm{i}n\frac{k}{2j}\pi },k=0,\dots ,4j-1,$$

$$\sum _{n=-2j}^{2j-1}{\tilde{c}}_n^{(j)}{e}^{\mathrm{i}n\frac{k}{2j}\pi }=f\left(\frac{k+\frac{1}{2}}{2j}\pi ,{\theta }_j\right),k=0,\dots ,4j-1,$$

$$x(\lambda ,\theta )=\cos \lambda \sin \theta ,y(\lambda ,\theta )=\sin \lambda \sin \theta z(\lambda ,\theta )=\cos \theta ,(\lambda ,\theta )\in [0,2\pi ]\times [0,\pi ],$$

$$\tilde{f}(\lambda ,\theta )=\{\begin{array}{ll}g(\lambda ,\theta ),& (\lambda ,\theta )\in [0,\pi ]\times [0,\pi ],\\ h(\lambda -\pi ,\theta ),& (\lambda ,\theta )\in [\pi ,2\pi ]\times [0,\pi ],\\ h(\lambda ,2\pi -\theta ),& (\lambda ,\theta )\in [0,\pi ]\times [\pi ,2\pi ],\\ g(\lambda -\pi ,2\pi -\theta ),& (\lambda ,\theta )\in [\pi ,2\pi ]\times [\pi ,2\pi ],\end{array}$$

$$\tilde{f}(\lambda ,\theta )\approx \sum _{j=-\lfloor \frac{m}{2}\rfloor }^{\lceil \frac{m}{2}\rceil -1}\sum _{k=-\lfloor \frac{n}{2}\rfloor }^{\lceil \frac{n}{2}\rceil -1}{C}_{jk}{e}^{\mathrm{i}j\theta }{e}^{\mathrm{i}k\lambda },$$

$$f_j=\sum _{k=0}^{n-1}{c}_k{e}^{-2\pi i{x}_j{\omega }_k},0\le j\le n-1,$$

$$\mathbf{f}={\hat{F}}_2\mathbf{c},$$

$${\hat{F}}_2^{\ast }{\hat{F}}_2\mathbf{c}={\hat{F}}_2^{\ast }\mathbf{f},$$

$$\left|x_j-\frac{j}{n}\right|\le \frac{\gamma }{n},0\le j\le n-1,$$

$${\hat{F}}_2^{\ast }{\hat{F}}_2{\mathbf{e}}_1={\hat{F}}_2^{\ast }1=(1^T{\hat{F}}_2{)}^{\ast }=({\hat{F}}_2^{T}1{)}^{\ast }.$$

$$\widetilde{f}(\lambda,\theta)=\sum_{j=0}^{p}\sum_{k=-p}^{p}g_{j}^{k}\frac{e^{{\rm i}k\lambda}}{\sqrt{2\pi}}\left\{\begin{array}[]{c}\cos j\theta,\quad\quad\;k\,\text{even}\\ \sin(j+1)\theta,\,k\,\text{odd}\end{array}\right\}.$$

$$\tilde{f}(\lambda ,\theta )=\sum _{j=-p}^p\sum _{k=-p}^{p-1}{C}_{jk}{e}^{\mathrm{i}j\theta }{e}^{\mathrm{i}k\lambda },$$

$$X_{j,k}=\{\begin{array}{ll}{C}_{j,k}& \text{if −p+1≤k≤p−1}\\ \frac{1}{2}{C}_{j,-p}& \text{if k=±p}\end{array},-p\le j,k\le p.$$

$$\tilde{f}(\lambda ,\theta )$$

$$\displaystyle=\sum_{j=0}^{p}\sum_{k=-p}^{p}e^{{\rm i}k\lambda}\left\{\begin{array}[]{c}((X_{jk}+\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu_{-jk})\cos(j\theta),\;\text{$k$ even}\\ ((X_{jk}-\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu_{-jk})\sin(j\theta),\;\text{$k$ odd}\end{array}\right\},$$

$$\left(\begin{array}{cccccc}{g}_0^0& g_0^{-1}& g_0^1& \cdots & g_0^{-p}& g_0^p\\ g_1^0& g_1^{-1}& g_1^1& \cdots & g_1^{-p}& g_1^p\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ g_{p-1}^0& g_{p-1}^{-1}& g_{p-1}^1& \cdots & g_{p-1}^{-p}& g_{p-1}^p\\ g_p^0& 0& 0& \cdots & g_p^{-p}& g_p^p\end{array}\right).$$

$$\left(\begin{array}{cccccccc}{a}_0^0& a_1^{-1}& a_1^1& a_2^{-2}& a_2^2& \cdots & a_p^{-p}& a_p^p\\ a_1^0& a_2^{-1}& a_2^1& a_3^{-2}& a_3^2& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{p-2}^0& a_{p-1}^{-1}& a_{p-1}^1& a_p^{-2}& a_p^2& & ⋮& ⋮\\ a_{p-1}^0& a_p^{-1}& a_p^1& 0& 0& \cdots & 0& 0\\ a_p^0& 0& 0& 0& 0& \cdots & 0& 0\end{array}\right).$$

$$f(\lambda ,\theta )=\sum _{j=1}^3{c}_j(2-2\mathbf{x}(\lambda ,\theta )\cdot \mathbf{x}({\lambda }_j,{\theta }_j){)}^{3/2},$$

$$\{c_1,c_2,c_3\}=\{5,-3,8\},$$

Full text

Keywords: HEALPix, Fast Spherical Harmonic Transform, Non-uniform Fast Fourier Transform, Double Fourier Sphere, Cosmic Microwave Background Radiation

1 Introduction

About 379,000 years after the universe began, the dense plasma of matter cooled enough for neutral hydrogen to form. During this epoch of recombination, the universe was becoming increasingly transparent to photons, which eventually began to move freely through space. Now faintly glowing at the edge of the observable universe, these photons form the Cosmic Microwave Background (CMB) radiation, which has become the strongest evidence for the Big Bang Theory to date [3]. While the CMB has been deemed “the most perfect black body ever measured in nature” [42], there are temperature and polarization fluctuations that give insight into the primordial universe [28]. These anisotropies are consequences of the initial density distribution of matter, and analyzing them can provide a better understanding of the geometry and composition of the universe [3, 15].

Using ground-based telescopes, balloons, and satellites which probe the sky in the microwave and infra-red frequencies, scientists have measured the CMB temperature differences at small angular scales. These measurements are quantized and stored as a high resolution sky map of the CMB using the Hierarchical Equal Area isoLatitude Pixelation (HEALPix) scheme [11] for the sphere; see Figure 1a) for an example sky map. Once a sky map is composed, it can then be analyzed by its angular power spectrum. This quantity measures the amplitude of the CMB temperature fluctuations as a function of angular scale and is used to estimate parameters of the cosmological model for the universe [42]. For example, the confirmation of the first peak in the temperature angular power spectrum affirmed that the universe is spatially flat [17]. The values of the temperature angular spectrum at higher frequencies are crucial to many aspects of modern cosmology, including the density of dark matter and dark energy in the universe. The CMB power spectrum (Figure 1b) is calculated from the spherical harmonic coefficients, $a_{\ell}^{m}$, of the sky map as follows:

[TABLE]

The spherical harmonic conventions used in this work are detailed in Appendix A.

The HEALPix scheme [11] and the associated eponymous software [10] have a number of desirable properties for data collection on the sphere. First, each pixel has the same surface area and the pixel centers (points) are quasi-uniformly distributed over the sphere. This is important since any white noise produced by the microwave receivers is exactly integrated into white noise in the pixel area. Second, the pixels produced by the scheme are based on a hierarchical subdivision of the sphere, which allows for data locality in computer memory and fast search procedures. Finally, the pixel centers are isolatitudinal, allowing for a significant reduction in the computational cost of performing discrete spherical harmonic transforms—a central component to computing and analyzing the angular power spectrum of the CMB data sets, which from the Planck mission consist of millions of pixels [30]. These properties have made the HEALPix scheme popular for other applications in astrophysics/astronomy [35, 21, 29], and to several other disciplines, including geophysics [41], planetary science [25], nuclear engineering [32], and computer vision [16].

In this paper, we focus on an aspect of the HEALPix scheme that has received very little attention in the literature: accuracy and computational complexity improvements of the discrete spherical harmonic transform. We first review the current techniques used in the HEALPix software [10], which are based on equal-weight quadrature, ring-weight quadrature, and pixel-weight quadrature. We then introduce a new algorithm for computing spherical harmonic coefficients for data collected on HEALPix grids. The main motivation for the method is Slevinsky’s recently developed fast spherical harmonic transform (FSHT) [38], which converts bivariate Fourier coefficients for data on the sphere to spherical harmonic coefficients of the data with near optimal complexity. By combining the nonuniform fast Fourier transform (NUFFT) [33] and the double Fourier sphere (DFS) [40] methods, we give a fast and accurate method for obtaining the bivariate Fourier coefficients for functions sampled on the HEALPix grid, which we then use with the FSHT to obtain the spherical harmonic coefficients. For a HEALPix grid with $N$ pixels (points), the computational complexity of the method is $\mathcal{O}(N\log^{2}N)$, with an initial set-up cost of $\mathcal{O}(N^{3/2}\log N)$, which compares favorably with the complexity of the current methods available in the HEALPix software when multiple maps need to be analyzed at the same time. Using numerical experiments, we demonstrate that the new method also appears to be more accurate than the current methods for synthetic data over the whole spectrum, with a convergence rate at least two times higher. We believe this new scheme will be useful not only for CMB analysis, but also for the many applications of the HEALPix scheme given above that require a spherical harmonic analysis. Additionally, the algorithm presented here has natural generalizations for other “equal-area” isolatitudinal sampling strategies for sphere that do not have a natural way to do fast and accurate spherical harmonic transforms [34, 6, 20, 22].

The remainder of the paper is structured in the following manner. In section 2, we offer supporting information on the HEALPix grid as well as details and analysis of the current methods used in the HEALPix software for computing the spherical harmonic coefficients of CMB maps. We present the new algorithm for fast spherical harmonic analysis of data collected on the HEALPix grid in section 3. Numerical results comparing the presented method with that of the methods in the HEALPix software for calculating the angular power spectrum of functions on the sphere are given in section 4. Finally, in section 5, we give some brief conclusions and remarks on future directions of the work.

2 Background and Current Approach

2.1 HEALPix Scheme

The HEALPix scheme333The HEALPix scheme produces a grid consisting of a collection of pixels of different shapes but the same area. However, for our method we do not exploit this fact and simply treat the center of each pixel as a point with the given value of the pixel. was created to discretize functions on the sphere at high resolutions. In addition to creating an equal area pixelization of the sphere, one of the primary motivations behind the scheme was to allow for more computationally efficient spherical harmonic transforms on increasingly large CMB datasets [11]. While there are many options for discretizing the sphere, there is no known deterministic method that gives a set of quasiuniform points and allows for exact spherical harmonic decompositions of band-limited functions using equal-weight quadrature. While the HEALPix scheme does not offer optimal complexity for spherical harmonic analyses, it does achieve some efficiency gains over existing schemes for discretizing the sphere. This improvement is accomplished primarily by the isolatitudinal distribution of pixels.

The HEALPix grid resolution is defined using the parameter $N_{side}=2^{t},\,t\in\mathbb{N}$, which creates $N_{side}^{2}$ equal area divisions of each base pixel. Figure 2 illustrates the base resolution grid, $t=0$, and the increasing levels of refinement $t=1,2,3$, where each base pixel is subdivided further into four equal area pixels. A HEALPix map therefore has $N=12N_{side}^{2}$ equal area (but differently shaped) pixels, with the centers placed on $4N_{side}-1$ rings of constant latitude. For any $N_{side}$, the HEALPix centers, which we henceforth call the HEALPix points, are defined algebraically using three regions of the sphere, two polar ($\mathcal{N}$ and $\mathcal{S}$) and one equatorial ($\mathcal{E}$) [19]. In spherical coordinates, the points in these regions are given as

[TABLE]

The final HEALPix point set is $\mathcal{X}=\mathcal{N}\bigcup\nolimits\mathcal{E}\bigcup\nolimits\mathcal{S}$. The number of points on each ring varies in the polar regions, with only four points on the rings closest to the north and south poles of the sphere, whereas the rings in the equatorial region have the same number of points. This point distribution is illustrated more clearly in Figure 3, where the HEALPix points are displayed to a latitude-longitude grid.

The biggest computational advantage for spherical harmonic analysis in the HEALPix scheme lies in the equally-spaced points on each ring of constant latitude. While this aides computation in the longitude direction with FFTs, the misaligned and unequally spaced points in latitude make fast bivariate Fourier analysis impossible without modification. We address this in the new algorithm presented in section 3.

2.2 HEALPix Software Spherical Harmonic Analysis

The standard method in the HEALPix software [10] for estimating the angular power spectrum (1) of data at the HEALPix points approximates the exact spherical harmonic coefficients ($\widetilde{a}_{\ell}^{m}$) of the data as

[TABLE]

where $(\lambda_{i},\theta_{i})$ are HEALPix points in longitude-latitude, $f$ is the data, and $Y_{\ell}^{m}$ is a spherical harmonic of degree $\ell$ and order $m$ (see Appendix A for a discussion of the spherical harmonic conventions used in this paper). While the user can input any band limit $\ell_{max}$ for this approximation, the software default is $\ell_{max}=3N_{side}-1$. Due to the isolatitudinal nature of the HEALPix points, this computation is done with $\mathcal{O}(N^{3/2})$ complexity as opposed to $\mathcal{O}(N^{2})$ [11]. Note that $N=\mathcal{O}(\ell_{max}^{2})$, so the complexity of the $a_{\ell}^{m}$ computation is equivalent to $\mathcal{O}(\ell_{max}^{3})$. The expression (3) is a low-order approximation to the continuous inner product (23) which defines the coefficients, since it uses a simple equal weight quadrature. To improve this approximation, the software employs an iterative procedure, which is referred to as a “Jacobi iteration” [11]. In order to illustrate how the iterative method converges, we explain it below in the language of linear algebra.

The analysis operation, defined in (3), produces an approximation to the spherical harmonic coefficients from the data $f$ on the sphere, whereas the synthesis operation reconstructs the data given the spherical harmonic coefficients:

[TABLE]

Note that we use a hat on $f$ to indicate that computing the spherical harmonic coefficients according to (3) and using them in (4) gives different function values in general. In matrix-vector notation, we denote (3) and (4) as

[TABLE]

where $\mathbf{a}$ is the vector of spherical harmonic coefficients and $\mathbf{f}$ and $\widehat{\mathbf{f}}$ are the vectors of data values at the HEALPix points. Using this notation, the iterative procedure in the HEALPix software can be written as

[TABLE]

where $\mathbf{r}$ is the residual vector and $\mathbf{a}^{(0)}=\mathbf{A}\mathbf{f}$. Substituting the first equation of (5) into the last and using the fact that the analysis matrix satisfies $\mathbf{A}=\frac{4\pi}{N}\mathbf{S}^{*}$, gives the iteration

[TABLE]

This is just stationary Richardson iteration (or Gradient Decent) with relaxation parameter $\frac{4\pi}{N}$ applied to the normal equations $\mathbf{S}^{*}\mathbf{S}\mathbf{a}=\mathbf{S}^{*}\mathbf{f}$ [4, pp. 44–45]. Thus, the iterative procedure converges to the least squares solution to (3), provided the spectral radius of $\mathbf{I}-\frac{4\pi}{N}\mathbf{S}^{*}\mathbf{S}$ is less than one. The spectral radius also determines the convergence rate. Since the HEALPIx points are equidistributed, we know that (3) converges to the integral (23) as $N\rightarrow\infty$ (in exact arithmetic) [14]. Thus, the spectral radius of $\mathbf{I}-\frac{4\pi}{N}\mathbf{S}^{*}\mathbf{S}$ converges to [math] as $N\rightarrow\infty$ and we expect the iteration (6) to converge more rapidly as $N$ increases. Table 1 gives evidence of this result by displaying the spectral radius of $\mathbf{I}-\frac{4\pi}{N}\mathbf{S}^{*}\mathbf{S}$ for increasing values of $N$.

The default option in the HEALPix software sets the number of iterations of (6) to 3. While this does improve the accuracy of computing the spherical harmonic coefficients, it adds to the cost, as each iteration requires doing an analysis and synthesis ((3) and (4)) at a cost of $\mathcal{O}(\ell_{max}^{3})$ operations each. Since the solution converges to the least squares solution, one could improve the convergence of the Richardson iteration method by using algorithms like LSQR or conjugate gradient on the normal equations [27].

2.2.1 Pixel Weights and Ring Weights

As an alternative to the iterative scheme, the HEALPix software also has the option of using quadrature weights to improve the accuracy of the computation of the spherical harmonic coefficients. In this case, the equal weight quadrature approximation (3) is generalized to

[TABLE]

where $w_{i}$ are the quadrature weights. There are two options for using quadrature weights. The first is “pixel weights”, which uses different weights for each HEALPix point. These weights are computed using the positive quadrature weight algorithm from [19], which consists of solving a system involving a Gram matrix containing the spherical harmonics whose size is proportional to $N$ [31]. For large $N$, the weights are computed once and stored. The second option is to use “ring weights”, which use different weights for each ring of the HEALPix point sets. The computation of the ring weights is done using similar ideas to the pixel weights, but a much smaller system has to be solved [31]. The new method introduced in this paper does not use quadrature weights directly, but instead computes the bivariate Fourier coefficients of the HEALPix data and then uses these to obtain the spherical harmonic coefficients.

3 HP2SPH

The algorithm presented here, named HP2SPH, introduces a new way to calculate the spherical harmonic coefficients of data sampled at the HEALPix points (2). The outline for the algorithm is given in Algorithm 1, and each of the pieces are described below.

3.1 Step 1: Transform the data to a tensor product latitude-longitude grid

As described in Section 2.1, the HEALPix grid has an unequal number of points on the rings in the northern ($\mathcal{N}$) and southern ($\mathcal{S}$) sets (2), and the points on the rings in the equatorial ($\mathcal{E}$) set are shifted on every other ring. This structure leads to the pixels being misaligned in latitude. By upsampling the data on the northern and southern points in longitude so that there are an equivalent samples of the data on each ring and shifting the data at equatorial points in longitude, we can use fast algorithms to obtain the bivariate Fourier coefficients of the data as discussed in the next section. On the two polar point sets, we upsample the data using the trigonometric interpolant of the data on each ring of these sets to the non-shifted equally spaced longitude points on the equatorial rings, i.e.,

[TABLE]

We also interpolate the data on the rings in the equatorial point set with shifted longitude points, to these $\lambda$ values. Figure 4(b) illustrates the upsampling procedure leading to a tensor product latitude-longitude grid of data of size $(4N_{side}-1)\times 4N_{side}$.

We describe the interpolation procedure here for the data in the northern point set $\mathcal{N}$. Consider the latitude values for the northern rings, $\theta_{j}=\arccos\left(1-\frac{j^{2}}{3N_{side}^{2}}\right),\;j=1,\dots,N_{side}$. We approximate the data in each ring using a trigonometric expansion of the form

[TABLE]

The coefficients in the expansion are determined by enforcing interpolation of the given data values

[TABLE]

With the minor algebraic manipulation of (9),

[TABLE]

we see the interpolation conditions yield the system

[TABLE]

which can be computed using the inverse FFT. We can then obtain the Fourier coefficients $c_{n}^{(j)}$ in (9) for the data at the non-shifted values through simple multiplication444Horner’s rule (and Estrin’s scheme for higher accuracy at small angles [7]) could also be used to implement the shift, avoiding loss of accuracy due to evaluation of high frequency complex exponentials.. Finally, we pad the vector containing the coefficients $c_{n}^{(j)}$ with the appropriate number of zeros to get to a total of $4N_{side}$, so that the expansion in longitude in each ring has the same number of Fourier coefficients. The values of the interpolant on each ring can then be obtained at the upsampled values (8) by applying the FFT on these padded vectors. A similar procedure can be applied to the data on the rings in the southern point set $\mathcal{S}$.

On the rings in the equatorial set $\mathcal{E}$ where the longitude values are shifted by $\pi(k+\frac{1}{2})/(2N_{side})$, we compute the Fourier coefficients of the data using (10) with $j=N_{side}$. We then obtain the coefficients in (9) at the non-shifted points from which the values can be computed using the FFT. No padding or upsampling is needed in this case.

3.2 Step 2: Compute Bivariate Fourier Coefficients

Bivariate Fourier analysis for data on the sphere requires the application of the DFS method to obtain periodicity of the data in latitude and to retain spherical symmetry. When we apply this method to the upsampled HEALPix data, there is an issue that the points in latitude are not equally-spaced, making the standard FFT unsuitable. To bypass this issue we use an NUFFT. Both the DFS technique and NUFFT method we use are discussed below for completeness.

3.2.1 Double Fourier Sphere (DFS) Method

A natural approach to representing a function on the sphere is to use a latitude-longitude coordinate transform, defined by

[TABLE]

which maps the sphere to a rectangular domain. While this transformation allows for performing computations with the function $f(\lambda,\theta)=f(x(\lambda,\theta),y(\lambda,\theta),z(\lambda,\theta))$, it also introduces artificial boundaries at the north and south poles. In addition, the change of variables does not maintain the symmetry of functions on the sphere. Specifically, the transform described in (11) does not preserve the periodicity in the latitude direction, which is necessary for bivariate Fourier analysis to be applicable and for results to make physical sense. These problems are solved by the DFS method.

Originally introduced by Merilees in [23] (see also [40]) the DFS method transforms a function on the sphere to a rectangular grid while maintaining bi-periodicity. This can be thought of as “doubling up” the function $f(\lambda,\theta)$ to form a new function that preserves periodicity in both the latitude and longitude directions. Algebraically, this new function, $\widetilde{f}(\lambda,\theta)$, is defined on $[0,2\pi]\times[0,2\pi]$ as follows [40]

[TABLE]

where $g(\lambda,\theta)=f(\lambda,\theta)$ and $h(\lambda,\theta)=f(\lambda+\pi,\theta)$ for $(\lambda,\theta)\in[0,\pi]\times[0,\pi]$. Figure 5 illustrates the DFS method applied to the surface of the Earth, which shows the preservation of bi-periodicity in (c). We note that the DFS method can also be easily applied to discrete data sampled at a tensor product latitude-longitude grid using (12), which is what we do for the upsampled HEALPix data. In this case, (12) corresponds to flipping and shifting the data matrix appropriately.

Once the DFS method is applied to a function on the sphere, it can be approximated using a 2D bivariate Fourier expansion:

[TABLE]

where $m$ and $n$ represent the number of frequencies in (doubled-up) latitude and longitude, respectively.

Note that the HEALPix grid does not include points at the north and south poles. When applying the DFS to the upsampled data from Step 1, this leads to a relatively large gap in the points in latitude over the poles compared to the other points, which can lead to issues with the inverse NUFFT (see below). To bypass this issue, we construct values at the two poles by using a weighted quadratic least squares fit [8] that combines the data from the three rings closest to each pole.

Remark**.**

The Fourier coefficients of the upsampled data in longitude are computed in Step 1. These can be used directly in the DFS procedure by applying (12) in Fourier space in the $\lambda$ variable, which amounts to appending the (padded) coefficient matrix from Step 1 with a flipped version of itself with all odd wave numbers multiplied by $-1$. It then only remains to compute the Fourier coefficients in latitude $\theta$ to obtain the full bivariate Fourier coefficients. This is the focus of the next step.

3.2.2 Nonuniform Fast Fourier Transform (NUFFT)

The use of the nonuniform discrete Fourier transform (NUDFT) in many domain sciences has led to the development of algorithms for computing it efficiently. If these algorithms are quasi-optimal requiring $\mathcal{O}(n\log n)$, then they are referred to as a nonuniform fast Fourier transform (NUFFT). Given a vector $\mathbf{c}\in\mathbb{C}^{n}$, the one-dimensional NUDFT computes the vector $\mathbf{f}\in\mathbb{C}^{n}$ defined by

[TABLE]

where $x_{j}\in[0,1]$ are the samples and $\omega_{k}\in[0,n]$ are the frequencies. If the samples are equispaced ($x_{j}=j/n$) and the frequencies are integer ($\omega_{k}=k$), then the the transform is a uniform DFT, which can be computed by the FFT in $\mathcal{O}(n\log n)$ operations [5]. When either the samples are nonequispaced or the frequencies are noninteger, the FFT does not directly apply without some careful manipulation [2].

Ruiz and Townsend [33] contributed to the collection of NUFFT algorithms by utilizing low rank matrix approximations to relate the NUDFT to the uniform DFT. There are three types of NUDFTs and inverse NUDFTs that they account for in their algorithm: NUDFT-I, which has uniform samples but noninteger frequencies; NUDFT-II, which has nonuniform samples and integer frequencies; NUDFT-III, which has both nonuniform samples and nonuniform frequencies [12]. Since our HP2SPH method only uses the one-dimensional inverse NUFFT of the second type, we discuss the NUFFT-II algorithm [33].

Given Fourier coefficients, $\mathbf{c}\in\mathbb{C}^{n\times 1}$, the NUFFT-II attempts to approximate the vector

[TABLE]

to machine precision in quasi-optimal complexity. Here $(\widehat{F}_{2})_{jk}=e^{-2\pi ix_{j}k}$, $x_{j}$ are nonuniform samples (restricted to be in $[0,1]$), and $k$ are integer frequencies for $0\leq j,k\leq n-1$. Notice that the DFT matrix for uniform samples and integer frequencies is similarly $F_{jk}=e^{-2\pi ijk/n}$. The key to the NUFFT-II algorithm described in [33] is that if the samples are nearly equispaced, then $\widehat{F}_{2}$ can be related to the Hadamard product of $F$ and a low rank matrix. This means that given a rank $K$ approximation which relates $\widehat{F}_{2}$ and $F$, the NUFFT-II can then be computed using $K$ FFTs with $\mathcal{O}(Kn\log n)$ cost. In practice, machine (double) precision can be achieved with $K=14$ [33].

In the case of the inverse NUFFT-II, Ruiz and Townsend advocate for the use of the conjugate gradient (CG) method in order to solve the linear system $\widehat{F}_{2}\mathbf{c}=\mathbf{f}$ for $\mathbf{c}$. Since $\widehat{F}_{2}$ is not hermitian, the CG method is applied to the normal equations:

[TABLE]

where $\widehat{F}_{2}^{*}\widehat{F}_{2}$ is simply a Toeplitz matrix. Therefore, the inverse NUFFT-II can be calculated using the CG method and a fast Toeplitz multiplication [9] in $\mathcal{O}(R_{\scriptscriptstyle CG}n\log n)$ operations, where $R_{\scriptscriptstyle CG}$ is the number of CG iterations. The following suggestion is placed on the nonuniform function samples to avoid ill-conditioning in the system (16) [33]:

[TABLE]

where $0\leq\gamma<1/4$. When this condition is satisfied, it has been experimentally observed that $R_{\scriptscriptstyle CG}\leq 10$ for a large range of $n$.

For the method proposed in this paper, we apply the inverse NUFFT-II in latitude to the DFS upsampled HEALPix data from Step 2. Unfortunately, the HEALPix points in latitude direction do not meet the condition (17). To bypass this issue, we instead use a least squares solution to compute fewer coefficients at the higher wave numbers than what the given data may support. We describe this procedure below since it not discussed in [33].

The inverse NUFFT-II method computes first column of the symmetric Toeplitz matrix $\widehat{F}_{2}^{*}\widehat{F}_{2}$ in (16) in the following manner:

[TABLE]

The last expression above can be computed efficiently by the NUFFT-I algorithm, since the NUDFT-I matrix is simply the transpose of the NUDFT-II matrix [33]. To compute a least squares solution to (15) with fewer coefficients, we simply truncate the first column obtained from the above method to $m<n$ terms and form the resulting $m\times m$ Toeplitz matrix $\widehat{F}_{2}^{*}\widehat{F}_{2}$. The right hand side for the least squares solution is obtained by similarly computing $\widehat{F}_{2}^{*}\mathbf{f}$ and truncating this to $m$ terms.

For the DFS upsampled HEALPix data from Step 2, there are $8N_{side}$ coefficients in latitude, but only $4N_{side}$ coefficients in longitude. To keep the number of Fourier modes in both directions (nearly) the same, we choose $m=4N_{side}+1$ as the truncation parameter for the least squares solution for computing the Fourier coefficients in latitude. This is also a convenient choice since the method in step three for converting bivariate Fourier coefficients of data on the sphere to spherical harmonic coefficients requires the number of coefficients in each direction is the same and an odd number (we explain how to convert the coefficients in longitude to $m=4N_{side}+1$ in the next section).

Remark**.**

For problems where the data may contain noise (e.g., for the CMB application), there could be an issue with this noise being amplified in steps 1 and 2. For step 1, we should not expect any additional noise to be introduced, since we are simply computing the Fourier coefficients on each ring using the original data and then shifting the coefficients and padding them with zeros. Step 2 has two areas where there could be an issue with noisy data. The first is in constructing values at the poles and the second is in the application of the NUFFT in latitude. However, both of these steps apply a least squares procedure, which provides some smoothing. In our tests on CMB data, we did not observe any amplification of noise that was present in the data. The HP2SPH method has a further benefit of using a backward stable algorithm for computing the spherical harmonic coefficients (as discussed next), which ensures that the resulting uncertainty in the spherical harmonic coefficients has only a low algebraic growth with respect to degree and is always proportional to the norm of the noise in the data.

3.3 Step 3: Obtain the spherical harmonic coefficients via the fast spherical harmonic transform (FSHT)

In [38], Slevinsky derives a fast, backward stable method for the transformation between spherical harmonic expansions and their bivariate Fourier series (given in (13)) by viewing it as a change of basis. This relation is defined as a connection problem, and the matrices that arise in the present case are well-conditioned, making them ideal for fast computations. Slevinsky describes the change of basis in two steps: converting expansions in normalized associated Legendre functions to those of only order zero and one, and then re-expressing these in trigonometric form. In other words, it uses spherical harmonic expansions of order zero and one as intermediate expressions between higher-order spherical harmonics and their corresponding bivariate Fourier coefficients.

The first step of the algorithm takes advantage of the fact that the matrix of connection coefficients between the associated Legendre functions of all orders and those of order zero and one can be represented by a product of Givens rotation matrices. This enables the use of the butterfly algorithm, which can be thought of as an abstraction of the algebraic properties of fast Fourier transforms. The term butterfly was introduced in [24], where it was used for analyzing scattering from electrically large surfaces, and then further developed in [26] for use in special function transforms. Slevinsky uses the butterfly algorithm to perform a factorization of the well-conditioned matrices of connection coefficients.

The second step of this method exploits the hierarchical decompositions of the connection coefficient matrices between the associated Legendre functions of order zero and one to the Chebyshev polynomials of the first and second kind, respectively. This step quickly computes the fast orthogonal polynomial transforms using an adaptation of the Fast Multipole Method [13] and low-rank matrix approximations. The total pre-computation time for both steps is $\mathcal{O}(\ell_{max}^{3}\log\ell_{max})$, and execution time is asymptotically optimal with $\mathcal{O}(\ell_{max}^{2}\log^{2}\ell_{max})$ operations. This FSHT is implemented in Julia with the package FastTransforms [37] (as are the NUFFT methods from [33] used in Step 2).

The FSHT in FastTransforms assumes the input function has a bivariate Fourier expansion of the form

[TABLE]

Any function on the sphere is required to have these even/odd conditions on its bivariate Fourier coefficients [23]. At the end of step 2 we have obtained the bivariate Fourier expansion of the data of the form

[TABLE]

where $p=N_{side}/2$. Since we are dealing with real-valued data, we can expand Fourier coefficients array in $\lambda$ to an odd number of points. The expanded array is defined by

[TABLE]

Using the array $X$, we can write (19) as

[TABLE]

from which we can obtain the coefficients $g_{j}^{k}$ in (18).

The FSHT software takes bivariate Fourier coefficients $g_{j}^{k}$ as input in an array organized as follows:

[TABLE]

The output of the software is the approximate spherical harmonic coefficients of the data arranged in an array of the form

[TABLE]

The angular power spectrum (1) can then be computed from this array.

4 Numerical Results

In this section we present a few numerical tests comparing the spherical harmonics and angular power spectrum (1) computed by our new method HP2SPH to the values computed by the HEALPix software employing the iterative scheme (6), pixel weights, and ring weights (7). The first test compares the rate at which the two methods converge to the spherical harmonic coefficients by applying them to deterministic (i.e. non-noisy) functions sampled at the HEALPix points with known coefficients. The second test compares the accuracy of the methods using deterministic functions, which have a known power spectrum. In the third test, we compare the methods after calculating the angular power spectrum for a real CMB data map, which contains noise.

4.1 Convergence of Spherical Harmonic Coefficients

We choose the test function

[TABLE]

where $\mathbf{x}(\lambda,\theta)=[x(\lambda,\theta)\;y(\lambda,\theta)\;z(\lambda,\theta)]$ from (11) and the parameters, which we picked randomly, are given by

[TABLE]

The function $(2-2\mathbf{x}(\lambda,\theta)\cdot\mathbf{x}(\lambda_{\rm c},\theta_{\rm c}))^{3/2}$ is a called a potential spline of order $3/2$ centered at $\mathbf{x}(\lambda_{\rm c},\theta_{\rm c})$ and its exact spherical harmonic coefficients are given by [18]

[TABLE]

These values are used to compare the convergence rates of the methods to the exact spherical harmonic coefficients of $f$. We do this by plotting in Figure 6 the maximum absolute errors of the coefficients against the parameter $t$, which is used to determine the grid resolution parameter ($N_{side}=2^{t}$). Note that the spherical harmonic coefficients of the CMB decay at a rate between $\mathcal{O}(\ell^{-2})$ and $\mathcal{O}(\ell^{-3})$ [39], which is slower than the decay rate of the coefficients of the test function (20) (since, for all $-\ell\leq m\leq\ell$, $|Y_{\ell}^{m}(\lambda_{\rm c},\theta_{\rm c})|\leq\sqrt{(2\ell+1)/4\pi}$, and the remaining terms in (21) decay at a rate of $\mathcal{O}(\ell^{-5})$). This means that the test function has more smoothness than an actual CMB data set.

Figure 6(a) compares the four methods and shows that the HP2SPH method converges to the spherical harmonic coefficients of (20) at a rate at least twice as fast as any of the HEALPix methods. Although the consecutive iterative refinement steps used in the HEALPix method produce progressively better errors, Figure 6(b) illustrates that this does not improve the convergence rate (as discussed in Section 2.2). It is also important to note that after 8 iterative steps, there are no further improvements in the accuracy, indicating the algorithm has nearly converged to the least squares solution to (3). The results the HEALPix method with pixel weights look pretty erratic with convergence achieved to 8-10 digits around $t=7$, but no further reductions. This could be because of potential errors in the computed quadrature weights.

Next, we test how the convergence rates of the methods are affected by high frequencies. To do this we add 15 spherical harmonics to (20) with the following randomly chosen degrees and orders:

The results from this test are displayed in Figure 7. We note that in order to ensure the calculation of asymptotic convergence rates for all methods, we excluded the errors for $t=8$ in the lines of best fit. This test shows that the convergence rates of the new method as well as the iterative HEALPix scheme and the HEALPix method with ring weights are not affected by the addition of high degree spherical harmonic terms to (20). The HEALPix method with pixel weights shows a similar erratic behavior to Figure 6(a), with no steady reductions in the errors after $t=9$. The new HP2SPH method has the lowest errors of all the methods.

4.2 Errors in the Angular Power Spectrum

In this test, we first explore the accuracy of all the methods for computing the angular power spectrum of (20). These results are compared to the exact spectrum, which is calculated using the exact spherical harmonic coefficients (21) in (1).

The angular power spectrum of the four methods are displayed in Figure 8(a). We see that the algorithms produce similar results for lower degrees $\ell$, but the HEALPix method with the iterative refinement scheme (6) diverges for degrees greater than approximately $\ell=100$, whereas the ring weight and pixel weight quadrature (7) results diverge for degrees greater than $\ell=200$. To compare the methods more directly, Figure 8(b) plots the absolute errors in the angular power spectrum for each degree $\ell$. While the HEALPix method with ring weights performs the best out of all of the HEALPix methods for $\ell\leq 50$, the errors increase rapidly for higher $\ell$. Conversely, the HEALPix method with pixel weights does not perform as well for smaller $\ell$, yet it performs better than the other HEALPix methods for larger $\ell$. The HP2SPH method offers comparable results for low $\ell$ as the pixel weights method while still maintaining accuracy for high $\ell$.

Similar to Test 1, we test the accuracy of the methods when computing the power spectrum of data with high frequencies, as occur in real CMB data maps. As before, we add several spherical harmonic functions of high degree to the function (20). The new test function appends the following degrees and orders:

The power spectrum of this function is the same as (20), but with the value at each degree $\ell$ of appended spherical harmonics increased by $\frac{1}{2\ell+1}$.

Figure 9(a) displays the errors in the angular power spectrum of this new function for all the methods over all $\ell$, while Figure 9(b) displays the errors only over the range of $\ell$ that were appended to the base function. We again see that the ring weights provide the highest accuracy for low $\ell$, but then the errors increase rapidly, while the HEALPix iterative method has the largest errors for low $\ell$, but some of the smallest errors for the $\ell$ corresponding to the appended spherical harmonics. In contrast, the new HP2SPH method provides small errors over the entire angular power spectrum, clearly showing its benefits over the HEALPix methods for determining the angular power spectrum of deterministic functions on the sphere.

4.3 Application to Real CMB Map

Our final numerical test compares the methods on the real CMB map shown in Figure 1. Figure 10 (a) shows the angular power spectrum for this map computed with the methods, while (b) shows the errors in the three HEALPix methods compared to the new HP2SPH method. We see from the figure that the new method (in blue) produces visually the same results as the HEALPix methods (red, magenta, and cyan), indicating that the new method is not anymore susceptible to noise than the HEALPix methods.

5 Conclusions and Remarks

We have presented a new method, HP2SPH, for performing discrete spherical harmonic analysis on data collected using the HEALPix scheme. The method utilizes the FFT, NUFFT, and the FSHT to compute the spherical harmonic transform in near optimal computational complexity ($\mathcal{O}(N\log^{2}N)$ complexity for $N$ total HEALPix points). Several numerical tests were presented to demonstrate the improved convergence and accuracy of the new method relative to the various HEALPix approaches for problems involving synthetic data with no noise, except that introduced by roundoff errors. In the case of a real CMB map with additional types of noise, the power spectra computed by the methods show good agreement. The new HP2SPH benefits further from the backward stability of the FSHT, which ensures that the resulting uncertainty in the spherical harmonic coefficients has only a low algebraic growth with respect to degree and is always proportional to the norm of the noise in the data. We anticipate the new method will be applicable to the many other areas where the HEALPIx scheme is used and is naturally generalizable to other equal-area isolatitudinal sampling schemes for the sphere.

For our next steps, we will work to optimize the implementation of the method, which is currently in Julia, to improve its actual run-time. This will include transcribing our code into a lower-level language like C; efforts in this direction are already underway for the FSHT [36]. In addition to this, we will include the ability to perform Fourier synthesis on a CMB map, i.e. given an angular power spectrum, we will return the corresponding CMB map values. For this purpose, our method has another advantage over HEALPix in that we will have the bivariate Fourier coefficients, which will simply make the synthesis an application of the FFT and NUFFT. Finally, we plan to add functionality for analyzing the polarization of CMB maps.

Acknowledgements

We are grateful to Ann Almgren for her comments and suggestions as well as Richard Mikael Slevinsky for assisting us with using the FastTransforms package. We also thank the entire HEALPix team, especially Kris Górski and Martin Reinecke, for their feedback and assistance on this work. Finally, we thank the anonymous reviewers whose comments helped improve the paper. The first author’s work was supported by the SMART Scholarship, which is funded by The Under Secretary of Defense-Research and Engineering, National Defense Education Program / BA-1, Basic Research. The second author’s work was supported by National Science Foundation grant CCF 1717556.

Appendix A Spherical Harmonic Conventions

We denote a scalar spherical harmonic of degree $\ell\geq 0$ and order $-\ell\leq m\leq\ell$ as $Y_{\ell}^{m}(\lambda,\theta)$, where $\lambda$ is the azimuth angle and $\theta$ is the zenith angle. We define these functions as

[TABLE]

where $Y_{\ell}^{m}=(-1)^{m}Y_{\ell}^{|m|}$ for $m<0$ and $P_{\ell}^{m}(\cos\theta)$ are the associate Legendre functions. As eigenfunctions of the Laplace-Beltrami operator, spherical harmonics are the natural basis for square integrable functions on the sphere [1]. In other words, any $L^{2}$-integrable function $f$ on the sphere can be uniquely represented as

[TABLE]

where the spherical harmonic coefficients, $\widetilde{a}_{\ell}^{m}$, are found using the usual $L^{2}$-inner product for scalar functions on the sphere:

[TABLE]