## Abstract

We construct near-optimal quadratures for the sphere that are invariant under the icosahedral rotation group. These quadratures integrate all (*N*+1)^{2} linearly independent functions in a rotationally invariant subspace of maximal order and degree *N*. The nodes of these quadratures are nearly uniformly distributed, and the number of nodes is only marginally more than the optimal (*N*+1)^{2}/3 nodes. Using these quadratures, we discretize the reproducing kernel on a rotationally invariant subspace to construct an analogue of Lagrange interpolation on the sphere. This representation uses function values at the quadrature nodes. In addition, the representation yields an expansion that uses a single function centred and mostly concentrated at nodes of the quadrature, thus providing a much better localization than spherical harmonic expansions. We show that this representation may be localized even further. We also describe two algorithms of complexity for using these grids and representations. Finally, we note that our approach is also applicable to other discrete rotation groups.

## 1. Introduction

Many problems in physics, mathematics and engineering involve integration and interpolation on the sphere in . Of particular importance are discretizations of rotationally invariant subspaces of that integrate all spherical harmonics up to a fixed order and degree. A typical approach to discretizing the sphere is that of equally spaced discretization in azimuthal angle and Gauss–Legendre discretization in polar angle, leading to an unreasonably dense concentration of nodes near the poles. It is well known that, in a variety of applications, such concentration of nodes may lead to problems when using these grids.

Alternatively, Sobolev (1962) (see also McLaren 1963; Sobolev 1992) suggested the use of grids that are invariant under finite rotation groups. In such constructions, there is no clustering of nodes and, moreover, the number of nodes necessary to integrate a particular subspace is close to optimal. In fact, in some cases it is optimal, a notion we make precise later. There have been several constructions of this type of grids using algebraic approaches (see McLaren 1963; Lebedev 1976; Lebedev & Laikov 1999; Popov 2008; and references therein), but the technique limits the results to a few particular orders and degrees. An attempt to construct grids invariant under the icosahedral group may be found in Konyaev (1977, 1982) (with some negative weights in the early construction). We also refer to and Stroud (1971) for a review and further references.

In this paper, we develop a systematic numerical approach for constructing nearly optimal quadratures invariant under the icosahedral group to integrate rotationally invariant subspaces of up to a fixed order and degree. Using these grids and a reproducing kernel, we show how to replace the standard basis of spherical harmonics on a rotationally invariant subspace by a representation formed using a single function centred at the quadrature nodes. The reproducing kernel is mostly concentrated near the corresponding grid point. In the resulting representation, the coefficients, up to a factor, are the values on the grid of the function being represented. We may interpret this construction as an analogue of Lagrange interpolation on the sphere (see also Sloan 1995; Sloan & Womersley 2000) and note that it allows us to develop well-conditioned linear systems for interpolation in contrast to some earlier constructions (Fernandez 2005; Fernandez & Prestin 2005).

An alternative to spherical harmonics has long been sought, especially for numerical purposes. Spherical harmonics provide an efficient orthonormal basis, nicely subdivided into rotationally invariant subspaces. However, the global support of these functions poses a serious difficulty in problems where physical effects are localized. In fact, the global nature of spherical harmonics is a consequence of their optimality. Therefore, if we want localized functions to represent the same subspaces, we necessarily must have a less efficient representation.

Work on alternative representations has accelerated since the development of wavelets. Several proposals for local and multi-resolution representations have been suggested over the years, resulting in numerous publications (e.g. Taylor *et al.* 1997; Freeden *et al.* 1998; Beylkin & Cramer 2002; Guilloux *et al.* 2009; and references therein). Several attempts to generalize the wavelet transform on Euclidean spaces to the sphere have been made (e.g. Holschneider 1996; Freeden *et al.* 1998; Antoine *et al.* 2002; bib20; and references therein). One particular generalization (Antoine *et al.* 2002) uses a doubly periodic discretization of the sphere to provide algorithmic efficiency. Such discretization, however, does not assure rotational invariance, as the nodes will necessarily concentrate near the polar regions. Discrete multi-resolution transforms on the sphere have been suggested using geodesic grids (Schröder & Sweldens 1995), equal-area subdivisions (Lessig 2007) and irregular grids (Daubechies *et al.* 1999; Le Gia & Mhaskar 2008), mostly for the purpose of graphics. These transforms are also useful for analysis of data measured on the sphere, e.g. for geophysical or astrophysical applications. Such transforms may have a formal extension to high-order versions, but only low-order constructions appear to be practical.

