Spectral Decomposition of Discrepancy Kernels on the Euclidean Ball, the Special Orthogonal Group, and the Grassmannian Manifold

To numerically approximate Borel probability measures by finite atomic measures, we study the spectral decomposition of discrepancy kernels when restricted to compact subsets of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^d$$\end{document}Rd. For restrictions to the Euclidean ball in odd dimensions, to the rotation group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{SO}(3)$$\end{document}SO(3), and to the Grassmannian manifold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {G}_{2,4}$$\end{document}G2,4, we compute the kernels’ Fourier coefficients and determine their asymptotics. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}L2-discrepancy is then expressed in the Fourier domain that enables efficient numerical minimization based on the nonequispaced fast Fourier transform. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{SO}(3)$$\end{document}SO(3), the nonequispaced fast Fourier transform is publicly available, and, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {G}_{2,4}$$\end{document}G2,4, the transform is derived here. We also provide numerical experiments for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{SO}(3)$$\end{document}SO(3) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {G}_{2,4}$$\end{document}G2,4.


Mathematics Subject Classification
cf. [40,42,43], see Sect. 2 for explicit examples 1 . For fixed n ∈ N, we aim to minimize D β (μ, ν n ) among all n-point sets {x 1 , . . . , x n } ⊂ R d . The present manuscript is concerned with discretizations of (1.2) that facilitate numerical minimization. The associated discrepancy kernel K β : R d × R d → R is defined by yield that (1.2) is identical to If a compact set X ⊂ R d is known in advance such that supp(μ) ⊂ X, then we shall restrict the minimization to {x 1 , . . . , x n } ⊂ X, so that only the restricted kernel K β | X×X matters. By endowing X with a finite Borel measure σ X having full support, Mercer's Theorem yields an orthonormal basis {φ l } ∞ l=0 for L 2 (X, σ X ) and coefficients (a l ) ∞ l=0 such that the spectral decomposition a l φ l (x)φ l (y), x, y ∈ X, (1.5) holds with absolute and uniform convergence. We call (a l ) ∞ l=0 the Fourier coefficients of the kernel K β | X×X . If supp(μ), supp(ν n ) ⊂ X, then the Fourier expansion of the φ l (x j ), (1.6) where the Fourier coefficientsμ l andν n,l of the measures μ and ν n , respectively, are well-defined if a l = 0. Truncation of the discretization (1.6) enables the use of the nonequispaced fast Fourier transform, thereby offering more efficient minimization of D β (μ, ν n ), cf. [31,33]. Thus, we aim to A) compute (a l ) ∞ l=0 and (φ l ) ∞ l=0 in the Fourier expansion (1.5) of K β | X×X .
Fourier decay properties generally quantify Sobolev smoothness. To accomplish (B), we aim to determine the asymptotics of K β | X×X 's Fourier coefficients (a l ) ∞ l=0 in (1.5).
For X = S d−1 and a particular choice of β, the kernel K β | S d−1 ×S d−1 essentially coincides with the Euclidean distance, see [12,13]. The Fourier expansion is determined in [10], and the decay of the Fourier coefficients yields that K β | S d−1 ×S d−1 reproduces the Sobolev space H β (S d−1 ) = H d This manuscript is dedicated to derive analogous results for other compact sets X. We focus on the unit ball, the special orthogonal group, and the Grassmannian manifold, We achieve goal (A) for X = B d with odd d. Both goals, (A) and (B), are achieved for SO (3) and G 2,4 . We also provide numerical experiments. For SO (3), the computations are based on the nonequispaced fast Fourier transform designed in [32,44]. For G 2,4 , we derive the nonequispaced fast Fourier transform by parametrization via the double covering S 2 × S 2 and developing the respective transform there.
For SO(d) and G k,d with fixed k and d, in principle, one could still be able to compute the Fourier expansion in goal (A). However, one may be faced with rather complicated expressions. In our present computations for SO (3) and G 2,4 , the structural relations to the unit sphere enabled the use of Chebychev and Legendre polynomials, which reduced the complexity. Nonetheless, we do accomplish (B) for the general cases SO(d) and G k,d .

