Thinplate Splines on the Sphere

In this paper we give explicit closed forms for the semi-reproducing kernels associated with thinplate spline interpolation on the sphere. Polyharmonic or thinplate splines for ${\mathbb R}^d$ were introduced by Duchon and have become a widely used tool in myriad applications. The analogues for ${\mathbb S}^{d-1}$ are the thin plate splines for the sphere. The topic was first discussed by Wahba in the early 1980's, for the ${\mathbb S}^2$ case. Wahba presented the associated semi-reproducing kernels as infinite series. These semi-reproducing kernels play a central role in expressions for the solution of the associated spline interpolation and smoothing problems. The main aims of the current paper are to give a recurrence for the semi-reproducing kernels, and also to use the recurrence to obtain explicit closed form expressions for many of these kernels. The closed form expressions will in many cases be significantly faster to evaluate than the series expansions. This will enhance the practicality of using these thinplate splines for the sphere in computations.


Introduction
In this paper we give explicit closed forms for the semi-reproducing kernels associated with thinplate spline interpolation on the sphere. Polyharmonic or thinplate splines for R d were introduced by Duchon in his classic papers [10,11] and have become a widely used tool in myriad applications. The analogues for S d−1 ⊂ R d are the thinplate splines for the sphere. The topic was first discussed by Wahba [22,23] in the early 1980's, for the S 2 case. Wahba presented the associated semi-reproducing kernels as infinite series. These semi-reproducing kernels play a central role in expressions for the solution of the associated spline interpolation and smoothing problems.
The main aims of the current paper are to give a recurrence for these semi-reproducing kernels, and also to use the recurrence to obtain explicit closed form expressions. Here we are building on previous work of Martinez-Morales [16]. Unfortunately, there are errors in the theory presented in [16] and consequently many of the expressions given there for the kernels are incorrect. The closed form expressions given here will usually be significantly faster to evaluate than the series expansions. This will enhance the practicality of using the thinplate splines for the sphere in computations. arXiv:1801.01313v2 [math.CA] 12 Aug 2018 The paper is laid out as follows. Section 2 discusses the central role played by semi-reproducing kernels in the solution of both interpolation and penalized least squares fitting problems. Section 3 develops semi-reproducing kernels associated with the thinplate splines on the sphere, that is, semi-reproducing kernels associated with minimum energy interpolation and penalized least squares fitting problems with a particular choice of energy. The energy chosen being that naturally associated with iterated Laplace-Beltrami operators. These semi-reproducing kernels are given in this section as infinite series. Section 3 also recalls known results concerning Fourier-Gegenbauer expansions that will be needed later. Section 4 motivates the construction of an operator T and its adjoint T * . It also presents some fundamental properties of these operators. These operators were initially developed in Martinez-Morales [16]. Section 5 gives a recurrence for the various thinplate spline kernels K d,m (x, y), where d indicates the dimension and m is the power of the associated differential operator. More precisely, it gives a recurrence for the related functions k d,m : [−1, 1] → R where k d,m ( x, y ) = K d,m (x, y). Sections 6 and 7 give short closed form expressions for many of the functions k d,m .
For ease of access to the relevant background we will base our notation on that used in Dai and Xu [8]. Occasionally, when the value for the dimension d is particularly important, we will supplement the symbol they use with a d.