A practical high-order local construction based on a ‘cubed sphere’ has been successfully used in global atmospheric modelling (Taylor *et al.* 1997) and geodesy (Beylkin & Cramer 2002). As a replacement for spherical harmonics, such ‘cubed-sphere’ representations proved to be practical in the development of high-order numerical techniques but, at the same time, they lack a number of useful properties.

We view our approach as the first step in constructing a local and multi-resolution representation of functions on the sphere that respects rotationally invariant subspaces. We note that the high efficiency of quadratures constructed in this paper implies a near-uniform distribution of nodes on the sphere. On the other hand, the nodes maintain a regular organization visually similar to that of geodesic or equal-area grids. Moreover, our grids are associated with rotationally invariant subspaces, an important property in a number of numerical applications, e.g. geodesy. To date, we have constructed grids that integrate subspaces of maximum order and degree ranging from 5 up to 210. As an example, we illustrate in §3*b* a grid with 7212 nodes integrating subspaces of maximal order and degree 145.

Our method for constructing these grids is based on Newton’s method, and hence the key to its success is a good initial distribution of nodes. At first, as the initial distribution, we used spherical codes with icosahedral symmetry constructed by Hardin *et al.* (2000). We then developed a simple approach to map a two-dimensional lattice to the icosahedron and then to the sphere, yielding the initial distribution of nodes structurally similar to that of the resulting quadrature nodes. We describe this approach in §3*c*.

The rotationally invariant spherical grids constructed here have many applications. Let us mention a few specific problems in some detail. First, due to the central role played by spherical harmonics in the theory of gravity and magnetic fields, solutions to many geodetic problems use them as a basis (e.g. Heiskanen & Moritz 1967; Backus *et al.* 1996). Yet, their global support is inconsistent with the physical nature of the problem, leading to many difficulties in, for example, constructing the gravity models. The grids developed in this paper provide a first step towards replacing spherical harmonics with localized functions. We plan to continue work in this direction. Second, the equations used in global atmospheric modelling are typically posed on the sphere (Williamson *et al.* 1992). Current spectral methods that use spherical harmonics suffer from the above-mentioned problems of nodal clustering and require additional steps to alleviate the problem (Jakob-Chien & Alpert 1997). An approach in Taylor *et al.* (1997) not only provides a practical alternative, but also introduces artificial singularities associated with vertices of the cube. The new representations developed in this paper eliminate clustering and singularities due to the coordinate system and should provide efficient solution methods. Third, acoustic and electromagnetic scattering problems posed as integral equations involve integration over spherical domains Colton & Kress (1992). New algorithms for the numerical solution of these integral equations may be based on the results of this paper. Finally, we mention a numerical technique used in molecular dynamics calculations known as discrete variable representation (DVR; Light *et al.* 1985). In the DVR method, integrals over spherical domains in are computed, and some implementations of the DVR method (Haxton 2007) use quadratures developed by Lebedev (1976) and Lebedev & Laikov (1999). Our quadratures should extend such methods by effectively allowing an arbitrary order and degree.

We start by reviewing necessary mathematics, outline a general method for constructing nodes invariant under the icosahedral group and illustrate resulting quadratures through several examples. Based on these new quadratures, we also develop a representation of functions in rotationally invariant subspaces as linear combinations of a single function centred at the quadrature nodes. We then discuss properties of this representation that may be interpreted as an analogue on the sphere of the usual Lagrange interpolation. We present two algorithms to work with these new grids (here, *N* is the maximum order and degree of the subspace). We also consider a representation via localized functions and briefly describe its properties. We conclude by listing directions for future research in developing practical methods for applications dealing with the sphere.

## 2. Preliminaries

Here, we establish notation and state some well-known results about spherical harmonics, reproducing kernels and quaternionic representation of rotations. We also formulate two theorems of Sobolev (1962) for quadratures that are invariant under a finite rotation group and a theorem of Molien (1897) on the number of invariants for a finite group.

### (a) Spherical harmonics and reproducing kernels