Two Introductory Examples
We first present a well-known elementary example on the interval [0, s], for which both aims A) and B) are achieved. Second, to support our perspective on discrepancy, we prove that the so-called Askey function is a discrepancy kernel of the form (1.3).
where the inner product between f and g is given by f , The associated discrepancy kernel is In order to additionally integrate over r , recall the (generalized) hypergeometric functions where f 1 , . . . , f k , g 1 , . . . , g l , z ∈ R and ( f ) n := f · ( f + 1) · · · ( f + n − 1) denotes the Pochhammer symbol with ( f ) 0 := 1. We consider G d : [0, ∞) → R given by Since d is odd, either d+1 4 or d−1 4 is a natural number, so that the series terminates and G d is a polynomial in r 2 on [0, 1]. By integration with respect to G d , we obtain the L 2 -discrepancy and the associated discrepancy kernel respectively. It turns out that K d coincides with Askey's function.
The proof is presented in Appendix A. Provided that d ≥ 3, Askey's kernel function reproduces the Sobolev space H d+1 2 (R d ) with an equivalent norm, see [49].

The Distance Kernel on S d−1
This section is dedicated to recall results on discrepancy kernels on the sphere S d−1 ⊂ R d , for d ≥ 2, from [12,13,31,46] that shall guide our subsequent investigations. Denote the geodesic ball of radius r centered at z ∈ S d−1 by and endow [0, π] with the weighted Lebesgue measure sin(r )dr , whereas S d−1 carries the normalized, orthogonal invariant surface measure The associated discrepancy kernel is According to [12,13,31], see also [2], K β d satisfies is the reproducing kernel Hilbert space associated with the reproducing kernel The coefficients in the Fourier expansion satisfy |c m | ∼ m −d , cf. [12]. This is the same asymptotics as the coefficients in (3.3) In order to determine the Fourier coefficients of kernels on the sphere that are polynomial in x − y , such as K β d | S d−1 ×S d−1 , we require the Fourier coefficients of the monomial terms x − y p . For any p ∈ N, the Fourier expansion (3.5) Note that (3.5) is well-defined for the entire range p > −(d − 1) and p is not required to be an integer. For p > 0, the following proposition is essentially due to [10], see also [12,14]. Simple continuation arguments cover the full range of p, and the asymptotics . (3.6) and the series (3.4) terminates if p ∈ 2N.
For p ∈ 2N, the term (− p 2 ) is not well-defined and (3.6) is to be understood with the convention It is noteworthy that the kernel K d,r in (2.2) for d = 3 is a discrepancy kernel that does not generate a Sobolev space on R d but its restriction does. The proof of the following proposition is presented in Appendix B. Proposition 3.2 Let r ≥ 1. The reproducing kernel Hilbert space of K 3,r , given by (2.2) with d = 3, is continuously embedded into H 2 (R 3 ), but the reverse embedding does not hold. In contrast, K 3,r | S 2 ×S 2 reproduces H For supp(μ), supp(ν n ) ⊂ S 2 , the L 2 -discrepancy (1.6) for K β 3 with X = S 2 becomes (3.8) Fig. 1 The target measure μ is supported on two circles on the sphere S 2 with weight ratio 9/1. Numerical minimization of (3.9) splits 50 points into 45 points equally distributed on one and 5 points on the other circle whereμ m l denotes the Fourier coefficient of μ with respect to Y m l , cf. (1.6). By truncating this series, the nonequispaced fast Fourier transform on S 2 , cf. [33,39,41], enables efficient minimization of

Discrepancy Kernels on Compact Sets
Here we discuss discrepancy kernels that extend the kernels of the previous section in a natural way. For d ≥ 1, let us define the half-space For fixed s > 0, we consider the mapping h : The associated discrepancy kernel is with K β d as in (3.1). In contrast to K β d , the kernel K β d,s is not identically zero outside of S d−1 × S d−1 and makes also sense for d = 1.

