Quadrature rules on manifolds

Quadrature rules on manifolds play an important role for numerical integration with many applications in science. A quadrature rule on a manifold \( \mathcal M \subset \mathbb R^n \) with \(M\) quadrature points \( \boldsymbol p_1,\dots, \boldsymbol p_M \in \mathcal M \) and associated quadrature weights \( w_1,\dots,w_M \in \mathbb R\) approximates the integral of a continuous function \(f:\mathcal M \to \mathbb R\) over the manifold \(\mathcal M\) by a finite sum, i.e., \[ \int_{\mathcal M} f(\boldsymbol x) \,\mathrm d \boldsymbol x \approx \sum_{i=1}^M w_i f(\boldsymbol p_i).\]

Of particular importance are Gauß-type quadrature rules which integrate exactly all multivariate polynomials of maximal degree \(N\), i.e., \[ \int_{\mathcal M} f(\boldsymbol x) \,\mathrm d \boldsymbol x = \sum_{i=1}^M w_i f(\boldsymbol p_i), \qquad f(\boldsymbol x)=\prod_{i=1}^n x_i^{\alpha_i}, \; \sum_{i=1}^n \alpha_i \le N. \tag{1} \] If additionally all quadrature weights are equal, i.e., \(w_1=\dots=w_M= \frac1M \int_{\mathcal M} \mathrm d \boldsymbol x\), such quadrature rules are of Chebyshev-type. Gauß- and Chebyshev-type quadrature rules are called optimal if for prescribed degree of exactness \(N\) the identity (1) is achieved for the minimal number \(M\) of quadrature points.

The tables below list several quadrature rules for the two-dimensional sphere \[\mathrm S^2 := \left\{ \boldsymbol x \in \mathbb R^3 \;:\; \| \boldsymbol x \|_2 = 1 \right\}\] and the three-dimensional rotation group \[\mathrm{SO(3)}:=\left\{ \boldsymbol R \in \mathbb R^{3\times3} \;:\; \boldsymbol R^\top \boldsymbol R = \boldsymbol I,\; \det \boldsymbol R = 1 \right\}.\] This page was originally created by Manuel Gräf. The algorithms for the computation of such quadrature rules together with the general mathematical framework can be found in his PhD thesis [Gr13Diss]; the optimization software is available as the LORM library.

Putatively optimal quadrature rules on the sphere \(\mathrm S^2\)

The quadrature rules on the sphere \(\mathrm S^2\) have been computed numerically by minimizing the worst-case quadrature error \[ E_N( \boldsymbol p_1, w_1, \dots, \boldsymbol p_M, w_M ) := \left( \sum_{n=0}^N \sum_{k=-n}^n \left| \frac{1}{\sqrt{4\pi}} \delta_{0,n} - \sum_{i=1}^M w_i \overline{Y_n^k(\boldsymbol p_i)} \right|^2 \right)^{\frac12}\] over the product manifold \( (\mathrm S^2 \times \mathbb R)^M \), where \(\delta_{0,n}\) is the Kronecker delta and \(Y_n^k:\mathrm S^2 \to \mathbb C\) are the orthonormalized spherical harmonics.

For the efficient evaluation of the worst-case quadrature error \(E_N\) in the optimization routine we utilized the nonequispaced fast spherical Fourier transforms (NFSFT), which are implemented in the NFFT library.

Group symmetries: Quadrature rules on the sphere \(\mathrm S^2\) which are invariant under finite rotation groups \(\mathcal G \subset \mathrm{SO(3)}\) may be written in the form \[ \sum_{i=1}^M w_i f(\boldsymbol p_i) = \sum_{i=1}^{M_{\mathrm{gen}}} w_i \sum_{\boldsymbol p \in \mathcal G \boldsymbol p_i} f(\boldsymbol p),\qquad M= \sum_{i=1}^{M_{\mathrm{gen}}} \left|\mathcal G \boldsymbol p_i\right|, \] where the first \(M_{\mathrm{gen}}\) quadrature points are the generators of the disjoint orbits \( \mathcal G \boldsymbol p_i := \{ \boldsymbol G \boldsymbol p_i \;:\; \boldsymbol G \in \mathcal G\} \subset \mathrm S^2,\; i=1,\dots,M_{\mathrm{gen}}\).

The idea of taking group symmetries for quadrature rules into account dates back to Sobolev [So62] and McLaren [McL63]. In particular, we computed optimal Gauß- and Chebyshev-type quadrature rules which are invariant under the inversion group \(\overline{\mathrm C}_1\), the cyclic groups \(\mathrm C_k \; (k=1,\dots,5)\), the tetrahedral group \(\mathrm T\), the octahedral group \(\mathrm O\), or the icosahedral group \(\mathrm I\).