We denote the sphere in as An orthonormal basis for is given by the spherical harmonics
2.1
where the polar angle *θ*∈[0,*π*], the azimuthal angle *ϕ*∈[0,2*π*) and are the normalized associated Legendre functions given by
Any may be expanded as
2.2
with the coefficients given by
2.3
where 〈·,·〉 is the inner product and * denotes complex conjugation. We define a subspace of spherical harmonics with fixed degree *n* as
2.4
The spherical harmonics are linearly independent and hence the dimension of is 2*n*+1. The subspace of maximum degree *N* is then the direct sum
2.5
and has dimension (*N*+1)^{2}. We make use of the addition theorem (e.g. Gradshteyn & Ryzhik 1994), which states that, for ,
2.6
where *P*_{n} is the Legendre polynomial of degree *n*. We also use the reproducing kernel for that yields
2.7
where
2.8
The identity in equation (2.7) may be verified by using equation (2.6). We rely on equation (2.7) to develop a new representation of functions in and an analogue of Lagrange interpolation on the sphere.

### (b) Quaternionic representation of rotations and points on the sphere

We use quaternions to represent points on the sphere and to perform rotations. A quaternion **q** is defined as
where and the symbols **i**,**j** and **k** satisfy **i**^{2}=**j**^{2}=**k**^{2}=**ijk**≡−1 and **ij**=**k**, **j k**=

**i**and

**=**

**k**i**j**. Here,

*w*is the scalar part of the quaternion and (

*x*,

*y*,

*z*) is its vector part, which may be associated with a vector in . Thus, vectors in may be represented by quaternions by simply setting the scalar part to zero. The norm of

**q**is defined as Quaternions form an algebra, , and the multiplicative inverse of

**q**is given by We represent points on the sphere as quaternions with zero scalar part and unit norm, and ∥

**q**∥=1.

To rotate a point on the sphere represented by by an angle *θ* around the unit vector , we construct the quaternion, ,
The rotated point is given by the quaternion , and it is not difficult to show that **v**′ has zero scalar part and unit norm (e.g. Kuipers 1999). We use quaternions to represent actions of the icosahedral group (see electronic supplementary material).

### (c) Sobolev’s theorems

Sobolev’s (1962) paper contains two key results that we summarize now. For clarity, we specialize these results to the case of the icosahedral rotation group, noting that our approach may be applied without essential modifications to other discrete rotation groups.

Consider evaluating the integral , where *f* belongs to the rotationally invariant subspace , using a quadrature rule with nodes and weights so that
for all . As both the integration domain and the subspace have rotational symmetry, it is natural to look for rotationally invariant quadrature rules. Following Sobolev (1962), we define invariant quadrature rules as follows: given a finite rotation group *G*, a quadrature rule *Q* is invariant under *G* if, for all *g*∈*G*,
We then have the following theorems.

### Theorem 2.1

*Let Q be a quadrature rule invariant under the group G. Then, Q is exact for all functions* *if and only if Q is exact for functions f invariant under G*.

This theorem effectively reduces the size of the system of nonlinear equations that must be solved to determine a quadrature invariant under the group *G*. The next result gives a formula to calculate the number of invariant functions under the group *G* in a subspace of spherical harmonics of a given degree *n*. Let *q*_{1}=5 be the number of edges meeting at a vertex of an icosahedron, *q*_{2}=3 be the number of sides of its (triangular) face and *q*_{3}=2 denote the order of rotation about mid-points of opposing edges.

### Theorem 2.2

*For a given degree n, the number of functions invariant under the icosahedral rotation group in a subspace of spherical harmonics* is given by
*where* ⌊ ⌋ *denotes the integer part*.1

This result may also be obtained using a theorem of Molien (1897). In that approach, the number of invariants for a finite group may be obtained as coefficients of a generating function. In our case, we (see Meyer 1954) obtain the following theorem.

### Theorem 2.3

*For a given degree n, the number of functions invariant under the icosahedral rotation group in a subspace of spherical harmonics* *is given by coefficients S*(*n*) *of the series expansion of the generating function*

It is not difficult to see that both theorems 2.2 and 2.3 yield the same result. We use theorems of this section to determine the number of equations contributed by each subspace to the nonlinear system of equations determining the quadrature nodes and weights.

## 3. Quadratures for the sphere

The main difficulty in constructing quadratures comes from the need to solve a large system of nonlinear equations. Without using a special structure of these equations, general root-finding or optimization methods typically fail. The essence of our approach is to develop and use such a structure within a root-finding method.

We start by noting four different types of orbits of the icosahedral rotation group. In general, with three exceptions described below, a point on the sphere under the action of the group generates a total of 60 points (which is the group’s order). However, if a point is a vertex of the icosahedron, then it generates a total of only 12 distinct points. Also, if a point is the projection of the centre of an icosahedron face onto the sphere, it generates 20 distinct points in total. Finally, if a point is the projection onto the sphere of the midpoint of an edge, it generates a total of 30 distinct points. When describing these different types of orbits, it is sufficient to consider a single point, a generator of its orbit. The orbit of a point with spherical coordinates (*θ*,*ϕ*) is the set and, depending on the type of orbit, has cardinality 12, 20, 30 or 60.

