{'fraction_non_alphanumeric': 0.0642797734563745, 'fraction_numerical': 0.02277317913580736, 'mean_word_length': 3.806777079763944, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasi Monte Carlo integration and kernel-based function approximation on Grassmannians

Numerical integration and function approximation on compact Riemannian manifolds based on eigenfunctions of the Laplace-Beltrami operator have been widely studied in the recent literature. The standard example in numerical experiments is the Euclidean sphere. Here, we derive numerically feasible expressions for the approximation schemes on the Grassmannian manifold, and we present the associated numerical experiments on the Grassmannian. Indeed, our experiments illustrate and match the corresponding theoretical results in the literature.


Introduction
The present paper is dedicated to numerical experiments concerning two classical problems in numerical analysis, numerical integration and function approximation. The novelty of our experiments is that we work on the Grassmannian manifold as an example of a compact Riemannian manifold illustrating theoretical results in the recent literature.
Indeed, recent data analysis methodologies involve kernel based approximation of functions on manifolds and other measure spaces, cf. [24,25] and [16,18,17]. The kernels are build up by what is known as diffusion poly-nomials, which are eigenfunctions of elliptic differential operators, commonly chosen as the Laplace-Beltrami operator when dealing with compact Riemannian manifolds.
Numerical implementations of the approximation schemes require pointwise evaluation of the eigenfunctions. However, explicit formulas for eigenfunctions are only known in few special cases. If the manifold is the unit sphere S d−1 , for instance, then the eigenfunctions of the spherical Laplacian are the spherical harmonics, which are indeed polynomials in the usual sense. The corresponding kernels for the sphere have been computed explicitly in [19,20].
The kernel based approximation requires the computation of the corresponding integral operator, see (7) in Section 3. In the realm of numerical integration, the integral itself is usually approximated by a weighted sum over sample values, see also [14,15]. The latter fits well to the common scenario when the target function needs to be approximated from a finite sample in the first place.
Numerical integration on Euclidean spaces is a classical problem in numerical analysis. Recently, Quasi Monte Carlo (QMC) numerical integration on compact Riemannian manifolds has been studied in [6] from a theoretical point of view, see also [27]. If more and more samples are used, then the smoothness parameter of Bessel potential spaces steers the decay of the integration error. QMC integration has been introduced for the sphere in [7], where many explicit examples are provided and extensive numerical experiments illustrate the theoretical claims.
The major aim of the present paper is to provide numerical experiments for the above integration and approximation schemes when the manifold is the Grassmannian, i.e., the collection of k-dimensional subspaces in R d , naturally identified with the collection G k,d of rank-k orthogonal projectors on R d , cf. [8,Chapter 1].
Therefore, we require explicit formulas of the kernels used in [24,25]. Indeed, the degree of a diffusion polynomial, by definition, relates to the magnitude of the corresponding eigenvalues. We check that diffusion polynomials of degree at most 2t/ √ k are indeed usual multivariate polynomials of degree t restricted to the Grassmannian. The explicit formula for the kernel is derived through generalized Jacobi polynomials. By computing cubature formulas on Grassmannians through some numerical minimization process, we are able to provide numerical experiments for the approximation of functions on Grassmannians and for the QMC integration on Grassmannians supporting the theoretical results in [6,24,25].
The outline is as follows: In Section 2 we recall QMC integration from [6,7], and we recall the approximation scheme from [24,25] in the special case of the Grassmannian in Section 3. Section 4 provides feasible formulations for numerical experiments. Indeed, Section 4.1 is devoted to derive explicit formulas for the involved kernel by means of generalized Jacobi polynomials. We check the relations between diffusion polynomials and ordinary polynomials restricted to the Grassmannian in Section 4.2, and we provide the framework for numerically computing cubatures in Grassmannians in Section 4.3. The numerical experiments are provided in Section 5.

Quasi Monte Carlo integration
We identify the Grassmannian, the collection of k-dimensional subspaces in R d , with the set of orthogonal projectors on R d of rank k, denoted by G k,d := {P ∈ R d×d sym : P 2 = P ; tr(P ) = k}.
Here, R d×d sym is the set of symmetric matrices in R d×d and tr(P ) denotes the trace of P . The dimension of the Grassmannian is dim(G k,d ) = k(d − k). The canonical Riemannian measure on G k,d is denoted by µ k,d , which we assume to be normalized to one. Without loss of generality, we assume k ≤ d 2 throughout since G d−k,d can be identified with G k,d .
As a classical problem in numerical analysis, we aim to approximate the integral over a continuous function f : G k,d → C by a finite sum over weighted samples, i.e., we consider points {P j } n j=1 ⊂ G k,d and nonnegative weights In order to quantify the error by means of the smoothness of f , we shall define Bessel potential spaces on G k,d , for which we need some preparation. Let {ϕ } ∞ =0 be the collection of orthonormal eigenfunctions of the Laplace-Beltrami operator ∆ on G k,d , and {−λ } ∞ =0 are the corresponding eigenvalues arranged, so that 0 = λ 0 ≤ λ 1 ≤ . . .. Without loss of generality, we choose each ϕ to be real-valued, in particular, ϕ 0 ≡ 1. The Fourier transform of Essentially following [6,25], we formally define (I − ∆) s/2 f to be the distribution on G k,d , such that ( i.e., f ∈ H s p (G k,d ) if and only if f ∈ L p (G k,d ) and (I − ∆) s f ∈ L p (G k,d ). Note that this definition is indeed consistent with [6,25], see [6, Theorem 2.1, Definition 2.2] in particular. For s > k(d − k)/p with 1 ≤ p ≤ ∞, the space H s p (G k,d ) is embedded into the space of continuous functions on G k,d , see, for instance, [6]. For 1 < p < ∞, this embedding also follows from results on Bessel potential spaces on general Riemannian manifolds with bounded geometry, cf. [31,Theorem 7.4.5,Section 7.4.2], and on R d with 1 ≤ p ≤ ∞, see [30,Chapter V,6.11]. According to [6], for any sequence of points {P t j } nt j=1 ⊂ G k,d , t = 0, 1, 2, . . ., and positive weights where the constant does not depend on t. Thus, we cannot do any better than the rate n In order to quantify the quality of weighted point sequences {(P t j , ω t j )} nt j=1 , t = 0, 1, 2, . . . , for numerical integration, we make the following definition, whose analogous formulation on the sphere (with constant weights) is due to [7].
holds for all f ∈ H s p (G k,d ). In case p = 2, given s > k(d − k)/2, any sequence of QMC systems {(P t j , ω t j )} nt j=1 for H s 2 (G k,d ) is also a sequence of QMC systems for W s 2 (G k,d ), for all s satisfying s ≥ s > k(d − k)/2, cf. [6].
Especially for the integration of smooth functions, random points lack quality when compared to QMC systems. Proposition 1. For s > k(d − k)/2, suppose P 1 , . . . , P n are random points on G k,d , independently identically distributed according to µ k,d then it holds 2 , so that on average QMC systems indeed perform better than random points for smooth functions. The proof of Proposition 1 is derived by following the lines in [7]. In fact, the result is already contained in [21, Corollary 2.8], see also [26], within a more general setting.
In order to derive QMC systems, we shall have a closer look at cubature points, for which we need the space of diffusion polynomials of degree at most t, defined by see [25] and references therein.
The number t refers to the strength of the cubature.
In the following result, cf. [6, Theorem 2.12], the cubature error is bounded by the cubature strength t, not the number of points.
Weyl's estimates on the spectrum of an elliptic operator yield cf. [22,Theorem 17.5.3]. This implies that any sequence of cubatures of strength t, respectively, must obey n t t k(d−k) asymptotically in t, cf. [10].
cf. [10]. In this case, Theorem 1 leads to so that we have settled that QMC systems do exist, for any s > k(d − k)/p, and can be derived via cubatures.
Remark 1. Cubature points {P j } n j=1 for Π t with constant weights ω j = 1 n are called t-designs. For all t = 1, 2, . . ., there exist t-designs, cf. [29]. The results in [5] imply that there are t-designs satisfying (4) provided that k = 1. However, for 2 ≤ k ≤ d 2 , it is still an open problem if the asymptotics (4) can be achieved by t-designs in place of cubatures.

Approximation by diffusion kernels
One example, where integrals over the Grassmannian are replaced with weighted finite sums, is the approximation of a function f : G k,d → C from finitely many samples. The approximation scheme developed in [24,25] works for manifolds and metric measure spaces in general, but we shall restrict the presentation to the Grassmannian. For where t ≥ 0. It is possible to quantify the best approximation error in dependence of the function's smoothness, see [25,Proposition 5.3] for the following result.
Given f ∈ H s p (G k,d ) we now construct a particular sequence of functions σ t (f ) ∈ Π t , t = 1, 2, . . ., that realizes this best approximation rate. Note that, since the collection {ϕ } ∞ =0 is an orthonormal basis for L 2 (G k,d ), any function f ∈ L 2 (G k,d ) can be expanded as a Fourier series by The approach in [24,25] makes use of a smoothly truncated Fourier expansion of f , where h : R ≥0 → R is an infinitely often differentiable and nonincreasing function with h(x) = 1, for x ≤ 1/2, and h(x) = 0, for x ≥ 1. Using the kernel K t on G k,d × G k,d defined by we arrive after interchanging summation and integration at the following alternative representation Note, that the function σ t (f ) is well-defined for general f ∈ L p (G k,d ), 1 ≤ p ≤ ∞, and it turns out that σ t (f ) approximates f up to a constant as good as the best approximation from Π t , cf. [25,Proposition 5.3].
If f needs to be approximated from a finite sample, then σ t (f ) in (7) cannot be determined directly and is replaced with a weighted finite sum in [25]. Indeed, for sample points {P j } n j=1 ⊂ G k,d and weights {ω j } n j=1 , we define Note that we must now consider functions f in Bessel potential spaces, for which point evaluation makes sense. We shall observe in the following that if samples and weights satisfy some cubature type property, then the approximation rate is still preserved when using σ t (f, {(P j , ω j )} n j=1 ) in place of σ t (f ). However, we need an additional technical assumption on the points {P j } n j=1 , for which we denote the geodesic distance between P, Q ∈ G k,d by where θ 1 , . . . θ k are the principal angles between the associated subspaces of P and Q respectively, i.e., and y 1 , . . . , y k are the k largest eigenvalues of the matrix P Q. We define the ball of radius r centered around P ∈ G k,d by The following approximation from finitely many sample points is due to [25,Proposition 5.3].
Theorem 4. For t = 1, 2, . . ., suppose we are given a sequence of point sets Then the approximation error, Note that the original result stated in [25] requires an additional regularity condition on the samples. This condition is satisfied since we restrict us to positive weights, cf. [15, Theorem 5.5 (a)]. The assumption (9) is a cubature type condition, for which our results in Section 4.2 shall provide further clarification. It indeed turns out that there are sequences {(P t j , ω t j )} nt j=1 satisfying (9) with n t t k(d−k) , in which case (10) becomes Note that the approximation rate in (11) matches the one in (5) for the integration error. The proof of Theorem 4 in [25] is indeed based on Theorem 3 and on the approximation of the integral σ t (f ) in (7) by the weighted finite sum σ t (f, {(P t j , ω t j )} n j=1 ) in (8). For related results on local smoothness and approximation, we refer to [12].

Numerically feasible formulations
This section is dedicated to turn the approximation schemes presented in the previous sections into numerically feasible expressions. In other words, we determine explicit expressions for the kernel K t in (6) and provide an optimization method for the numerical computation of cubature points or QMC systems on the Grassmannian.

Diffusion kernels on Grassmannians
The probability measure µ k,d is invariant under orthogonal conjugation and induced by the Haar (probability) measure µ O(d) on the orthogonal group O(d), i.e., for any Q ∈ G k,d and measurable function f , we have By the orthogonal invariance of the Laplace-Beltrami operator ∆ on G k,d it is convenient for the description of the eigenfunctions to recall the irreducible decomposition of L 2 (G k,d ) with respect to the orthogonal group. Given a nonnegative integer t, a partition of t is an integer vector π = (π 1 , . . . , π t ) with π 1 ≥ . . . ≥ π t ≥ 0 and |π| = t, where |π| := t i=1 π i is the size of π. The length l(π) is the number of nonzero parts of π. The space L 2 (G k,d ) decomposes into where H π (G k,d ) is equivalent to H d 2π , the irreducible representation of O(d) associated to the partition 2π = (2π 1 , . . . , 2π t ), cf. [4,23].
By orthogonal invariance the spaces H π (G k,d ) are eigenspaces of the Laplace-Beltrami operator ∆ on G k,d where, according to [23,Theorem 13.2], the associated eigenvalues are Note, for a given eigenvalue λ the corresponding eigenspace can decompose into more than one irreducible subspace H π (G k,d ).
Note that the following results are translations from representation theory used in [23], see also [4,2], into the terminology of reproducing kernels, where we have only adapted the scaling of the kernels. The space H π (G k,d ) equipped with the L 2 inner product is a finite dimensional reproducing kernel Hilbert space, and its reproducing kernel K π is given by Moreover, K π is zonal, i.e., the value K π (P, Q) only depends on the k largest eigenvalues y 1 (P, Q), . . . , y k (P, Q), of the matrix P Q counted with multiplicities, see [23]. It follows that the kernel K t in (6) is also zonal since it can be written as According to [23], the kernels K π are in one-to-one correspondence with generalized Jacobi polynomials. For parameters α, β ∈ R satisfying 1 2 (m − 1) < α < β − 1 2 (m − 1), the generalized Jacobi polynomials, J α,β π : [0, 1] m → R with l(π) ≤ m, are symmetric polynomials of degree |π| and form a complete orthogonal system with respect to the density where 0 < y 1 , . . . y m < 1, cf. [9,11]. For the special parameters α = k 2 with k ≤ d 2 and β = d 2 , and the normalization J (1, . . . , 1) = dim(H π (G k,d )) the generalized Jacobi polynomials in m = k variables can be identified with the reproducing kernels K π of H π (G k,d ), i.e., Now, (15) and (17) yield that the expression for the kernel K t in (6) can be computed explicitly by (y 1 (P Q), . . . , y k (P Q)).
Thus, avoiding the actual computation of {ϕ } ∞ =0 , we have derived the expression of K t by means of generalized Jacobi polynomials.

Diffusion polynomials on Grassmannians
This section is dedicated to investigate on the relations between diffusion polynomials Π t and multivariate polynomials of degree t restricted to the Grassmannian. Indeed, the space of polynomials on G k,d of degree at most t is defined as restrictions of polynomials by where C[X] t is the collection of multivariate polynomials of degree at most t with d 2 many variables arranged as a matrix X ∈ C d×d . Here, f | G k,d denotes the restriction of f to G k,d . It turns out that Pol t (G k,d ) is a direct sum of eigenspaces of the Laplace-Beltrami operator, i.e., cf. [23, Corollary in Section 11] and also [3, Section 2], which enables us to relate diffusion polynomials to regular polynomials restricted to the Grass-mannian. For k = 1, the eigenvalues (13) directly lead to For general k, the situation is more complicated and needs some preparation.
Lemma 1. Let d, k, t ∈ N with k ≤ d 2 be fixed. Then for any partition π with l(π) ≤ k and |π| ≥ t it holds Note that the right-hand side of (20), up to the square root and the ceiling function, is (13) with π i = t/k, for i = 1, . . . , k.
Proof. In view of (13), let us define For partitions π with |π| ≥ t, we obtain the lower bound (20) by solving the following convex optimization problem Indeed, we shall verify that the minimum is attained at x * := ( t k , . . . , t k ) with by checking the Karush-Kuhn-Tucker (KKT) conditions with µ 0 = 8 t k + 2d − 2(k + 1) and µ i = 2i(k − i), for i = 1, . . . , k. More precisely, denoting the canonical basis in R k by and conclude that the KKT-conditions are satisfied. Hence, (20) holds.
Asymptotically in t, diffusion polynomials of order 2 √ k t are indeed polynomials of degree at most t, and Theorem 3 yields, for f ∈ H s p (G k,d ), For related further studies on dist(f, Pol t (G k,d )), see [28]. In view of Theorem 5, we shall also define cubatures for Pol t (G k,d ).
Definition 3. For {P j } n j=1 ⊂ G k,d and positive weights {ω j } n j=1 , we say that We say that the points {P j } n j=1 ⊂ G k,d are a t-design for Pol t (G k,d ) if (21) holds for constant weights ω 1 = . . . = ω n = 1/n. It turns out that the numerical construction of cubature points and t-designs for Pol t (G k,d ) is somewhat easier than for Π t directly, which is outlined in the subsequent section.
Remark 2. Since Pol t (G k,d ) are restrictions of ordinary polynomials, we observe Pol t1 (G k,d ) · Pol t2 (G k,d ) ⊂ Pol t1+t2 (G k,d ). Thus, {(P t j , ω t j )} nt j=1 being cubatures for Pol 2t (G k,d ) yields that (9) is satisfied when Π t is replaced with Π s(t+1)− . The latter implies that we must then also replace For general k, the second set inclusion in Theorem 5 is sharp because λ(t, 0, . . . , 0) = 4t 2 + 2t(d − 2). The first set inclusion in Theorem 5 may only be optimal for t being a multiple of k. To prepare for our numerical experiments later, we shall investigate on G 2,d more closely.
Theorem 6. For k = 2, we obtain Note that s(t) in Theorem 6 satisfies s 2 (t) = λ( t 2 , t 2 ). It matches the definition in Theorem 5 provided that t is even. For odd t, s(t) in Theorem 6 is indeed larger than in Theorem 5, and the difference of the squares is 4.

Worst case error of integration on Grassmannians
Given some subspace H of continuous functions on G k,d the worst case error of integration (with respect to some norm · on H) for points {P j } n j=1 ⊂ G k,d and weights {ω j } n j=1 is defined by see also [21,26]. If H = H K is a reproducing kernel Hilbert space, whose reproducing kernel K is with r π ≥ 0, |π| ≥ 0, and sufficient decay of the coefficients, then the associated inner product is and the Riesz representation theorem yields Note that the worst case error is a weighted 2 -average of the integration errors of the basis functions ϕ 0 , ϕ 1 , ϕ 2 , . . .. Recall, for instance, H s 2 (G k,d ) is a Hilbert space with inner product and the Bessel kernel on the Grassmannian is K s (y 1 (P, Q), . . . , y k (P, Q)).
If s > k(d − k)/2, then it is easily checked that K s B is the reproducing kernel for H s 2 (G k,d ) with respect to the inner product (24), see also [6]. Note that the polynomial space Pol t (G k,d ) is also a reproducing kernel Hilbert space. Indeed, given a partition π with |π| ≤ t and l(π) ≤ k, the reproducing kernel of H π (G k,d ) with respect to the L 2 inner product is K π in (14). Due to (19), the reproducing kernels for Pol t (G k,d ) are exactly with r π > 0, |π| ≥ 0. Note that R t is indeed reproducing as a finite linear combination with nonnegative coefficients of reproducing kernels, and it reproduces Pol t (G k,d ) because of (19) and none of the coefficients vanish. Now, by Definition 3 any cubature for Pol t (G k,d ) has zero worst case error independent of the chosen norm, and thus independent of R t . A particularly simple reproducing kernel for Pol t (G k,d ) is see, for instance, [13]. Hence, formula (23) provides us with a simple method to numerically compute cubature points by some minimization method. In particular, t-designs {(P j , 1 n )} n j=1 are constructed by minimizing and checking for equality, which implies wce Polt, · R t ({(P j , ω j )} n j=1 ) = 0.

Numerical experiments
We now aim to illustrate theoretical results of the previous sections. The projective space G 1,d can be dealt with approaches for the sphere by identifying x and −x. The space G d−1,d can be identified with G 1,d , so that the first really new example to be considered here is G 2,4 . We computed points {P t j } nt j=1 ⊂ G 2,4 , for t = 1, . . . , 14, with worst case error wce Polt, · L 2 ({(P t j , 1/n t )} nt j=1 ) < 10 −7 by a nonlinear conjugate gradient method on manifolds, cf. [21, Section 3.3.1], see also [1,Section 8.3]. Although the worst case error may not be zero exactly, we shall refer to {P t j } nt j=1 in the following simply as t-designs. Note that G 2,4 has dimension dim(G 2,4 ) = 4, so that the number of cubature points must satisfy n t t 4 . Indeed we chose We emphasize that for t = 14 we computed n 14 = 8.475 projection matrices which almost perfectly integrate 25.425 polynomial basis functions.

Integration
In what follows we consider two positive definite kernels K 1 (P, Q) = (2 − tr(P Q)) 3 + 2 tr(P Q), It can be checked by comparison to the Bessel kernel, cf. [6], that the reproducing kernel Hilbert space H K1 equals the Bessel potential space H 7 2 2 (G 2,4 ), i.e., the corresponding norms are comparable. In contrast, the reproducing kernel Hilbert space H K2 is contained in the Bessel potential space H s 2 (G 2,4 ) for any s > 2. The worst case errors can be computed by dt is the hyperbolic sine integral. In view of illustrating Proposition 1, note that random P ∈ G k,d distributed according to µ k,d can be derived by P := Z(Z Z) −1 Z , where Z ∈ R d×k with entries that are independently and identically standard normally distributed, cf. [8, Theorem 2.2.2]. Figure 1 depicts clearly the superior integration quality of the computed cubature points over random sampling. Moreover, it can be seen that the theoretical results in Proposition 1 and Theorem 1 with (5) are in perfect accordance with the numerical experiment, i.e., the integration errors of the random points scatter around the expected integration error and cubature points achieve the optimal rate of n − 7 8 for functions in H 7 2 2 (G 2,4 ). In Figure 2 we aim to show the contrast between the integration of functions in H K1 and H K2 by using the computed t-designs. We know by Theorem 1 that the sequence of t-designs with a number of cubature points n t t 4 is a QMC system for any s > 2. Since H K2 is contained in any Bessel potential space H s 2 (G 2,4 ), for s > 2, we expect a super linear behavior in our logarithmic plots. Indeed, Figure 2 confirms our expectations. For t ≥ 11, the effect of the accuracy of the t-designs used becomes significant for integration of smooth functions. For that reason, we added the dashed red line, which represents the accuracy 10 −7 of the computed t-designs.

Approximation
Similar, as in the previous section we aim to approximate a smooth and a nonsmooth function, namely f 1 (P ) = K 1 (I 2 , P ), f 2 (P ) = K 2 (I 2 , P ), P ∈ G 2,4 , where K 1 , K 2 are from the previous section and I 2 is a projection matrix with 2 ones on the upper left diagonal. This time we observe that the function f 1 is contained in H 3 ∞ (G 2,4 ) but f ∈ H 3+ ∞ (G 2,4 ), for all > 0. For the smooth function f 2 , we have f 2 ∈ H s ∞ (G 2,4 ), for any s > 0. Since the computed t-designs are with respect to Pol t (G 2,4 ) and not Π t (G 2,4 ), we need an additional scaling of s(t) in σ t . According to Example 1, the choice s(t) = 2(t 2 + 3t + 3 + (−1) t ) − t √ 2, for small > 0, yields Π s(t) · Π s(t) ⊂ Pol 2t (G 2,4 ). For numerical experiments, we take to be smaller than the machine precision, so that it is effectively zero. Hence, in accordance with Theorem 4, we use the following kernel based approximation h(s(t) −2 λ(π))K π (P j , ·), where X 2t = {(P 2t j , 1/n 2t )} n2t j=1 and 1, x ≤ 1/2, 0, otherwise. The approximation error is determined by randomly sampling altogether 50000 points. The first 25000 are pseudo random according to µ 2,4 . Since f 1 has a nonsmooth point at I 2 the maximal error is expected around this point. Therefore, we sampled the other 25000 from normally distributed points around that point I 2 with variance 0.15 and 0.5 in the matrix entries, i.e., we choose Z ∈ R 4×4 with independent and identically distributed entries according to a normal distribution with mean zero and variance 0.15 and 0.5, respectively, and then project I 2 + Z onto G 2,4 , which we accomplished by a QR-decomposition in Matlab.
In Figure 3, we can observe the predicted decay in Theorem 4 for the function f 1 ∈ H 3 ∞ (G 2,4 ). Furthermore, as expected for the smooth function f 2 , the error appears to decrease super linearly.