Reproducing kernels and approximation
Let us begin with summarising the main ideas of reproducing kernels in indefinite inner product spaces and their role in both interpolation and penalized least squares fitting problems. These results show the central role of reproducing kernels in the solution of these problems. We focus on the specific case of reproducing kernels 1 for semi-Hilbert spaces, or simply semi-reproducing kernels, as this is the appropriate framework for thinplate spline approximation. For further definitions and basic properties of semi-reproducing kernels we refer to [3,4,7,20]. Our treatment of the relevant interpolation and penalized least squares problems is based on that of Strauss [21].
Let D be a subset of R d . Consider approximation from a vector space over the reals H ⊂ C(D). Assume the space H is endowed with a semi-inner product (·, ·), that is, the inner product is lacking definiteness. Thus, there are non-zero vectors f ∈ H with (f, f ) = 0. Further assume that the kernel H 0 of the semi-inner product (·, ·) is finite-dimensional, i.e., dim H 0 = m < ∞, and (f, f ) = 0 if and only if f ∈ H 0 .
A standard approach to deal with semi-inner products is to supplement the semi-inner product with an inner product on H 0 thereby obtaining a definite inner product on the space H (see [5]). Towards this aim, we need to decompose the space H into a direct sum H 0 ⊕ H 1 , such that the given semi-inner product provides a definite inner product on the subspace H 1 . Given m linearly independent functionals spanning the dual of H 0 , we define H 1 as the space of functions in H which are mapped onto zero by all these functionals.
Let us recall this approach using point evaluations. Nevertheless, it is important to note that there are many choices for such a set of functionals. Clearly, the decomposition obtained for the space H depends upon the choice made for the m functionals.
Definition 2.1. A set of distinct points X is said to be unisolvent for H 0 if the only function in H 0 which is zero at all points of X is the zero function.
Given a unisolvent set X = {z 1 , . . . , x m } for H 0 , where m = dim(H 0 ), the set of point evaluation functionals {δ x : x ∈ X } is linearly independent on H 0 . Hence, we can find a Lagrange basis u 1 , . . . , u m of H 0 with respect to X , i.e., u i (x j ) = δ ij , 1 ≤ i, j ≤ m. Then defines an inner product on H 0 . The basis u 1 , . . . , u m is orthogonal with respect to the inner product [·, ·] 0 . Furthermore, the mapping The definition of inner product for the full space H which follows will make P 0 the orthogonal projection onto H 0 . The subspace H 1 can then be defined via the projector P 1 = I − P 0 , i.e., Since H 0 is the kernel of the semi-inner product (·, ·), i.e., (f, f ) = 0 iff f ∈ H 0 , the semi-inner product is definite on H 1 via construction. Therefore, defines a definite inner product on H (see [5] for further details). We are interested in Hilbert spaces carrying the special property of being reproducing kernel spaces. There is a one-to-one correspondence between reproducing kernel Hilbert spaces and positive definite kernels. A similar relation holds true for semi-reproducing kernel Hilbert spaces.
Definition 2.2. Given n ∈ N, a pair (X , a) with X = {x 1 , . . . , x n } a set of distinct points from D and a = (a 1 , . . . , a n ) T ∈ R n is called an The set of all H 0 -increments is denoted by H ⊥ 0 .
Note that an H 0 -increment can naturally be identified with a linear functional vanishing on H 0 . The structure of semi-reproducing kernels already shows that we can expect the reproducing kernel for a semi-Hilbert space to be positive definite only on a suitable subspace.
for all H 0 -increments (X , a). K is said to be strictly conditionally positive definite with respect to H 0 if the inequality in (2.1) is strict whenever in addition a is nonzero.
There is a correspondence between conditionally positive definite kernels and semi-reproducing spaces (see [2,17,20]). The associated reproducing property can again be stated in terms of H 0 -increments. Definition 2.4. Let H ⊂ C(D) be a semi-Hilbert space, i.e., H is a semi-inner product space with semi-inner product (·, ·) the kernel H 0 of which is finite-dimensional, and H is complete with respect to the induced semi-norm. A symmetric kernel K : D × D → R is called a semireproducing kernel for H if (·, ·) reproduces H 0 -increments, i.e., for all λ X ,a annihilating H 0 the following two properties hold n j=1 a j K(·, x j ) ∈ H, Given m = dim H 0 linearly independent functionals λ 1 , . . . , λ m on H 0 and the corresponding Lagrange basis obviously provides a reproducing kernel for H 0 with respect to the inner product Using the orthogonal projection we can again define H 1 = (I − P 0 )H. By construction, Furthermore, if K is a semi-reproducing kernel of H the kernel is the reproducing kernel of H 1 . Here, the superindex in the last term indicates the functional operating on the first and second variable, respectively, Thus, the space H is a reproducing kernel Hilbert space itself with reproducing kernel See [7,20] for details. Note that if K is given such that P 0 K(·, x) = 0 for all x ∈ D, then the projection in the above expression vanishes, i.e. K 1 = K. The framework of reproducing kernel spaces allows us to consider regularized interpolation problems in broad mathematical generalities (see [3,Chapter 2.1]). The rather beautiful result of Strauss [21] concerning mixed interpolation and regularized least squares problems provides a special case of [3, Theorem 59].
The following notation will be used. For a set of points X = {x 1 , . . . , x n } ⊂ D and a function f on D we write f X for the vector (f (x 1 ), . . . , f (x n )) T , K X for the n × n matrix with ij-entry K(x i , x j ), and C X for the m × n matrix u i (x j ) , where u 1 , . . . , u m is a basis of H 0 . Furthermore, W † denotes the pseudo-inverse of the matrix W appearing in the assumptions of the following theorem. Under the assumptions of the theorem Theorem 2.5. Let H and (·, ·) be as in Definition 2.4, and let K be a semi-reproducing kernel for (H, (·, ·)) with respect to H 0 . Further suppose that K is strictly conditionally positive definite with respect to H 0 . Let µ > 0, n ≥ m = dim H 0 , 0 ≤ p ≤ n, and W be an n × n matrix of the form where R is p × p and symmetric positive definite. Then given any set X = {x 1 , . . . , x n } of n distinct points in D which is unisolvent for H 0 , and n corresponding values y i ∈ R, there is a unique member of the space H minimizing the quadratic functional over those functions in H which satisfy the interpolation conditions This function can be written in the form where the coefficients a = (a 1 , . . . , a n ) T and b = (b 1 , . . . , b m ) T are the solution of the system Note that the statement reduces to the well-known result concerning the solution of the smoothest interpolation problem when p = 0, and to a known expression for the smoothing spline when p = n.