With these types of orbits in mind, we consider four types of quadratures. The first type assumes that all generators, except for a vertex of the icosahedron, give rise to orbits of size 60, i.e.
3.1
where are the coordinates of the vertices of an icosahedron inscribed in the unit sphere, *w*_{v} are their associated weights and *N*_{g} is the number of generators with coordinates and weights . For each *g*_{i}∈*G*, we denote . For this quadrature rule, we fix the positions of vertices and hence do not consider them as unknowns. The subscript v indicates that the vertices are held fixed. Therefore, for this type of quadrature, there are 3*N*_{g}+1 unknown generator coordinates and weights, although there are a total of 60*N*_{g}+12 nodes.

The second type of quadrature has nodes at the fixed vertices of an icosahedron and one generator with a fixed position at the projection of the centre of an icosahedron face onto the sphere, giving rise to 20 points on the sphere. We assume that all other generators give rise to orbits of length 60. Thus, we have
3.2
where are the coordinates of the face centres of the icosahedron projected onto the sphere and *w*_{f} is the associated weight. For this quadrature, there are 3*N*_{g}+2 unknown generator coordinates and weights, with a total of 60*N*_{g}+32 nodes.

The third type of quadrature has nodes at the fixed vertices of an icosahedron and one generator with a fixed position at the projection onto the sphere of the midpoint of an edge, giving rise to 30 fixed nodes of the quadrature. We have
3.3
where are the coordinates of the projection onto the sphere of the midpoints of edges of the icosahedron and *w*_{e} is the associated weight. For this quadrature, there are 3*N*_{g}+2 unknown generator coordinates and weights, with a total of 60*N*_{g}+42 nodes.

Finally, the fourth type of quadrature has, as fixed nodes, vertices and projections of both face centres and edge midpoints of the icosahedron. We have
3.4
For this quadrature, we have 3*N*_{g}+3 unknown generator coordinates and weights, with a total of 60*N*_{g}+62 nodes.

Using theorems 2.2 or 2.3, we determine the number of invariant functions in the subspace . Theorem 2.1 allows us to limit the number of equations to not exceed the number of invariant functions. Ideally, when constructing a system of equations for the generator coordinates and weights, we would like to match the number of equations to the number of unknowns. For example, to construct a quadrature that integrates exactly the subspace , we must integrate the 10 invariant functions in . Therefore, we have 10 equations and, using quadrature (3.1), we have 3*N*_{g}+1 unknowns, where *N*_{g} is the number of generators. Setting 3*N*_{g}+1=10 gives *N*_{g}=3 and thus we look for a quadrature with 192=3×60+12 nodes. We solve the corresponding system (see §3*b*) and obtain a solution, thus verifying its existence directly. We illustrate a set of 7212 nodes that exactly integrates in §3*b*. This set of nodes was found using the quadrature *Q*_{v} in equation (3.1).

However, it appears that, in some cases, so constructed system of equations may not have a solution. Our conclusion is based on the behaviour of Newton’s iteration and so far has not been verified analytically. In such cases, we reduce the number of equations by removing those from the subspace of the highest degree and solve a formally under-determined system. We note that, in this situation, Newton’s iteration may not converge quadratically, and the resulting quadratures may be less efficient than the count of invariants suggests. However, the reduction in efficiency is negligible for practical purposes. We also note that there may be more than one solution of these equations with different positive weights. In all cases, we directly verify the resulting quadratures for all spherical harmonics.

It is also possible to set two or more of the generator weights equal, thereby reducing the number of unknowns. We have confirmed numerically the existence of quadratures where some of the generator weights are identical. We note that there has been interest in constructing quadratures for the sphere where all weights are equal, i.e. so-called spherical t-designs (e.g. Conway & Sloane 1999 and references therein). The approach outlined above allows us to address this problem as well. We note, however, that it might be more natural to pose a problem of finding quadratures with two or more distinct weights depending on the type of quadrature in equations (3.1)–(3.4), where we set *w*_{j}=*w*_{1} for all .

## Remark 3.1

For quadratures *Q*_{ve} and *Q*_{vfe}, the number of unknowns may be slightly different if the node is on an edge, but not at its centre. In this case, the orbit has length 60, as in the general case, but there are only two degrees of freedom because the node is constrained to lie on the great circle joining two vertices of the icosahedron.