Its reproducing kernel Hilbert space is
where the inner product between f and g is given by The identity (4.2) for x, y ∈ S d−1 with s = 1 has been established in [13], see also (3.2). Essentially, the same proof still works for the more general situation. Theorem 4.3 provides a simple form of K β d,s | X×X with X = B d s , which may facilitate further computations. An immediate consequence is x − y , for x, y ∈ B d s .

The Euclidean Ball B d
This section is dedicated to derive the Fourier expansion of the discrepancy kernel K β d,1 in (4.1) on B d . Proposition 4.2 has covered d = 1, and we now derive the spectral decomposition of The case d = 3 with p = −1 is discussed in [37]. Let {C α m : m ∈ N, α > −1/2} denote the family of Gegenbauer polynomials with the standard normalization , α = 0.
By α = d 2 − 1, the addition theorem for spherical harmonics yields For d ≥ 3 and arbitrary real p > 1 − d, we deduce from [18] that For x = 0 or y = 0, the right-hand side of (5.3) is well-defined by analytic continuation.
Using the addition formula (5.1) we obtain The Fourier expansion of K d, p m with respect to the measure r d−1 dr satisfies as well as the integral operator in Mercer's theorem induced by K d, p , so that direct computations yield . This leads to the Fourier expansion Thus, the original problem is reduced to the spectral decomposition of the sequence of kernels K m, j and eigenfunctions ϕ d, p m, j . We now specify these eigenvalues and eigenfunctions. In the following theorem, J ν denotes the Bessel function of the first kind of order ν and ζ k := e 2π i/k is the k-th root of unity.

Remark 5.2 Computer experiments seem to indicate that the nullspace of
In that case, the function where A [1, ] (ω) denotes the (1, ) minor of A(ω), spans the eigenspace associated with λ.
Appendix D is dedicated to the proof of Theorem 5.1. The proof reveals strong ties with polyharmonic operators on the unit ball and higher order differential operators on the interval [0, 1]. We refer to [1] for structurally related spectral decompositions of polyharmonic operators on [0, 1] with homogeneous Neumann boundary conditions.
The corresponding eigenspaces are 1-dimensional with the representative The formulas in Corollary 5.3 are derived from Theorem 5.

The Rotation Group SO(3)
In this section we derive the Fourier expansion of the discrepancy kernel on the special orthogonal group SO(3). The eigenfunctions turn out to be classical functions but the coefficients and their decay rates need to be determined. We also provide numerical experiments by using the nonequispaced fast Fourier transform on SO(3).

Fourier Expansion on SO(3)
By identifying R d×d with R d 2 , Theorem 4.3 applies to subsets of R d×d endowed with the trace inner product x, y F := trace(x y), x, y ∈ R d×d , and the induced Frobenius norm · F on R d×d . In this way, SO(3) is contained in , and it is natural to consider s = √ 3. We endow SO (3) (3)), cf. [48]. For p > 0, the Fourier expansion holds and, analogous to (3.5), the coefficients are We now compute these coefficients for the entire range p > −3.
The proof is given in Appendix E. For p ∈ 2N, we again apply the convention (3)) is the reproducing kernel Hilbert space associated with the reproducing kernel The choice p = 1 in Proposition 6.1 implies that the kernel , with an equivalent norm. Indeed, Theorem 4.3 We use the standard parametrization of SO(3) by S 3 via unit quaternions, which is then mapped into R 3 by stereographic projection and R 3 is parametrized by The target measure μ is supported on two disjoint parts with weight ratio 9/1 colored in darker blue by the cylindrical surface and the great circle. Numerical minimization of (6.5) splits 30 points in SO(3) into 27 points on the inner surface and 3 points on the great circle. We plotted 6 points on the great circle but antipodal points correspond to the same point in SO (3) Since