Tables: The listed quadrature rules have a worst-case quadrature error \(E_N < 10^{-12}\). The columns indicate the degree of exactness \(N\), the number of quadrature points \(M\), the symmetry group \(\mathcal G\), and provide a link to the text files containing the Cartesian coordinates of the quadrature points and weights. For Gauß-type quadrature rules the four Cartesian coordinates \((p_1,p_2,p_3,w)\) of a quadrature point \(\boldsymbol p:=(p_1,p_2,p_3)^\top \in \mathrm S^2\) with weight \(w\) are provided row-wise, whereas for Chebyshev-type quadrature rules the weight \(w\) is omitted. The first two numbers in the files give the number of coordinates (4 or 3) and the number \(M\) of quadrature points.

The accuracy of the computed quadrature rules can keep up with algebraically constructed quadrature rules. In fact, we recomputed with a precision of 11 digits the well-known Gauß-type quadrature rules with degree of exactness \(N = 1,2,3,5\) and the ones constructed in [Po94] with \(N=4,6,7,17\), in [Po95] with \(N=8\), in [Po98] with \(N=11\), and in [McL63] with \(N=9,14\).

Gauß-type quadratures on \(\mathrm S^2\)

Note that the set of quadrature points of a Chebyshev-type quadrature rule on the sphere \(\mathrm S^2\) with degree of exactness \(N\) is also known as spherical \(N\)-design. For further collections of spherical designs see the site of R. H. Hardin and N. J. A. Sloane under related websites.

Chebyshev-type quadratures on \(\mathrm S^2\)

(Non-optimal) Chebyshev-type quadratures on \(\mathrm S^2\) of high degree \(N\) (spherical \(N\)-designs)

The following table contains some examples of computed Chebyshev-type quadrature rules (spherical \(N\)-designs) for a high degree of exactness \(N\) and shows the practicability of the CG method proposed in [GrPo11] for such high-dimensional optimization problems. The text file for each point set consists of two columns which contain the spherical coordinates \((\varphi,\theta) \in [0,2\pi] \times [0,\pi) \) of the quadrature points \( \boldsymbol p := (\sin\theta \cos\varphi, \sin\theta \sin\varphi, \cos \theta)^{\top} \in \mathrm S^2 \).

Putatively optimal quadrature rules on the rotation group \(\mathrm{SO}(3)\)

The quadrature rules on the rotation group \(\mathrm{SO}(3)\) have been computed numerically by minimizing the worst-case quadrature error \[ E_N( \boldsymbol p_1, w_1, \dots, \boldsymbol p_M, w_M ) := \left( \sum_{n=0}^N \frac{2n+1}{8\pi^2} \sum_{k,l=-n}^n \left| \delta_{0,n} - \sum_{i=1}^M w_i \overline{D_{k,l}^n(\boldsymbol p_i)} \right|^2 \right)^{\frac12}\] over the product manifold \( (\mathrm{SO}(3) \times \mathbb R)^M \), where \(\delta_{0,n}\) is the Kronecker delta and \(D_{k,l}^n:\mathrm{SO}(3) \to \mathbb C\) are the Wigner D-functions.

For the efficient evaluation of the worst-case quadrature error \(E_N\) in the optimization routine we utilized the nonequispaced fast Fourier transforms on the rotation group (NFSOFT), which are implemented in the NFFT library.

Group symmetries: Quadrature rules which are invariant under finite rotation groups \( \mathcal H \subset \mathrm{SO(3)}\) and \(\mathcal G := \mathcal G_1 \times \mathcal G_2 \subset \mathrm{SO(3)} \times \mathrm{SO(3)}\) may be written in the form \[ \sum_{i=1}^M w_i f(\boldsymbol p_i) = \sum_{i=1}^{M_{\mathrm{gen}}} w_i \sum_{\boldsymbol p \in \mathcal H \boldsymbol p_i} f(\boldsymbol p), \qquad M=\sum_{i=1}^{M_{\mathrm{gen}}} \left| \mathcal H \boldsymbol p_i \right|,\] and \[\sum_{i=1}^M w_i f(\boldsymbol p_i) = \sum_{i=1}^{M_{\mathrm{gen}}} w_i \sum_{\boldsymbol p \in (\mathcal G_1 \times \mathcal G_2) \boldsymbol p_i} f(\boldsymbol p),\qquad M= \sum_{i=1}^{M_{\mathrm{gen}}} \left|(\mathcal G_1 \times \mathcal G_2) \boldsymbol p_i\right|,\] where the first \(M_{\mathrm{gen}}\) quadrature points are the generators of the disjoint orbits \( \mathcal H \boldsymbol p_i := \{ \boldsymbol H \boldsymbol p_i \;:\; \boldsymbol H \in \mathcal H \} \subset \mathrm{SO(3)},\; i=1,\dots, M_{\mathrm{gen}}\) and \( (\mathcal G_1 \times \mathcal G_2) \boldsymbol p_i := \{ \boldsymbol G \boldsymbol p_i \boldsymbol H^\top \;:\; \boldsymbol G \in \mathcal G_1, \boldsymbol H \in \mathcal G_2 \} \subset \mathrm{SO}(3),\; i=1,\dots,M_{\mathrm{gen}}\), respectively.