We also note that additional quadratures may be constructed by using the full icosahedral group that includes reflections. In this paper, we focus only on quadratures (3.1)–(3.4), and will report results of developing other quadratures in a future paper.

### (a) Efficiency of quadratures on the sphere

A useful way of measuring the efficiency of a quadrature was suggested in McLaren (1963) using the ratio of the dimension of the subspace to be integrated to the (maximum) number of degrees of freedom in the quadrature rule (two coordinates and a weight for each node),
3.5
where *M* is the number of nodes. If *η*=1, we call the quadrature optimal. Remarkably, there are a few special cases where *η* is slightly greater than 1 (see McLaren 1963). However, in general, the efficiency is less than 1. As the number of nodes *M* increases, the efficiency of each of the four types of quadratures, *Q*_{v}, *Q*_{vf}, *Q*_{ve} and *Q*_{vfe}, improves and approaches 1. In figure 1, we display the potential efficiency of quadrature *Q*_{v} as a function of the degree *N* of subspace using theorems 2.2 or 2.3 to count invariants and superimpose the actual efficiencies of computed quadratures. The behaviour of efficiency of other quadratures in equations (3.2)–(3.4) is similar.

For comparison, the efficiency of the standard discretization of the sphere, using equally spaced nodes in the azimuthal direction and Gauss–Legendre nodes in the polar direction, is *η*_{st}=2/3. Specifically, for a maximum order and degree *N*, we need at least *N*+1 equally spaced nodes to integrate in the azimuthal direction over the interval [0,2*π*) and (*N*+1)/2 nodes to integrate in the polar direction with a Gauss–Legendre quadrature over the interval [0,*π*], so that *M*=(*N*+1)^{2}/2. The efficiency of the standard quadrature is also illustrated in figure 1.

Although it is beneficial to have a highly efficient quadrature, it is the near-uniform distribution of nodes that is the key property of quadratures with icosahedral symmetry. As mentioned in §1, clustering of nodes is detrimental for many applications. For example, the smallest distance between the nodes determines the maximum time-step size when solving partial differential equations on the sphere. Hence, when using the standard discretization of the sphere, dynamics near the poles forces unnecessarily small time steps to be taken. Although there are procedures to deal with this problem (e.g. Jacob-Chien & Alpert 1997), with our new quadratures, the maximum time-step size is determined by the dynamics rather than the grid.

### (b) Iterative method for computing quadratures via Newton’s iteration

We are now in a position to describe our algorithm for finding nodes and weights for quadratures invariant under the icosahedral group. As the algorithm is the same for all quadratures (3.1)–(3.4), we use the generic notation *Q* to denote any one of the four types. For a given *N*, we first use theorems 2.2 or 2.3 to determine those *n*, 0≤*n*≤*N*, for which *S*(*n*)≠0. This gives a set of values of *n* that we denote as . For example, we obtain . We now define the set of orders *m* associated with a fixed degree *n* as
3.6
Note that, for odd *m*, vanishes identically by symmetry. Using indices with , we now construct the equations
3.7
where and *δ* is the Kronecker delta. The collection of equations in (3.7) is a system of nonlinear equations for the generator coordinates and weights and the weights associated with the fixed nodes. The group actions in (see equations (3.1)–(3.4)) are carried out using quaternionic multiplication, as described in §2*b*. Our formulation implicitly assumes that the invariant functions generated numerically from equation (3.6) for are independent. It has been our observation that this is the case because any missing invariant functions would have been discovered during *a posteriori* verification of quadratures.

Recall that Newton’s method for solving the system of nonlinear equations **F(x)**=**0** with , where *n*′≤*n*, is the iteration
3.8
starting with some initial vector **x**_{0}, where **J** is the Jacobian of **F**,
3.9
and is its generalized inverse.

Writing equation (3.7) in vector notation **F= 0**, we use equation (3.8) to solve for the generator coordinates and weights and the weight associated with the fixed nodes. The Jacobian is calculated analytically, as described in appendix A. As a way of finding initial weights, we use the initial node distribution to form and solve the linear system
where the matrix has entries , vector

**I**=(1,1,…,1)

^{t}, and

**w**is the vector of initial weights. Construction of the initial node distribution is described in §3

*c*.

In figure 2, we illustrate a quadrature with 7212 nodes that integrates functions in . The light coloured nodes correspond to the vertices of the icosahedron, whereas the darker ones are those constructed by group action on the generators. Note that each of the dark coloured nodes has six nearest neighbours, whereas the vertices of the icosahedron have only five nearest neighbours. So far, the largest subspace we constructed quadratures for is .