The Grassmannian G 2,4
First, the Fourier expansion of the discrepancy kernel on G 2,4 is computed. To prepare for developing the nonequispaced fast Fourier transform on G 2,4 , we then explicitly parametrize the Grassmannian G 2,4 by its double covering S 2 ×S 2 . Next, we derive the nonequispaced fast Fourier transform on S 2 × S 2 and provide numerical minimization experiments on G 2,4 .
The proof of this theorem is contained in Section F.1 of Appendix F. If p ∈ 2N, then a λ ( p, G 2,4 ) = 0 for all |λ| > p 2 . For τ > 2, the Sobolev space H τ (G 2,4 ) is the reproducing kernel Hilbert space with associated reproducing kernel cf. [11,16]. Since the coefficients in (7.7) behave asymptotically as λ −2τ , the choice p = 1 in Theorem 7.1 implies that the kernel K β 16 with an equivalent norm.

Parametrization of G 2,4 by (S 2 × S 2 ) ±1
To derive the nonequispaced fast Fourier transform on G 2,4 , we shall first explicitly construct the parametrization of G 2,4 by its double covering S 2 × S 2 . We denote the d × d-identity matrix by I d , and the cross-product between two vectors x, y ∈ S 2 is denoted by x × y ∈ R 3 . The mapping P : S 2 × S 2 → G 2,4 given by is surjective and, for all x, y, u, v ∈ S 2 , and direct computations lead to The right-hand side determines x and y up to the ambiguity (7.9). Under the Frobenius norm, P is distance preserving in the sense The latter follows from (F.20) in Lemma F.6 in Section F.2 of Appendix F. We shall now check how the spherical harmonics Y m l on S 2 relate to the eigenfunctions ϕ λ,l ∈ H λ (G 2,4 ) of the Laplace-Beltrami operator on G 2,4 , cf. (7.2). The functions Y m,n k,l : G 2,4 → C given by Y m,n k,l (P(x, y)) := Y m k (x) · Y n l (y) (7.13) are well-defined for m + n ∈ 2N, the latter taking into account the ambiguity (7.9). (7.14) The proof is presented at the end of Section F.2 of Appendix F. Note that the geodesic distance on G 2,4 is dist G 2,4 (P, Q) = √ 2 θ 2 1 + θ 2 2 , where θ 1 , θ 2 ∈ [0, π/2] are the principal angles determined by the two largest eigenvalues cos 2 (θ 1 ) and cos 2 (θ 2 ) of the matrix P Q. Aside from (7.12), P is also distance-preserving with respect to the respective geodesic distances, i.e., 4 The geodesic distance on S 2 induces the geodesic distance on S 2 × S 2 ±1 by dist 2 This equality follows from (F.23) in Lemma F.6 in the appendix via further direct calculations.
The identity (7.14) provides explicit expressions for the orthonormal basis of H λ (G 2,4 ) that is used to construct the reproducing kernel Q λ in (7.2). It also provides a fast Fourier transform on G 2,4 from the respective transform on S 2 × S 2 that is developed in the subsequent section.

Nonequispaced Fast Fourier Transform on G 2,4
The nonequispaced fast (spherical) Fourier transform on S 2 has been developed in [39,41] under the acronym nfsft. Here, we shall derive the analogous transform on S 2 × S 2 , which induces the nonequispaced fast Fourier transform on G 2,4 via the mapping P and (7.14) with (7.13).
at n scattered locations (x j , y j ) n j=1 ⊂ S 2 × S 2 . Direct evaluation of (7.15) leads to O(nM 4 ) operations. We shall now derive an approximative algorithm that is more efficient for n M. By following the ideas in [39,41], switching to spherical coordinates reveals that (7.15) is a 4-dimensional trigonometric polynomial. This enables the use of the 4dimensional nonequispaced fast Fourier transform nfft to significantly reduce the complexity. In spherical coordinates the spherical harmonics are trigonometric polynomials such that where 0 ≤ θ ≤ π , 0 ≤ ϕ ≤ 2π , and c m k,k ∈ C are suitable coefficients that we assume to be given or precomputed. Thus, for x = z(θ 1 , ϕ 1 ) and y = z(θ 2 , ϕ 2 ), there are coefficients b k,l k ,l ∈ C such that