Series representations for thinplate spline kernels on the sphere
For functions on R d interpolating and smoothing with thinplate/polyharmonic splines associated with the energy are very successful approximation methods. For sufficiently smooth functions f , decaying sufficiently fast at infinity, integration by parts gives where is the Laplacian. Analogues for the sphere come from considering instead of the Laplacian the Laplace-Beltrami operator , and working on the "Fourier" side since the spherical harmonics {Y nj } are both a complete orthonormal system for L 2 S d−1 and also eigenfunctions of the Laplace-Beltrami operator. More precisely, For the explicit value of the dimension N d,n see (3.1), below. This section will consider corresponding spaces of functions on the sphere and the relevant semi-reproducing kernels K d,m . The central role played by these semi-reproducing kernels in minimal energy interpolation and regularised least squares fitting is clear from the discussion in Section 2, and in particular Theorem 2.5. This topic was first considered by Wahba [22] in the S 2 case. See also Wahba's monograph [24]. The reader can find valuable additional relevant material in Freeden, Gervens and Schreiner [12, Chapter 5], Levesley, Light, Ragozin and Sun [15], and Cheney and Light [7,Chapter 32]. The material in these references differs somewhat from what appears here. Usually this is due to a treatment based on reproducing kernels rather than semi-reproducing kernels, or to a different choice of the energy. Gneiting [13] gives an excellent survey of recent work concerning kernels for the sphere.
Let H d n denote the space of real harmonic polynomials homogeneous of degree n on R d . The spherical harmonics are the restrictions of these to the sphere S d−1 . In a slight abuse of notation the space of spherical harmonics of degree n on S d−1 is also written H d n . The dimension of the space is Spherical harmonics are a complete orthogonal system on L 2 S d−1 with respect to the inner product For the set {Y n j : 1 ≤ j ≤ N d,n }, being an orthonormal basis for H d n , the addition formula shows that the reproducing kernels of the spaces H d n are the zonal polynomials W λ n which are Gegenbauer polynomials normalized so that W λ n (1) = 1. The Gegenbauer polynomial of order λ ≥ 0 and degree n ∈ N 0 is defined as the hypergeometric polynomial The Gegenbauer polynomials are orthogonal with respect to the inner product Note that since the convolution of zonal functions on the sphere remains zonal, the inner product (3.2) reduces to (3.4) for zonal functions, where λ = d−2 2 . The remaining weight function has the integral Although, Gegenbauer polynomials provide a complete orthogonal system for all λ > − 1 2 , we will fix λ = d−2 2 throughout this paper. The constant C λ n (1) relates to the dimension N d,n given in (3.1); indeed, Since {Y n j : 1 ≤ j ≤ N d,n , n ∈ N 0 } is a complete orthonormal system for L 2 S d−1 , we can consider Fourier series where Further, F d, m is the space of all functions f ∈ F d m with Fourier coefficients a nj = 0, for all 1 ≤ j ≤ N d,n and 0 ≤ n ≤ . In what follows we will consider approximations s to f whose smoothness is measured by an inner product on F d, m with an additional spherical polynomial part of degree viewed as a trend. The case most frequently occurring in the literature is that of F d,0 m .
b nj Y n j and m even, by the extended Parseval identity in the space L 2 S d−1 . In view of the equality between expressions (3.6) and (3.7) define an "energy" semi-inner product for F d m by It is clear from (3.6) and (3.7) that this is an analogue of the usual semi-inner product associated with smoothing splines on R d .
by the addition formula (3.3). K d,m, is clearly a zonal kernel since it depends only on the cosine of the angle between x and y. We will use the notation K d,m for the kernels K d,m,0 , and call these kernels the thinplate spline kernels for the sphere S d−1 . The K 2,m and K 3,m kernels were initially introduced by Wahba [22]. Since K d,m, (x, y) is zonal we can sensibly define functions k d,m, by Explicitly, We will refer to the functions k d,m, as semi-reproducing functions, and the functions k d,m = k d,m,0 as thinplate spline functions. To show property (i) let y be some fixed point in S d−1 and define k y (·) = K(·, y). From the definition of K, k y has Fourier coefficients a nj = [n(n + d − 2)] −m Y n j (y), n ≥ + 1 and 1 ≤ j ≤ N d,n . (3.12) Hence, where in the second to last step the addition formula ( That is the reproducing property holds. Proof of part (b). In view of the characterisations of strict positive definiteness of zonal kernels given by Chen, Menegatto and Sun [6] for d > 2, and by Menegatto [18] for d = 2, part (b) follows from the signs of the Gegenbauer coefficients of K d,m displayed in equation (3.9).
Above we have shown K d,m, is the reproducing kernel for the space F d, m which arises from using a Fourier projection onto spherical polynomials to split the space F d m into a direct sum H 0 ⊕ F d, m , where H 0 = ∪ n=0 H d n is the space of spherical polynomials of degree at most . For interpolation problems it is natural to use instead a projection onto polynomials interpolating at a certain finite set of points, and this results in a different direct sum decomposition. Fortunately, here the semi-reproducing kernel approach becomes especially convenient as K d,m, is a semi-reproducing kernel for the space F d, m with respect to the space of polynomials H 0 . This is the content of the following easily shown lemma whose proof is included for the sake of completeness. Proof . F d m = H is considered as a semi-inner product space with semi-inner product (·, ·) m, . Choose P 0 as Fourier projection onto the kernel H 0 of the semi-inner product, which is the space of spherical polynomials of degree not exceeding . Set the first property, i.e., property (2.2), of a semi-reproducing kernel. Also, given any f ∈ K d m the direct sum splitting allows us to write which follows from (h 0 , h 1 ) m, = 0 for all h 0 ∈ H 0 and h 1 ∈ F d, m , and also from K d,m, being the reproducing kernel for F d, m . Continuing, using the vanishing property of H 0 -increments, the second property, i.e., property (2.3), of a semi-reproducing kernel. Therefore, K d,m, is a semi-reproducing kernel for F d m with semi-inner product (·, ·) m, , as required.
In view of Lemma 3.2, Theorem 2.5 concerning the solution of interpolation and penalized least squares fitting problems applies to interpolation and smoothing problems on the sphere. Applying the theorem the semi-reproducing kernel K = K d,m, , defined in equation (3.9), plays a central role in interpolation and penalized least squares fitting problems posed in the space H = F d m , with semi-inner product (·, ·) m, with respect to the finite dimensional subspace H 0 = ∪ n=0 H d n . In this section we have given series expansions for these kernels. In Section 5, Theorem 5.1, below, we will provide a recurrence relation for the particularly important thinplate spline kernels, K d,m, , via a recurrence for the corresponding functions k d,m . In Sections 6 and 7 the recurrence relation will be used to give short explicit expressions for many of the thinplate spline kernels.

The operator T and its adjoint T *
In this section we discuss an operator T , and its adjoint T * , which will be crucial parts of the recurrence for the thinplate spline functions k d,m . These operators were defined by Martinez-Morales in [16].
In view of the series expansions for the kernels k d,m given in equation (3.11) a multiplier operator with Fourier multiplier of (n(n + 2λ)) −1 would transform k d,m into k d,m+1 . The operators T and T * , are discussed below have some, but not quite all, the desired properties.
Note that the differential equation for the Gegenbauer polynomials is given by [1, equation (22.6.5)] or [19, Call the operator on the left hand side of the equation D λ . Then, for λ > 0 and n ∈ N, Let us formally invert the equation D λ f = g using averages centered at 1. Dealing with the outer derivative in D λ we obtain that x .
If f (y) is continuous at y = 1, (4.1) the term on the right hand side for y = 1 vanishes. We can then proceed obtaining Setting f (x) = C λ n (x) and g(x) = −n(n + 2λ)C λ n (x), condition (4.1) is clearly satisfied. We therefore obtain the following statement.
Proof . We would like to give a direct proof of the statement. Towards this goal, we use an integral given in [ Decomposing the integral over [0, 1] into two integrals over [0, x] and [x, 1], respectively, we can use the formula to obtain that Note that the function on the right hand side vanishes at y = 1. Towards the claim, it remains to integrate the polynomial C λ+1 n−1 which readily follows from [19, equation (18.9.19)] completing the proof.
Similarly, we could have treated the average centered at −1 giving the following result.
Proposition 4.2. For λ > 0 and n ∈ N we have that The operators have been defined in [16], showing that T * λ is the adjoint of T λ with respect to the inner product (3.4). To be precise, the following statement holds true (cf. [16,Theorem 3]). and The proof exploits the fact that with the weight 1 − z 2 λ− 1 2 in the inner integral, continuity of f suffices to cope with the singularity introduced by the weight 1 − y 2 −λ− 1 2 in the outer integral of T λ . Based on this observation, the theorem follows from Fubini's theorem.
Propositions 4.1 and 4.2 above thus show that the polynomials C λ n are basically -up to a constant -eigenfunctions of T * and T , respectively. This is somewhat obvious from the fact that both T and T * invert D λ . Both propositions have been derived in [16, Lemma 1] via a different proof employing the Rodriguez formula of the Gegenbauer polynomials.

A recurrence for the thinplate spline functions for the sphere
This section concerns a recurrence for the thinplate spline functions k d,m for S d−1 . The recurrence will be used in Sections 6 and 7 to give short explicit forms for many of the functions k d,m . It is important to note that the corresponding recurrence given in [16,Theorem 4] is incorrect and does not yield the thinplate spline functions k d,m . Consequently, many of the explicit formulas claimed for the functions k d,m in the paper [16] are also incorrect. 1], and e 0 (x) = 1 for all x ∈ [−1, 1]. The thinplate spline functions, k d,m , m ∈ N, for the sphere S d−1 , defined via the series (3.11), are alternatively generated by the recurrence Proof . Fix d ≥ 2 and define a sequence of functions (f 1 , f 2 , . . .) by the recurrence g n C λ n in terms of the Gegenbauer polynomials C λ n . The uniqueness theorem tells us that functions with the same coefficients are identical.
First, observe that the constant term in the definition of f m ensures that the zeroth Fourier coefficient (f m ) 0 is zero. Therefore, consider in what follows Fourier coefficients of index n ≥ 1.
Clearly, e 0 ∈ which by Proposition 4.1 yields In the special case of f = e 0 , the first integral vanishes due to orthogonality (3.5). Furthermore, We therefore obtain that for n ≥ 1, (T λ e 0 ) n = − 1 n(n + 2λ) .
Thus, the function f 1 generated as specified in equation (5.1) has the same Fourier coefficients as the thinplate spline function k d,1 of equation (3.11). Thus, by uniqueness, it is k d,1 . Now returning to general functions f we can rewrite (5.2) as , when n ≥ 1.
Hence f m has the same Fourier coefficients as k d,m . Therefore, by uniqueness, it is k d,m . That is the thinplate spline functions, k d,m , are generated by the recursion of the theorem.
6 Explicit forms for some of the thinplate spline functions k d,m The reader will recall the correspondence between the zonal thinplate spline kernels K d,m and the associated functions k d,m , see (3.10). The recurrences of the previous section yield explicit formulas for many of the thinplate spline functions k d,m . A sample of these explicit expressions is presented below, thereby correcting the expressions given in [16].
6.1 Functions for S 1 See Wahba [24, p. 22] for explicit forms of these functions in terms of Bernoulli polynomials.
7 Explicit formulas for the thinplate spline functions k d,1 , d > 2 In the section simple explicit formulas will be obtained for the thinplate spline functions k d, 1 .
As explained at the end of Section 3 the function k d,1 is associated with approximation problems on S d−1 . Theorem 5.1 and the formula for the operator T λ given in Proposition 4.2 lead to a method of calculating k d,1 . Define Then, with λ = (d − 2)/2, Actually it is somewhat easier to deal with the indefinite integral where here we mean any fixed representative value of the indefinite integral, since the second term in (7.2) deletes the constant term. Then, from equation (7.2),

The functions k d,1 when d is even
This subsection considers the functions k d,1 when d is even. Thus in this subsection λ is a positive integer.
Analogously, define the indefinite integral H β (y) = 1 − y 2 −(2β+1)/2 dy, β > −1/2. Families of representatives may be generated by the recurrence Explicitly, for λ ∈ N, starting from H 1 , as given by equation (7.9), the recurrence generates representatives of the form Proof . Equations (7.9) and (7.10) follow immediately from the definition (7.8). The recurrence (7.11) follows from the definition (7.8) via an easy integration by parts. The general form of the explicit expression for H λ (y), λ ∈ N, given in equation (7.12), is clear from the expression for H 1 (y) and the recurrence (7.11). Consider now the expression for the coefficients b λ j ocurring in equation (7.12). From the recurrence and the formula for H 1 the term involving y 1 − y 2 −(2j−1)/2 first appears for λ = j, where it has the value b j j = 1/(2j − 1). This term is then propagated to the functions H λ , with λ > j, via the recurrence. Hence, Given the formula (7.7) for G λ (y) when λ is a positive integer, the definition (7.8) of H λ and the definition (7.3) of S λ , S λ (y) = 1 − y 2 −(2λ+1)/2 a λ 0 π 2 + arcsin(y) + a λ 1 y 1 − y 2 1/2 + · · · + a λ λ y 1 − y 2 (2λ−1)/2 dy = a λ 0 π 2 + arcsin(y) dH λ (y) where I 1 and I 2 are the first and second indefinite integrals, respectively. Ignoring the constant parts in the indefinite integrals a representative value of I 2 is A representative value of I 1 is Now note that Hence, the terms in I 1 and I 2 involving ln 1 − y 2 have coefficients of equal magnitude and opposite sign. Thus, we conclude that a representative value of S λ (y) is Now recall from equation (7.4) that where λ = (d − 2)/2. To calculate this quantity first define f µ as the Beta integral Then Therefore, substituting the various quantities into equation (7.4)

The functions k d,1 when d > 1 is odd
Now turn to the calculation of the functions k d,1 when d > 1 is odd. Then the Gegenbauer parameter λ = (d − 2)/2 = κ + 1 2 , for some nonnegative integer κ. We will need the following technical lemmas. Lemma 7.3. Let α be an integer and β a nonnegative integer. If α is nonnegative further suppose β < a. Then Proof . The identity will be proven by induction on β. For β = 0 the result is immediate. Now assume that the identity is true for β = κ − 1 where κ ∈ N. Then J κ,a = κ 0 Applying the induction hypothesis twice showing that the identity also holds for β = κ. Substituting into equation (7.13) yields the result.