### (c) Generation of initial grid

At first, as the initial node distribution, we used spherical codes with icosahedral symmetry generated by Hardin *et al.* (2000). In using such initial distributions, we observed that the key to convergence of Newton’s method has been the structural similarity of the initial grid with the final quadrature grid. This led us to a simple construction that assures such structural similarity.

Let us start by generating a lattice in the plane formed by equilateral triangles
with the lattice vectors **e**_{1}=(1,0) and . We then select three points,
on the lattice to form an equilateral triangle with vertices *P*_{0}, *P*_{1} and *P*_{2}, where are positive integers. Indeed, it is easy to check that |*P*_{0}*P*_{1}|^{2}=|*P*_{0}*P*_{2}|^{2}=|*P*_{1}*P*_{2}|^{2}=*s*^{2}+*t**s*+*t*^{2}. We consider *t* and *s* as parameters that uniquely determine the number and type of generators for which the triangle *P*_{0}*P*_{1}*P*_{2} serves as a template. Thus, *t* and *s* uniquely determine the total number of unknown parameters in the nonlinear system (3.7). We illustrate this construction in figure 3.

The centre of the triangle *P*_{0}*P*_{1}*P*_{2} is at the point
which may or may not be a point of the lattice. If the centre coincides with a lattice point, then the type of quadrature is one with a face-centred node. In the same manner, for some choices of parameters *t* and *s*, there might be points on the sides of the triangle and, in such case, it leads to a quadrature with edge-centred nodes. As we pointed out already, the type of quadrature affects the count of independent variables.

To find the generators, we select nodes in the interior of a fundamental region of the icosahedral rotation group, defined by connecting the centre of the triangle with any two of its vertices (Grove & Benson 1985). We use the triangle *P*_{0}*P*_{1}*P*_{c} (see figure 3 for an illustration of this particular fundamental region). We note that, by definition, acting with the icosahedral rotation group on points inside the fundamental region will produce points on all faces of the icosahedron and thus it is sufficient to consider only points inside that region and, in addition, points on the sides of the triangle *P*_{0}*P*_{1}*P*_{c}.

In placing points from the fundamental region to the surface of the icosahedron inscribed into the unit sphere, it is convenient to use the barycentric system of coordinates associated with the triangle *P*_{0}*P*_{1}*P*_{c}. The barycentric coordinates (*τ*_{0},*τ*_{1},*τ*_{c}) of a lattice point {*m***e**_{1}+*n***e**_{2}} are found by explicitly solving the 3×3 linear system
which yields

We recall that a point is inside the triangle *P*_{0}*P*_{1}*P*_{c} if *τ*_{0},*τ*_{1},*τ*_{c}>0, and it is on a side if one of the barycentric coordinates is equal to 0. To position the generators onto the sphere, we first map the points from the fundamental region formed by the triangle *P*_{0}*P*_{1}*P*_{c} to a face of the icosahedron. Let vectors define a face of the icosahedron and the centre of the face. A point with the barycentric coordinates (*τ*_{0},*τ*_{1},*τ*_{c}) with respect to the triangle *P*_{0}*P*_{1}*P*_{c} is then mapped to the point . We then project this point radially outward to the surface of the unit sphere. Using the resulting generators, we apply group actions to produce all the nodes on the surface of the sphere.

Finally, we note that, as the number of generators grows, we use previously constructed quadrature nodes to map a grid template on the face of the icosahedron to the sphere. We use the barycentric coordinates of points within the corresponding triangles on the face of the icosahedron formed by pre-images of previous quadrature nodes on the sphere.

## 4. An analogue of Lagrange interpolation and an alternative representation for rotationally invariant subspaces of

Using the results in §§2 and 3, we now construct an alternative representation for functions for invariant subspaces of . Let us discretize equations (2.7) and (2.8) using quadratures developed in §3. The quadrature must be chosen so that it integrates exactly subspaces of order and degree at least 2*N* because this is the maximal order and degree of the product *K*(** ω**·

**′)**

*ω**f*(

**′) in equation (2.7). Using such a quadrature, we have the identity 4.1 where are unit vectors denoting the coordinates of the quadrature nodes and**

*ω**w*

_{j}are the associated weights.