Numerical example on G 2,4
By Theorem 7.1, we can calculate the coefficients of the kernel expansion The eigenfunctions ϕ λ,l are given by the tensor products of spherical harmonics in (7.13), cf. Theorem 7.2. For supp(μ), supp(ν n ) ⊂ G 2,4 , the L 2 -discrepancy (1.6) of the kernel K β 16, whereμ λ,l is the Fourier coefficient of μ with respect to ϕ λ,l , cf. (1.6). Let us consider μ = σ G 2,4 . According to [11] (see also [15,16]), the lower bound cf. [11,25]. Note that we can efficiently solve the least squares minimization (7.19) by using the nonequispaced fast Fourier transform on G 2,4 derived from the nonequispaced fast Fourier transform on S 2 × S 2 of Sect. 7.3 and applying Theorem 7.2. Figure 3 shows logarithmic plots of the number of points versus the L 2 -discrepancy. We observe a line with slope −5/4 as predicted by (7.20).

Funding Open access funding provided by Austrian Science Fund (FWF).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Proof of Theorem 2.1 Let h d : [0, ∞) → R denote Euclid's hat function given by
is the indicator function of B d 1/2 . For x, y ∈ R d and t = x − y ≤ 1, we derive where the last equality is due to partial integration. An explicit expression for h d is stated in [29,Equation (11)], so that we obtain .
We now compare coefficients of powers of t with those of the polynomial (1 − t) d+1 2 . In order to check the coefficient of t 2l , we first observe where the last equality makes use of Gauss' Theorem for the hypergeometric series evaluated at 1 (cf.
Hence, the coefficients for even powers of t match.
Thus, the coefficient of t 2 j+1 in K d (x, y) is nonzero if and only if j ≤ d−1 4 . Moreover, it is given by Further computations using the duplication formula eventually lead to − d+1 2 2l+1 . The case d+1 4 ∈ N is checked analogously. The observation K d (x, y) = 0 for x − y > 1 concludes the proof.

Appendix B: Proofs for Section 3
Proof of Proposition 3.2 By expressing the R d -Fourier transform 1 B d r (ξ ) in terms of the Bessel function of the first kind of order d/2 and using its asymptotics, we deduce so that K 3,r | S 2 ×S 2 is a polynomial of degree 3 in x − y . Its Fourier coefficients (a m ) m∈N are linear combinations of the Fourier coefficients of the monomial terms, so that (3.7) implies a m m −3 .
We have checked that there are no cancelations in these linear combinations. Therefore, the asymptotics (3.7) also imply the associated bound from below, which leads to a m ∼ m −3 . Thus, K 3,r | S 2 ×S 2 reproduces the Sobolev space H

Appendix C: Proofs for Section 4
Proof of Proposition 4.2 Let us consider s = 1. The general case follows from rescaling. In order to determine the Fourier expansion of K β 1,1 | [−1,1]×[−1,1] , we need to determine the eigenfunctions with respect to the positive eigenvalues of the integral operator For λ > 0, the equation T φ = λφ is equivalent to the associated Sturm-Liouville eigenvalue problem. Indeed, differentiating twice on both sides and carrying out a short calculation, we arrive at the second order homogeneous differential equation ). Direct calculations using the boundary conditions determine c 1 , c 2 and λ, so that Mercer's Theorem and normalization of the eigenfunctions provide the claimed Fourier expansion of the kernel.
Computations analogous to [43, Section 9.5.5] and [22] verify the claimed form of the reproducing kernel Hilbert space.

Appendix D: Proofs for Section 5
The two linearly independent eigenfunctions of the differential operator