In particular, we computed optimal Gauß- and Chebyshev-type quadrature rules which are invariant under the tetrahedral group \(\mathrm T \times \mathrm C_1\), the octahedral group \(\mathrm O\times \mathrm C_1\), or the icosahedral group \(\mathrm I \times \mathrm C_1\).

Tables: The listed quadrature rules have a worst-case quadrature error \(E_N < 10^{-12}\). The columns indicate the degree of exactness \(N\), the number of quadrature points \(M\), the symmetry group \(\mathcal G\), and provide a link to the text files containing the Cartesian coordinates of the quadrature points and weights. For Gauß-type quadrature rules the ten Cartesian coordinates \((p_{11},p_{12},p_{13},p_{21},p_{22},p_{23},p_{31},p_{32},p_{33},w)\) of a quadrature point \(\boldsymbol p := (p_{ij})_{i,j=1}^3 \in \mathrm{SO(3)}\) with weight \(w\) are provided row-wise, whereas for Chebyshev-type quadrature rules the weight \(w\) is omitted. The first two numbers in the files give the number of coordinates (10 or 9) and the number \(M\) of quadrature points.

Note that quadrature rules with degree of exactness \(N\) on the rotation group \(\mathrm{SO}(3)\) are in one-to-one correspondence to quadrature rules with degree of exactness \(2N+1\) and \(2M\) quadrature points on the three-dimensional sphere \( \mathrm S^3 := \{ \boldsymbol x \in \mathbb R^4 \;:\; \|\boldsymbol x \|_2 = 1 \} \) which are invariant under the inversion group.

For the illustration of the quadrature points we convert rotation matrices \( \boldsymbol R \in \mathrm{SO}(3) \) to unit quaternions \( \boldsymbol q := (q_0,q_1,q_2,q_3) \in \mathrm S^3 \) with real part \(q_0 \ge 0\) and project them stereographically by \( \boldsymbol q \mapsto \frac{1}{1+q_0}(q_1,q_2,q_3) \in \mathbb R^3 \) into the unit ball \( \{ \boldsymbol x \in \mathbb R^3 \;:\; \| \boldsymbol x \|_2 \le 1 \} \). Since antipodal quaternions \( \pm \boldsymbol q \) correspond to exactly one rotation matrix, we plot for symmetry reasons both antipodal quaternions if \( q_0 = 0 \).

Gauß-type quadratures on \(\mathrm{SO}(3)\)

Note that the set of quadrature points of a Chebyshev-type quadrature rule on the rotation group \(\mathrm{SO}(3)\) with degree of exactness \(N\) can be considered as spherical \((2N+1)\)-design on the sphere \(\mathrm S^3\).

Chebyshev-type quadratures on \(\mathrm{SO}(3)\)

C++ software library

The optimization algorithms are implemented in a C++ library for optimization on Riemannian manifolds (LORM), which is free software licensed under the MPL 2.0.

• lorm-1.0-alpha.tar.gz (August 22, 2013)

Related websites

• R. H. Hardin and N. J. A. Sloane: Spherical designs
• R. Womersley: Interpolation and cubature on the sphere
• Ch. F. F. Karney: Nearly optimal coverings of orientation space
• Wikipedia: Lebedev quadrature

Publications

[Gr13Diss]

M. Gräf.
Efficient algorithms for the computation of optimal quadrature points on Riemannian manifolds.
PhD thesis, Chemnitz University of Technology, Universitätsverlag Chemnitz, ISBN 978-3-941003-89-7, 2013. pdf

[GrPo11]

M. Gräf and D. Potts.
On the computation of spherical designs by a new optimization approach based on fast spherical Fourier transforms.
Numer. Math. 119, 699–724, 2011. pdf

Selected references

[McL63]

A. D. McLaren.
Optimal numerical integration on a sphere.
Math. Comput. 17, 361–383, 1963. JSTOR

[Po94]

A. S. Popov.
Cubature formulae for a sphere invariant under cyclic rotation groups.
Russian J. Numer. Anal. Math. Modelling 9, 535–546, 1994. doi

[Po95]

A. S. Popov.
Cubature formulae for a sphere which are invariant with respect to the tetrahedral group.
Comput. Math. Math. Phys. 35, 369–374, 1995.

[Po98]

A. S. Popov.
Cubature formulas on a sphere that are invariant with respect to octahedron rotation groups.
Comput. Math. Math. Phys. 38, 30–37, 1998.

[So62]

S. L. Sobolev.
Cubature formulas on the sphere invariant under finite groups of rotations.
Dokl. Akad. Nauk SSSR 146, 310–313, 1962. doi