We observe that if , then equation (4.1) provides an exact reconstruction of the function *f* from its values . Comparing equation (4.1) with standard Lagrange interpolation, we see that the functions play a role similar to that of Lagrange interpolating polynomials and therefore we may think of equation (4.1) as an analogue of Lagrange interpolation on the sphere. Formulas of this type have appeared in Sloan (1995) and Sloan & Womersley (2000) using the standard quadratures for the sphere that are not invariant under the action of a discrete rotation group. Due to concentration of nodes near the poles and lack of invariance, the properties of representations obtained in these papers differ from that in equation (4.1).

As every function may be written as a linear combination of functions , we may use such linear combinations to represent all functions in . We avoid calling a basis because the number of these functions exceeds the dimension of the subspace ; however, effectively they play the same role as a basis for . We elaborate on this further below and note that an attempt to use the same number of functions *K*(** ω**·

*ω*_{j}) as the dimension of leads to ill-conditioned systems (Sloan & Womersley 2002; Fernandez 2005).

Let us now consider and evaluate equation (4.1) at the quadrature nodes . We obtain
where we multiplied both sides by . Denoting and , we obtain the matrix identity
4.2
where the matrix with entries *K*_{ij} is a projector with eigenvalues of either 1 or 0. As the number of nodes *M* for the projector on a subspace of dimension (*N*+1)^{2} is at least (2*N*+1)^{2}/3, it asymptotically exceeds the dimension of the subspace by a factor approximately equal to 4/3. This results in zero eigenvalues to match the dimension of and the range of . For example, the dimension of is 64, and we need a quadrature with the minimum of *M*=72 nodes so that has 64 eigenvalues equal to 1 and eight eigenvalues equal to 0. We note that, given the projector, the null space does not cause problems in numerical computations because it allows us to work on the range coinciding with . The near-optimal number of nodes does reduce the number of unknowns in the problem of scattered data interpolation on the sphere, a problem we plan to consider elsewhere.

In figure 4, we show a cross section of *K*(** ω**·

**e**

_{z}) through the

*z*-axis, with the maximum value of

*K*normalized to 1.

## 5. Algorithms associated with icosahedral grids

For efficient use of the quadratures developed in this paper, we need fast algorithms for evaluation of sums on these grids. Currently, we only have algorithms for this purpose and are working on the development of faster methods. We note that although there are algorithms for evaluation of spherical harmonics, for most of them the break-even point in comparison with the usual versions is quite high (Spotz & Swarztrauber 2001). Here, we describe algorithms for the numerical implementation of the representation (4.1). We consider two approaches: the first is based on the ideas in Jakob-Chien & Alpert (1997) and the second on the unequally spaced fast Fourier transform (USFFT; Dutt & Rokhlin 1993; Beylkin 1995; Lee & Greengard 2005).

### (a) Application of kernel *K* using the Christoffel–Darboux formula

The key observation in Jakob-Chien & Alpert (1997) is that the numerical evaluation of equation (4.1) may be accelerated by using the addition theorem and the Christoffel–Darboux formula, provided that it is implemented via a fast algorithm. The kernel *K* may be written as
5.1
where are the normalized associated Legendre functions (see §2) and . With an appropriate rotation of coordinates, the quadrature nodes can be seen to lie on planes of constant polar angle, with five points equally spaced in azimuthal angle on each plane. We write the coordinates of the quadrature nodes in plane *i* as , where *δ*_{i} is the relative shift in azimuth from plane *i*−1 to plane *i*. We label the planes such that *θ*_{i−1}<*θ*_{i} and set *δ*_{1}=0. Substituting equation (5.1) into equation (4.1), we need to evaluate the sum
5.2
where *N*_{planes} is the number of planes. Following the key idea in Jakob-Chien & Alpert (1997), we evaluate the summation over the index *i* in equation (5.2) using a fast algorithm (e.g. fast multipole method (FMM) as in Jakob-Chien & Alpert (1997)). However, the quadrature used in Jakob-Chien & Alpert (1997) is the standard quadrature with Gauss–Legendre nodes in the polar direction (and thus *N*_{planes}∼*N*) and equally spaced nodes in the azimuthal direction. Using the fast Fourier transform to evaluate along the azimuthal direction and the FMM to evaluate in the polar direction yields the overall complexity of . Unfortunately, in our case, the number of planes *N*_{planes}∼*N*^{2}, thus yielding only an algorithm.

### (b) Application of kernel *K* using the unequally spaced fast Fourier transform

