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 O(Nlog2N), with an initial set-up cost of O(N3/2logN). This compares
favorably with O(N3/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