Numerical analysis of functions on the sphere involves (1) a class of mathematical operations, whose objects are (2) discretised maps, i.e. quantizations of arbitrary functions according to a chosen tessellation (exhaustive partition of the sphere into finite area elements). Hereafter we mostly specialise our discussion to CMB related applications of HEALPix, but all our statements hold true generally for any relevant deterministic and random functions on the sphere.
Considering point (1): Standard operations of numerical analysis which one might wish to execute on the sphere include convolutions with local and global kernels, Fourier analysis with spherical harmonics and power spectrum estimation, wavelet decomposition, nearestneighbour searches, topological analysis, including searches for extrema or zerocrossings, computing Minkowski functionals, extraction of patches and finite differencing for solving partial differential equations. Some of these operations become prohibitively slow if the sampling of functions on the sphere, and the related structure of the discrete data set, are not designed carefully.
Regarding point (2): Typically, a whole sky map rendered by a CMB experiment contains (i) signals coming from the sky, which are by design strongly bandwidth limited (in the sense of spatial Fourier decomposition) by the instrument's angular response function, and (ii) a projection into the elements of a discrete map, or pixels, of the observing instrument's noise; this pixel noise should be random, and white, at least near the discretisation scale, with a bandwidth significantly exceeding that of all the signals.
With these considerations in mind we propose the following list of desiderata for the mathematical structure of discretised full sky maps:
1. Hierarchical structure of the data base. This is recognised as essential for very large data bases, and was postulated in construction of the Quadrilateralized Spherical Cube (or quadsphere, see http://space.gsfc.nasa.gov/astro/cobe/skymap_info.html), which was used for the COBE data. An argument in favour of this proposition states that the data elements which are nearby in a multidimensional configuration space (here, on the surface of a sphere), are also nearby in the tree structure of the data base, hence the nearneighbour searches are conducted optimally in the data storage medium or computer RAM. This property, especially when implemented with a small number of base resolution elements, facilitates various topological methods of analysis, and allows easy construction of wavelet transforms on quadrilateral (and also triangular) grids. Figure 1 shows how a hierarchical partition with quadrilateral structure naturally allows for a binary vector indexing of the data base.

2. Equal areas of discrete elements of partition. This is advantageous because (i) white noise generated by the signal receiver gets integrated exactly into white noise in the pixel space, and (ii) sky signals are sampled without regional dependence, except for the dependence on pixel shapes, which is unavoidable with tessellations of the sphere. Hence, as much as possible given the experimental details, the pixel size should be made sufficiently small compared to the instrument's resolution to avoid any excessive, and pixel shape dependent, signal smoothing.
3. IsoLatitude distribution of discrete area elements on a sphere. This property is critical for computational speed of all operations involving evaluation of spherical harmonics. Since the associated Legendre polynomial components of spherical harmonics are evaluated via slow recursions, and can not be simply handled in an analogous way to the trigonometric Fast Fourier Transform, any deviations in the sampling grid from an isolatitude distribution result in a prohibitive loss of computational performance with the growing number of sampling points, or increasing map resolution. It is precisely this property that the COBE quadsphere is lacking, and this renders it impractical for applications to high resolution data.
A number of tessellations have been used for discretisation and analysis of functions on the sphere (for example, see Driscoll & Healy (1994), Muciaccia, Natoli & Vittorio (1998), Doroshkevich et al. (2005)  rectangular grids, Baumgardner & Frederickson (1985), Tegmark (1996)  icosahedral grids, Saff & Kuijlaars (1997), Crittenden & Turok (1998)  `igloo' grids, and Szalay & Brunner (1998)  a triangular grid), but none satisfies simultaneously all three stated requirements.
All three requirements formulated above are satisfied by construction with the Hierarchical Equal Area, isoLatitude Pixelation (HEALPix) of the sphere, which is shown in Figure 2. A more detailed description of HEALPix, its motivations, and applications can be found in Górski et al (2005).

Version 3.31, 20170106