As the kernel *K* may be written as in equation (2.8), it also has a finite Fourier series representation
Substituting this sum into equation (4.1), changing the order of summations and evaluating at arbitrary points on the sphere, , yields
5.3
where *M* is the total number of quadrature nodes. Recall that . The sum in equation (5.3) may be evaluated in operations using USFFT. We split equation (5.3) as
5.4
with
5.5
We may view equation (5.5) as evaluation of function *g* at points *n* *ω*_{j}, |*n*|≤*N* and (which may be interpreted as frequencies for the purposes of USFFT). As |*ω*_{i}|=1 in equation (5.5), fast evaluation of such sums is precisely the task for USFFT. Once we obtain values *g*(*n* *ω*_{j}), we compute the sum in equation (5.4) directly at the cost of .

## 6. Local and multi-resolution representations on the sphere

We expect that grids invariant under the icosahedral group (or other finite rotation groups) will play an important role in constructing efficient and practical local and multi-resolution representations of functions on the sphere. In this paper, we only point out that the rate of decay and the oscillations of the functions *K*(** ω**·

*ω*_{i}) may be altered by extending and bringing gradually to zero the spectrum of the kernel as a projector. We consider these results as preliminary because there might be representations with better localization.

If we were to measure the decay of a function on the sphere using variance defined as
6.1
where the mean is defined as
6.2
then we can show that, for large *N*, the variance of the kernel *K* decays only as .

To improve localization of the kernel, let us consider
6.3
where *p* is the over-sampling factor and the coefficients *a*_{n} are chosen to improve localization. Substituting equation (6.3) into equation (6.1) and minimizing the result with respect to *a*_{n}, we find that the resulting coefficients *a*_{n} decrease linearly with a particular (optimal) slope. Using these optimal coefficients, we achieve . As the total number of nodes is also proportional to *N*^{2}, this indicates that, for a given node, the number of neighbours needed to be taken into account to achieve a given accuracy remains constant as *N* becomes large.

We now show that may be used as a projector onto . For fixed ** ω**, we have , while . Thus, we obtain
Now, for we write
resulting in
so that is a projector onto . One of the benefits of using over

*K*is that we may consider fast algorithms exploiting the local nature of . We leave it as a subject for future research.

## 7. Conclusions

We introduce a general numerical method for constructing quadratures invariant under the icosahedral group that integrate rotationally invariant subspaces of . We note that, while the specific results presented here are for the icosahedral rotation group, quadratures invariant under other symmetry groups may be constructed following the same approach. Using these quadratures, we develop an exact representation of functions on rotationally invariant subspaces of , similar to Lagrange interpolation. We describe two algorithms of complexity for the numerical use of this representation. Furthermore, we develop a representation with localized functions and describe its properties.

In many ways, the results of this paper are just the first step in a programme to develop practical methods in a variety of applications that have to deal with the sphere. Further research should involve:

fast algorithms using these grids that have complexity, hopefully with a low break-even point in comparison with algorithms,

a greater variety of localized representations with associated fast algorithms,

multi-resolution representations that maintain a relation with the rotationally invariant subspaces, and

appropriate representations of integral and differential operators on the sphere and their embedding into a variety of applications.

## Acknowledgements

This research was partially supported by AFOSR grant FA9550-07-1-0135, NSF grant DMS-0612358 and DOE/ORNL grants 4000038129 and DE-FG02-03ER25583.

## Appendix A. Evaluation of the Jacobian

Here, we give details of computing the Jacobian of the nonlinear system (3.7). Let be spherical coordinates of a point on the unit sphere, which we may also write as a quaternion
Applying a rotation as quaternion multiplication, , we obtain a quaternion of the form
where the functions *f*, *g* and *h* describe the action of the group element. Hence, the rotated point has spherical coordinates
and
Thus, the spherical harmonic
after rotation, is written as
A 1
where is a normalization constant. We now compute the partial derivatives of equation (A 1) with respect to *θ* and *ϕ*. We obtain
and
Continuing with the second term, we have
and, furthermore,
As
and
we arrive at
and
We now calculate the factors , and
For the icosahedron group, these factors are of the form
where the numbers *α*_{1}, *β*_{1}, *γ*_{1}, *A*_{1}, *B*_{1}, *α*_{2}, *β*_{2}, *γ*_{2}, *A*_{2} and *B*_{2} are determined by direct computation using, for example, Mathematica and then tabulated.

## Footnotes

↵1 We note that an alternative formula for

*S*(*n*), eqn (17) in Sobolev's (1962) paper, is incorrect but does not affect anything else in his paper.- Received February 24, 2009.
- Accepted June 30, 2009.

- © 2009 The Royal Society