Proof of Lemma D.1 Up to a constant depending on d and p, K d, p is the Green's function of the polyharmonic equation d+ p
2 u = f on B d with certain nonlocal boundary conditions, cf. [38]. In particular and by specifying the constant, one deduces that any eigenfunction of T d, p in (5.5) with eigenvalueλ = 0 is an eigenfunction of

The Laplacian in polar coordinates is
The decomposition (5.7) (see also (5.4)) yields The linearly independent eigenfunctions of D where and we take As an eigenfunction of a positive integer power of the Laplacian, φ d, p m, j,l is real analytic on B d , cf. [5,36]. Hence, the radial part ϕ d, p m, j,l must be an analytic function on [0, 1] with even or odd parity for m even or odd, respectively, cf. [9]. The functions Y d,ω m do not have matching parity, which concludes the proof.
In order to identify the eigenvalues and the linear combination in Lemma D. 1 Proof of Lemma D. 2 The idea of the proof is to use the series expansion of the Bessel functions, apply the integral operator to each term, and eventually recover the righthand side of (D.4). We shall provide the skeleton of the proof and omit some lengthy computations.
Let d ≥ 1, p ≥ 1 − d and r ≥ s. If d + p is even, direct computations yield We obtain that For α > 0, k ∈ N 0 , ω ∈ C and r > 0, integration of each term of the power series of the Bessel function eventually yields which follows from direct computations and for α ∈ R, k ∈ N 0 , and n ∈ N with α / eventually lead to the claimed equality (D.4).
Both hypergeometric series can be evaluated by means of Gauss' Theorem. Thus, we obtain Now we reverse the order of summation in the second sum, i.e., we replace k by l − k there. Then both sums can be conveniently put together to yield The hypergeometric function can be evaluated by means of the classical terminating very-well-poised 5 F 4 -summation (cf. [47,Equation (2.3.4.6); Appendix (III.13)]). After some simplification one arrives at the desired expression (−1) m /(m! (m + α + 1)). If m < l, then the first sum in (D.8) does not contribute anything because of the term (m − l − 1)! in the denominator. In the second sum, the summation over i may be started at i = l − m, which can be evaluated by means of the binomial theorem. The result is zero except if k = l. Again, in the end one obtains (−1) m /(m! (m + α + 1)).
For i = 1, . . . , d+ p 2 , the hypergeometric functions in F are polynomials of exact degree d + p − 2i and thus linearly independent. Hence, for c = 0, the right-hand side of (D.9) vanishes if and only if A(ω) is singular and c is in its nullspace.
where s = arccos( x, y ). For d = 2, the addition theorem yields where T m are the Chebyshev polynomials of the first kind, i.e., T m cos(s) = cos(ms), for s ∈ [0, π]. The relation (3.4) for d = 2 with (E.1) and (E.2) leads to We now switch to SO (3). For x, y ∈ SO(3), the relation The Legendre polynomials satisfy T m = 1 2 (C 1 m − C 1 m−2 ) with T 0 = C 1 0 and C 0 −1 = 0, for m ∈ N, so that we obtain We derive (6.3) by calculating the differences 2 see [48] for (E.5). Analytic continuation and well-known relations for the gamma function cover the remaining values of p > −3.

Standard calculations yield
(m− p 2 ) (m+ p Appendix F: Proofs for Section 7

Symmetry arguments yield
The series expansion of (1 − x y) −q converges absolutely for q < 2, and the Legendre polynomials C 1 2 m are orthogonal to the monomials x k for k < m or m ≡ k mod 2. Therefore, the orthogonality relations yield The use of the generating function and further calculations lead to By reordering and making use of n! (n + 1) m−n+2k = (m + 2k)!, which cancels the identical term in the denominator, we are led to where the equality in the last line is due to the identities The relation ( 3 2 ) n (n + 3 2 ) m−n 2 = ( 3 2 ) m+n 2 and our choices m = |λ|, n = λ 1 − λ 2 yield ) λ , so that q = −p/2 and the definition of the 4 F 3 hypergeometric series conclude the proof of (7.5).
Proof of (7.6): The proof of the decay property for a λ ( p, G 2,4 ) requires some preparation and auxiliary results. For p / ∈ 2N, direct calculations yield In the following, we treat the case λ 1 = an and λ 2 = (1 − a)n for n ∈ N and n → ∞, with an arbitrary a ∈ [ 1 2 , 1], so that |λ| = n. Summarizing the proof of (7.6), we shall first verify that the summand in (F.4), as a sequence in k, is unimodal, i.e., it first increases until it has reached its maximum and then decreases, see Lemma F.2. Second, approximation of S(n, k) for k ≥ εn 2 with the help of Stirling's formula, where ε > 0 but fixed, leads to an asymptotic formula for k≥εn 2 S(n, k); see Lemma F.3. Third, we let ε → 0 to obtain the asymptotic behavior of the full sum k≥0 S(n, k) for n → ∞; see (F.8). If the result is substituted in (F.4), the claimed decay in (7.6) follows immediately upon observing To start with, we may consider S(n, k) as a function of real k.
Lemma F.1 Given n ∈ N, let k 0 = k 0 (n) ∈ (0, ∞) be such that ∂ ∂k S(n, k 0 ) = 0. Then . Then We postpone the proofs of Lemmas F.1, F.2, and F.3, and discuss their consequences first. If we let ε → 0 and substitute t = p+6 2(1+x) in the above integral, then the integral definition of the gamma function yields Lemma F.3 and (F.7) provide the asymptotic lower bound on k≥0 S(n, k) of the following two-sided claim: (1)) , for n → ∞. which is due to Lemma F.2, saying that S(n, k) grows until its maximum at k = k 0 ∼ An 2 . By Lemma F.3 and letting ε → 0, we obtain the upper bound in (F.8).
To complete the proof of (7.6) in Theorem 7.1, it remains to prove Lemmas F.1, F.2, and F.3.

Proof of Lemma F.1
The condition 0 = ∂ ∂k S(n, k 0 ) implies 0 = ∂ ∂k log S(n, k 0 ) . By using the digamma function ψ(z), the logarithmic derivative of (F.5) can be written as For k = o(n), the above expression is certainly positive for large n, hence nonzero. If the order of magnitude of k is at least the one of n (in symbols, n = O(k)), then we may estimate the logarithmic derivative by ∂ ∂k log S(n, k) = log n 2 where the asymptotics ψ(x) = log x − 1 2x + O x −2 , for x → ∞, cf. [24, Equation 1.18 (7)]), were used. By applying the exponential function on both sides of ∂ ∂k log S(n, k 0 ) = 0, together with the above estimation for ∂ ∂k log S(n, k) we obtain where Q(n, k) := .
so that we may assume that k 0 is of larger asymptotic order of magnitude than n. Define P(n, k) by Q(n, k) = 1 − P(n, k) n + 3 2 + k an + 3 2 + k (1 − a)n + 1 + k (k + 1) . (F.11) For k of larger asymptotic order of magnitude than n, the asymptotics of the digamma function implies that the error term O(n −1 ) in (F.9) and (F.10) may be replaced by O(k −2 ) and hence by o(k −1 ). The relations (F.10) and (F.11) lead to P(n, k 0 ) = o(k 3 0 ). Since direct computations yield asymptotically leading terms in P(n, k) must cancel each other. We have already excluded k 0 ∼ cn, so that we now consider k 0 ∼ cn 2 for an appropriate constant c. The terms k 3 and k 2 n 2 in P(n, k) must cancel each other, so that The solution c = 1 p+6 (2a 2 − 2a + 1) yields the leading term in (F.6). In order to derive the O(n) term in (F.6), we have to perform "bootstrap", i.e., we substitute k 0 = cn 2 +k 1 in (F.12) and apply analogous arguments to eventually conclude k 1 = O(n).
We are looking for P : Theorem F. 4 There are exactly two mappings S 2 × S 2 → G 2,4 satisfying (F.17). One is P as in (7.8), and the other is I 4 − P. In particular, P is surjective and, for all x, y, u, v ∈ S 2 , Proof of Theorem F. 4 The identity (F.17) for the specific choice of P is verified by expanding both sides of the equality and comparing the polynomial expressions. We omit the straightforward but lengthy computation.
Since the conjugate action of SO(4) on G 2,4 is transitive, the identity (F.17) also implies surjectivity. Since left and right eigenspaces of L(P(x, y)) in (7.10) and (7.11) are uniquely determined, we deduce that (F.18) holds.
Remark F. 5 The Grassmannian G 2,4 has been parametrized in [20] by means of angles but the explicit identity (F.17) that steers our present approach has not been considered there. Our parametrization is compatible with the one used in [19].
The following properties of P are useful.
Proof of Lemma F. 6 The product measure σ S 2 ⊗ σ S 2 is SO(3) × SO(3) invariant. According to (F.17), the pushforward measure of σ S 2 ⊗ σ S 2 on S 2 × S 2 under P is SO(4) invariant, so that the uniqueness of the Haar measure implies the first claim of the lemma. Each of the remaining claims is first proved for P as defined in Proposition F.4 and then argued that it also holds for I 4 − P. The identity (F.20) is easily observed for u, v = e 1 first and then (F.17) yields the general case. According to I 4 − P(x, y), I 4 − P(u, v) F = P(x, y), P(u, v) F , the identity (F.20) also holds for I 4 − P. The statement (F.21) follows from expanding both sides of the equality and comparing polynomial expressions in x and y. If it holds for P, then it must also hold for I 4 − P. One directly calculates (F.22).
In order to check (F.23), we first recognize that principal angles between I −P(x, y) and I − P(u, v) coincide with the ones between P(x, y) and P (u, v). By ξ ± as in (F.2) and as at the beginning of the proof of Theorem 7.1, we observe trace(P(x, y)P(u, v)) = 1 + ξ + ξ − . Hence, (F.20) leads to ξ + ξ − = x, u y, v . Theorem 7.2 implies that Q λ in (7.2) satisfies Q λ (P(x, y), P(u, v)) = c λ 2 C Comparison of this identity with (F.1) and few further calculations eventually lead to (F.23).
The unit quaternions provide a group structure on S 3 , so that the mapping a → S a between S 3 ±1 and SO(3) as well as (a, b) → L a R b between (S 3 × S 3 ) ±1 and SO(4) become group isomorphisms. Their combination induces a group isomorphism between SO(3) × SO(3) and SO(4) ±1 . Condition (F.17) requires that the respective actions on S 2 × S 2 and G 2,4 commute with P. Hence, the induced pullback is an intertwining isomorphism. In particular, P * maps one irreducible subspace into the other. Thus, the irreducible decomposition of L 2 (G 2,4 ) under the action of SO (4) is where m λ = λ 1 + λ 2 and n λ = λ 1 − λ 2 . By considering L(P) = x y , L(P)·L(P) = x x , and L(P) ·L(P) = yy , we observe that the homogeneous polynomials x i y j and x i x j , y i y j , for i, j = 1, . . . , 3, can be written as homogeneous polynomials in the matrix entries of P ∈ G 2,4 of degree 1 and 2, respectively. The monomial x α y β with |α| = m λ and |β| = n λ is composed of n λ factors of the form x i y j and (m λ − n λ )/2 factors of the form x i x j . Thus, x α y β is a homogeneous polynomial of degree n λ + 2((m λ − n λ )/2) = m λ in the matrix entries of P. We deduce for given coefficients (b i ) M i=1 ⊂ C and locations (x i , y i ) M i=1 ⊂ S 2 × S 2 . Similar arguments as above and in Sect. 7.3 provide a fast evaluation of (F.28) by using the adjoint nonequispaced fast Fourier transform adjoint nfft, cf. [39,41].