Intrinsic Random Functions on the sphere
Chunfeng Huang1 , Haimeng Zhang2 , Scott M. Robeson3 , Jacob Shields4
Abstract. Spatial stochastic processes that are modeled over the entire Earth’s surface require
statistical approaches that directly consider the spherical domain. Here, we extend the notion

arXiv:1606.01950v1 [math.ST] 6 Jun 2016

of intrinsic random functions (IRF) to model non-stationary processes on the sphere and show
that low-frequency truncation plays an essential role. Then, the universal kriging formula on the
sphere is derived. We show that all of these developments can be presented through the theory
of reproducing kernel Hilbert space. In addition, the link between universal kriging and splines is
carefully investigated, whereby we show that thin-plate splines are non-applicable for surface fitting
on the sphere.
Keywords. Kriging, Reproducing kernel Hilbert space, Splines, Stationary

1

Introduction

Global-scale phenomena, of which there are a multitude of important applications ranging from
climate science to epidemiology, can be viewed as random processes on the sphere. A common
assumption when modeling such processes is second-order stationarity (or stationarity for short in
this paper), where the mean is constant and the covariance function at two locations is assumed to
depend only on their distance (or, in the case of circular or spherical phenomena, angular separation;
see [26, 11, 9]). This assumption is difficult to evaluate and often deemed unrealistic in practice.
Several approaches have been proposed to relax this assumption such as axial symmetry [15, 32, 12]
and kernel convolution [42, 10]. Based on generalized random functions in [8], Matheron [20]
introduced a flexible family of non-stationary processes, named intrinsic random functions (IRFs),
where the process is assumed to have lower monomials as its mean and the transformed process
becomes stationary. The influence of this approach is profound, as it provides a basis for kriging
[5, 33, 2] in practice. Most developments and uses of IRFs have been in Euclidean spaces, where
1
2

Department of Statistics, Indiana University, Bloomington, IN. Email: huang48@indiana.edu.
Department of Mathematics and Statistics, University of North Carolina at Greensboro, Greensboro, NC.

Email:haimengzhanguncg@gmail.com
3
Department of Geography, Indiana University, Bloomington, IN. Email:srobeson@indiana.edu
4
Department of Statistics, Indiana University, Bloomington, IN. Email:shields@imail.iu.edu

1

the differencing or differential operation is used for the transformation. For example, IRFs on the
real line, also known as processes with random stationary increments, were discussed in [39] and
[41]. IRFs on other spaces are less explored. In Huang et al. [13], IRFs on the circle are developed,
where we showed that low-frequency truncation replaces differencing to achieve stationarity.
In this paper, we first define IRFs on the sphere. While lower monomials are used in Euclidean
spaces, the counterparts on the sphere are the lower order spherical harmonics. We formally show
that an IRF on the sphere is characterized by its frequency-truncated process (Theorem 1 in Section
2.1). Then, we derive the universal kriging formula for IRFs on the sphere, where we show that
the coefficients for lower spherical harmonics do not need to be estimated. In addition, we demonstrate that IRFs can be viewed through the theory of reproducing kernel Hilbert space (RKHS,
[1, 37, 34, 14, 13]). Based on this, kriging on the sphere is shown to be the same as the smoothing
formula in RKHS. We formally establish the equivalence between splines and kriging. Finally, we
carefully investigate splines on the sphere and find that the popularly used thin-plate spline approach is not applicable for surface fitting on the sphere.

2

Main results

Let Z(x) be a real-valued random process on a unit sphere S 2 , where x = (ψ, ζ) with the longitude
ψ ∈ [0, 2π) and the latitude ζ ∈ [0, π]. Assume Z(x) is continuous in quadratic mean and its
covariance function cov(Z(x), Z(y)) is strictly positive definite. One can expand Z(x) in spherical
harmonics which are convergent in quadratic mean [40, 15, 28, 29]:
Z(x) =

∞ X
l
X

Zl,m Ylm (x),

l=0 m=−l

where the real-valued spherical harmonic functions are

q
2l+1 (l−m)! m
m (x) =

P (cos ζ) cos(mψ), m = 1, . . . , l,
Y

l


q 2π (l+m)! l
0
Yl0 (x) = 2l+1
4π Pl (cos ζ),

q


 Y −m (x) = 2l+1 (l−m)! P m (cos ζ) sin(mψ), m = 1, . . . , l,
l
2π (l+m)! l

2

and Plm (·) are the associated Legendre polynomials and Pl0 (·) ≡ Pl (·) are the Legendre polynomials.
The coefficients Zl,m are random variables:
Z 2π Z π
Zl,m =
Z(ψ, ζ)Ylm (ψ, ζ) sin ζdζdψ.
0

0

We denote the group of rotations of the sphere by G. When the process Z(x) is assumed to be
stationary, that is, the mean is constant E(Z(x)) = µ and, for any two points x, y ∈ S 2 ,
cov(Z(x), Z(y)) = cov(Z(gx), Z(gy)),
for all g ∈ G. The stationarity implies that the coefficients Zl,m are uncorrelated [26, 40, 28]:
cov(Zl,m , Zl′ ,m′ ) = al I(l, l′ )I(m, m′ ),
where al ≥ 0,

P

l (2l

+ 1)al < ∞ and I(l, l′ ) is the indicator function taking value 1 if l = l′ , and

zero otherwise. By the addition theorem of Legendre polynomials, the covariance function
cov(Z(x), Z(y)) =

∞
X
l=0

al

l
X

Ylm (x)Ylm (y) =

∞
X
2l + 1
l=0

m=−l

4π

al Pl (cos(d(x, y))),

(1)

where d(x, y) is the spherical distance (angle) between x = (ψx , ζx ) and y = (ψy , ζy ):
d(x, y) = cos−1 {sin ζx sin ζy + cos ζx cos ζy cos(ψx − ψy )} .

(2)

The spectral representation of a positive definite function in the form of (1) is also given in [31].

2.1

Intrinsic random functions

Let C(S 2 ) be the space of real-valued continuous functions on the sphere. By the Riesz representation theorem ([30], Chapter 13), the dual of C(S 2 ) is the set of all finite regular signed Borel
measures on S 2 , denoted by Λ. Hence, for f ∈ C(S 2 ) and λ ∈ Λ, we define
Z
f (x)λ(dx),
f (λ) =
S2

and
f (gλ) =

Z

f (gx)λ(dx).
S2

for a rotated measure gλ, where g ∈ G.
3

Given an integer κ > 0, a measure λ ∈ Λ is called an allowable measure of order κ if it annihilates
the spherical harmonics of order less than κ, that is,
Z
Ylm (x)λ(dx) = 0, 0 ≤ l < κ,

|m| ≤ l.

S2

Let Λκ be the class of all allowable measures of order κ. It is obvious that Λκ+1 ⊂ Λκ ⊂ Λ. Next,
we establish the invariance of such allowable measures.

Proposition 1. For λ ∈ Λκ , the rotated measure gλ is invariant, that is, gλ ∈ Λκ .
Proof. For any g ∈ G, we have ([28] and [19], Ch. 3),
Ylm (gx)

=

l
X

′

S(l, m, m′ )Ylm (x),

x ∈ S 2,

m′ =−l

where S(l) = {S(l, m, m′ )}m,m′ =−l,...,l is a (2l + 1) × (2l + 1) orthogonal matrix that is independent
of x. Let λ ∈ Λκ , then Ylm (λ) = 0 for l < κ and |m| ≤ l. This leads to
Ylm (gλ) =

Z

S2

Ylm (gx)λ(dx) =

l
X

S(l, m, m′ )

m′ =−l

Z

′

S2

Ylm (x)λ(dx) = 0.

Proposition 1 shows that an allowable measure is invariant with respect to rotations of the
sphere, hence the set Λκ is closed for rotations. In short,
Ylm (λ) = 0 and Ylm (gλ) = 0,

for λ ∈ Λκ , l < κ, |m| ≤ l, and g ∈ G.

Parallel results on the circle are obtained in [13]. Such invariance is essential in developing IRFs,
see [20, 21].

Now, we extend the definitions of f (λ) and f (gλ) to the continuous random process Z(x). If
Z(x) is with the probability measure space (Ω, A, P ), we denote
Z
Z(x)λ(dx),
Z(λ) =
S2

that is, Z(λ) maps Λ to L2 (Ω, A, P ), and
Z(gλ) =

Z

Z(gx)λ(dx).
S2

4

One can define the inner product
hZ(λ1 ), Z(λ2 )i =

Z

S2

Z

S2

cov(Z(x), Z(y))λ1 (dx)λ2 (dy).

Then, λ → Z(λ) is a continuous mapping, and if λ1 = λ2 , Z(λ1 ) = Z(λ2 ) in the sense of norm
induced by this inner product.

Definition. For an integer κ, a continuous random process Z(x) on the sphere is called an Intrinsic
Random Function of order κ (IRFκ), if for any λ ∈ Λκ , the process is invariant with respect to any
g ∈ G. That is, E(Z(λ)) = E(Z(gλ)), and
cov(Z(λ1 ), Z(λ2 )) = cov(Z(gλ1 ), Z(gλ2 )),

λ1 , λ2 ∈ Λκ .

The following theorem characterizes the IRF on the sphere and reveals its connection to the
stationary process through low-frequency truncation.

Theorem 1. For an integer κ ≥ 1, a continuous random process Z(x) on the sphere is an IRFκ if
and only if its frequency-truncated process Zκ (x) is stationary, where
Zκ (x) =

∞ X
l
X

Zl,m Ylm (x).

l=κ m=−l

Proof. “⇐”, for an integer κ ≥ 1, assume the frequency-truncated process Zκ (x) is stationary.
Note that
Z(x) =
and Z(λ) =

R

l
X X

Zl,m Ylm (x) + Zκ (x),

l<κ m=−l

S2

Z(x)λ(dx) for λ ∈ Λκ . We have
Z(λ) =

Z

l
X X

S2

=

Z

S2

Zl,m Ylm (x) + Zκ (x) λ(dx)

l<κ m=−l

l
X X

!

!

Zl,m Ylm (x)

l<κ m=−l

λ(dx) +

Z

S2

Zκ (x)λ(dx).

The first term is with finite summation of Ylm (λ) and it vanishes since λ ∈ Λκ . Therefore,
Z(λ) =

Z

S2

Zκ (x)λ(dx) = Zκ (λ),

5

λ ∈ Λκ .

The expectation of Zκ (λ) is
E(Zκ (λ)) = E

Z

S2

Zκ (x)λ(dx).

Since λ is a finite measure and Zκ (x) is continuous in quadratic mean on the sphere, by Fubini’s
Theorem, we exchange the expectation and integral and obtain
Z
E(Zκ (x))λ(dx).
E(Zκ (λ)) =
S2

Similarly, for λ ∈ Λκ and g ∈ G,
Z(gλ) = Zκ (gλ),

and E(Zκ (gλ)) =

Z

S2

E(Zκ (gx))λ(dx).

Then, following that Zκ (x) is stationary with E(Zκ (x)) = E(Zκ (gx)), we have E(Zκ (gλ)) =
E(Zκ (λ)) and arrive at
E(Z(λ)) = E(Z(gλ)),

λ ∈ Λκ .

For the covariance part, when λ1 , λ2 ∈ Λκ , we have cov(Z(λ1 ), Z(λ2 )) = cov(Zκ (λ1 ), Zκ (λ2 ))
and cov(Z(gλ1 ), Z(gλ2 )) = cov(Zκ (gλ1 ), Zκ (gλ2 )). Similarly, Zκ (x) is stationary, which implies
cov(Zκ (λ1 ), Zκ (λ2 )) = cov(Zκ (gλ1 ), Zκ (gλ2 )). This leads to
cov(Z(λ1 ), Z(λ2 )) = cov(Z(gλ1 ), Z(gλ2 )),

λ1 , λ2 ∈ Λκ , g ∈ G.

Therefore, Z(x) is an IRFκ.
“⇒”, assume that Z(x) is an IRFκ, κ ≥ 1. For any y ∈ S 2 , consider the following measure
λy (dx) = δy (dx) −

l
X X

Ylm (y)Ylm (dx),

l<κ m=−l

where δy (dx) is the Dirac measure. First, for l0 < κ, |m0 | ≤ l0 ,
0
Ylm
(λy )
0

=

Z

S2

0
= Ylm
(y) −
0

0
Ylm
(x)
0

δy (dx) −

l
X X

l<κ m=−l

X

l
X

Ylm (y)

l<κ m=−l

Z

S2

!

Ylm (y)Ylm (dx)

0
Ylm
(x)Ylm (dx).
0

By orthogonality of the associated Legendre polynomials, the second term in the last equation
0
0
reduces to Ylm
(y). Therefore, Ylm
(λy ) = 0, that is, λy annihilates all associated Legendre poly0
0

nomials with order less than κ, which implies
λy (dx) ∈ Λκ .
6

For Z(x), we have
Z(λy ) = Z(y) −

l
X X

Zl,m Ylm (y),

l<κ m=−l

which is exactly the low-frequency truncated process at y, i.e.,
Z(λy ) = Zκ (y).
In the same way,
Z(gλy ) = Zκ (gy), g ∈ G.
Therefore, for any x, y ∈ S 2 , g ∈ G, we have cov(Zκ (x), Zκ (y)) = cov(Z(λx ), Z(λy )) and cov(Zκ (gx), Zκ (gy)) =
cov(Z(gλx ), Z(gλy )). Note that Z(x) is an IRFκ with cov(Z(λx ), Z(λy )) = cov(Z(gλx ), Z(gλy ))
for λx , λy ∈ Λκ , we obtain
cov(Zκ (x), Zκ (y)) = cov(Zκ (gx), Zκ (gy)).
Similarly, for λ ∈ Λκ ,
E(Zκ (x)) = E(Zκ (gx)).
That is, Zκ (x) is stationary.
When κ = 0, the IRF0 is the same as a stationary process. This is slightly different from [20],
where IRF(−1) is a stationary process. For notation simplification, we assume κ ≥ 1 throughout
the rest of this manuscript. Huang et al. [13] introduced IRFs on the circle and showed that an
IRF on the circle is characterized by its low-frequency truncated process. Theorem 1 shows that
the parallel result holds on the sphere and we further conjecture that this is true for hyper-spheres.

Remark 1. In Euclidean spaces, IRFs are associated with differential operations [20]. For example, the (κ + 1)th derivative of a differentiable IRFκ on the real line is stationary ([2], Chapter 4).
Theorem 1 indicates that low-frequency truncation operation of an IRFκ on the sphere becomes
stationary. The implication of this result for splines on the sphere is discussed in Section 2.4.

Let Z(x) be an IRFκ, then its low-frequency truncated process Zκ (x) is stationary (Theorem
1). We have [26, 40]
cov(Zl,m , Zl′ ,m′ ) = al I(l, l′ )I(m, m′ ),
7

l, l′ ≥ κ, |m| ≤ l, |m′ | ≤ l′ .

Hence, the covariance of Zκ (·) is
cov(Zκ (x), Zκ (y)) =

∞
X

al

l=κ

l
X

Ylm (x)Ylm (y).

m=−l

By the addition theorem of Legendre polynomials, we arrive at
cov(Zκ (x), Zκ (y)) =

∞
X
2l + 1
l=κ

4π

al Pl (cos(d(x, y))).

(3)

This covariance function plays an essential role in universal kriging on the sphere in Section 2.2.
In addition, it offers the base in constructing reproducing kernel Hilbert space (Section 2.3). We
denote this function by φ(·), and name it the intrinsic covariance function of an IRFκ on the sphere.

Remark 2. A more abstract view of IRFs in Euclidean spaces is based on equivalence class ([2],
Chapter 4). Here, we extend this discussion to IRFs on the sphere. Let Z(x) be an IRFκ and
Al,m , l < κ, |m| ≤ l be random variables, a new random process
Z ∗ (x) = Zκ (x) +

l
X X

Al,m Ylm (x),

(4)

l<κ m=−l

is also an IRFκ, where Zκ (x) is the low-frequency truncated process of Z(x). In fact, for λ ∈ Λκ ,
we have Ylm (λ) = 0 for l < κ, |m| ≤ l and
∗

Z (λ) = Zκ (λ) +

l
X X

Al,m Ylm (λ) = Zκ (λ).

l<κ m=−l

Therefore, Z ∗ (x) and Z(x) are both IRFκ and they share the same truncated process Zκ (x). All
random functions sharing the same truncated process Zκ (x) form an equivalence class. A natural
choice of the representation of this equivalence class is Zκ (x), and all others are of form of (4).

2.2

Universal kriging

One of the primary applications of developing IRFs is universal kriging. Assume that Z(x) is an
IRFκ, its intrinsic covariance function is φ(·), and the mean function is
E(Z(x)) =

l
X X

l<κ m=−l

8

βl,m Ylm (x),

(5)

where βl,m are the coefficients. Suppose that the data
{(xi , wi ), i = 1, . . . , n},

n > (2κ − 1),

x1 , . . . , xn ∈ S 2

(6)

are observed from this process with uncorrelated measurement errors
W (x) = Z(x) + ǫ(x),

x ∈ S2,

where ǫ(·) is white noise with mean zero and variance σ 2 .
To obtain the best linear unbiased estimator at x0 ∈ S 2 , following [5, 13], we start with a linear
estimator
Ẑ(x0 ) = η T w ,
e e
T
where w = (w1 , . . . , wn ) is the data vector and η = (η1 , . . . , ηn )T is to be determined. For notation
e
e
simplicity, we rearrange the lower spherical harmonics {Ylm (x), l < κ, |m| ≤ l} and denote them by

{qν (x), ν = 1, . . . , 2κ − 1}. The unbiasedness leads to

ηT Q = qT ,
e
e

where

Q = {qν (xi )}n×(2κ−1) ,
Hence,

n
X

q = (q1 (x0 ), . . . , q2κ−1 (x0 ))T .
e

ηi qj (xi ) = qj (x0 ),

j = 1, . . . , 2κ − 1,

i=1

which implies that the discrete measure
n
X

!

ηi δ(xi )

i=1

− δ(x0 ) ∈ Λκ .

The squared prediction error can be shown to be

where

E{Ẑ(x0 ) − Z(x0 )}2 = σ 2 η T η + η T Ψη − 2η T φ + φ(0),
e e
e e e e
Ψ = {φ(d(xi , xj ))}n×n ,

φ = (φ(d(x1 , x0 )), . . . , φ(d(xn , x0 )))T .
e
9

The goal of universal kriging is to minimize the squared prediction error, subject to the unbiasedness
constraints. Let ρ be a Lagrange-multiplier vector of size (2κ − 1), we minimize the following with
e
respect to both ρ and η
e
e
M (η , ρ) = σ 2 η T η + η T Ψη − 2η T φ + φ(0) + 2(η T Q − q T )ρ.
e e
e
e e
e e
e e e e

One, then, can derive the universal kriging formula as

 (Ψ + σ 2 I)η + Qρ = φ,
(7)
e
e e
 QT η = q .
e e
It is clear from this kriging formula that there is no need to estimate the coefficients {βl,m } in (5),

showing the advantage of universal kriging [20]. Note that when κ = 1, universal kriging reduces to
ordinary kriging [5] and the intrinsic covariance function φ(·) relates to the semi-variogram directly,
see [11, 13].

2.3

Splines and kriging

The goal of kriging is to find a linear unbiased predictor at an un-sampled location. It can also be
viewed as surface prediction ([5], Chapter 3), which has close ties to smoothing splines [37]. The
connection between splines and kriging have been extensively discussed in the literature, including
[22, 38, 18, 3, 4, 5, 37, 36, 16, 17, 7, 13]. Furrer and Nychka [7] showed that, given a covariance
function, one can construct a reproducking kernel and obtain a general spline estimate in Euclidean
spaces.
Based on [14], Huang et al. [13] established the connection between IRF and RKHS, and
provided the corresponding smoothing formula for universal kriging on the circle. Now, we extend
this discussion to the sphere. Note that, Levesley et al. [14] constructed the RKHS based on a
∗
∗
non-negative integer κ∗ and a sequence {βl∗ }∞
l=κ . While both κ and βl play essential roles in [14],

their motivation and interpretation are unclear. In this section, we demonstrate that κ∗ is exactly
the order of an IRF, and βl∗ are closely related to the IRF’s intrinsic covariance fucntion φ(·).
Now, following [14], we formally introduce RKHS. Based on the intrinsic covariance function
φ(θ) of an IRFκ, if φ(θ) is of the form (3), we define a function space
(
)
∞ X
l
∞
l
X
X
X
1
Xκ = f (x) =
al,m Ylm (x), x ∈ S 2 :
|al,m |2 < ∞ .
al
l=0 m=−l

l=κ

10

m=−l

Then, for f, g ∈ Xκ , a semi-inner product is defined as
∞
l
X
1 X
al,m,f al,m,g .
al

hf, giκ =

l=κ

(8)

m=−l

There is a nil space for this semi-inner product
N = span{Ylm (x), 0 ≤ l < κ, m = −l, . . . , l},

(9)

with dim(N ) = dN = 2κ − 1. Let {τ1 , . . . , τdN } ∈ S 2 be a set of distinct points. Then, the inner
product
hf, gi =

dN
X

f (τν )g(τν ) + hf, giκ

ν=1

is well defined and Xκ can be shown to be complete [14] with respect to the norm induced by the
inner product (8). In addition, there exist p1 (·), . . . , pdN (·) ∈ N such that pν (τµ ) = I(ν, µ) for
1 ≤ ν, µ ≤ dN . As discussed in [14], the space Xκ is a Hilbert space in which point evaluations are
continuous linear functionals. Therefore, for x, y ∈ S 2 , there exists a reproducing kernel
H(x, y) = φ(d(x, y)) −

dN
X
{φ(d(x, τν ))pν (y) + φ(d(y, τν ))pν (x)}
ν=1

+

dN X
dN
X

φ(d(τν , τµ ))pν (x)pµ (y) +

ν=1 µ=1

dN
X

pν (x)pν (y).

ν=1

Hence, given the intrinsic covariance function of an IRFκ, we construct a RKHS. Conversely,
it is obvious that H(x, y), being a reproducing kernel, is positive definite. Given this property and
from the Kolmogorov existence theorem, there exists a Gaussian random process with H(x, y) as
its covariance function. Suppose we construct a random process Z(x) whose covariance function is
H(x, y). We now show that Z(x) is an IRFκ with the intrinsic covariance function φ(d(x, y)) for
x, y ∈ S 2 . Note that measures in Λκ annihilate any functions in nil space N , that is, for λ ∈ Λκ ,
pν (λ) = 0, ν = 1, . . . , dN . Therefore,
cov(Z(λ1 ), Z(λ2 )) =

Z

S2

Z

S2

φ(d(x, y))λ1 (dx)λ2 (dx),

λ1 , λ 2 ∈ S 2 ,

This covariance is rotation invariant, hence, Z(x) is an IRFκ with φ(·) as its intrinsic covariance
function.
We can see from this development that κ∗ in [14] is exactly the order of an IRF, and βl∗2 = 1/al .
Such interpolation helps understand the RKHS theory on the sphere in [14]. Furthermore, we will
11

show that the smoothing formula in [14] finds its counterpart in universal kriging. Given observed
data {(xi , wi )}ni=1 in (6), the smoothing problem is to find a function fα (x) ∈ Xκ such that it
minimizes [34, 14]

n
X
(wi − f (xi ))2 + αkf k2κ ,

(10)

i=1

where α > 0 is the smoothing parameter and k · kκ is induced by the semi-inner product (8). The
minimizer can be shown to be
fα (x) =

dN
X

bν qν (x) +

ν=1

n
X

ci φ(d(xi , x)),

(11)

i=1

where b = (b1 , . . . , bdN )T and c = (c1 , . . . , cn )T satisfying
e
e

 (Ψ + αI)c + Qb = w,
e
e e
 QT c = 0,
e e

(12)

where Ψ and Q are the same as given in (7). To show the connection between this smoothing formula
and universal kriging, note that the smoothing formula (11) at an unobserved point x0 ∈ S 2 can
be re-written in the following manner

φ
fα (x0 ) = (cT , bT )  e  ,
e e
q
e
where φ and q are the same in (7). Solve for (cT , bT ) in (12) and we have
e e
e
e

−1 

Ψ
+
αI
Q
φ
  e .
fα (x0 ) = (w T , 0TdN ×1 ) 
T
e e
Q
0
q
e
Therefore, fα (x0 ) can be re-written as


∗
η
fα(x0 ) = (w T , 0T )  e  ,
e e
ρ∗
e
where

 
 



−1

or

η∗
Ψ + αI Q
 e =
 
ρ∗
QT
0
e

 (Ψ + αI)η ∗ + Qρ∗ = φ,
e
e
e
 QT η ∗ = q .
e
e
12

φ
e ,
q
e

(13)

Remark 3. These two equations (7) and (13) are exactly the dual formulas of universal kriging
[5, 2], but expressed using the intrinsic covariance function of an IRF on the sphere introduced
in this paper. Usually, universal kriging is viewed as a linear estimator of observed data and the
smoothing formula as the sum of lower spherical harmonic trends and the linear combination of
the intrinsic covariance. Based on equations (7) and (13), these two views are essentially the same.
The connection between universal kriging and the smoothing formula on the sphere is, thus, obvious. In so doing, we find that the lower monomials in Euclidean spaces need to be replaced by
lower spherical harmonics. This echoes the low-frequency truncation that replaced the differential
operations in Euclidean spaces, see Remark 1.

Remark 4. The smoothing parameter α in (13) equals the noise variance σ 2 in (7), and plays
the same role in prediction. For example, in the smoothing formula (11), when α increases to
infinity, the minimization procedure demands kf kκ approaches zero, which shows that c → 0, and
e
the smoothing formula reduces to a spherical harmonic regression. In kriging, when σ 2 → ∞,

the noise overwhelms the spatial dependency, the process becomes essentially uncorrelated. The

squared prediction error is dominated by σ 2 η T η . Then, universal kriging reduces to minimize η T η
e e
e e
subject to the unbiasedness restriction, which also leads to the same spherical harmonic regression
prediction. When both α and σ 2 decrease to zero, both smoothing and kriging result in exact
interpolation.

2.4

Splines on the sphere and thin-plate splines

Splines on the sphere are commonly used for surface fitting [35, 37, 27]. Given the data {(xi , wi ), i =
1, . . . , n} in (6), Wahba [35] introduced splines on the sphere as a function to minimize the following
n
X
(wi − f (xi ))2 + αJ(f ),
i=1

where the penalty term is through the Laplacian operator ∆ on the sphere and for m even,
Z 
2
∆m/2 f sin ζdζdψ.
Jm (f ) =
S2

The spline estimator is shown to be [35]

fα (x) = d +

n
X
i=1

13

ci K(x, xi ),

(14)

where
K(x, y) =

∞
1 X 2l + 1
Pl (cos(d(x, y))),
4π
lm (l + 1)

x, y ∈ S 2 .

l=1

In Euclidean spaces, the integer m plays a significant role and relates to the order of IRFs in kriging
[16]. While K(x, y) here takes the place of φ(·) in (11), the integer m only alters the form of K(x, y)
and loses its connection to the order κ in IRF on the sphere. Furthermore, while the nil space in
(10) is spanned by the lower spherical harmonics of orders up to κ, splines on the sphere in (14)
only include a constant term. By viewing splines on the sphere through the smoothing formula
(11), we now see that splines on the sphere act simply as ordinary kriging. This observation echoes
Remark 1 and reveals that splines with derivative penalty have limited applications on the sphere
(Huang et al. [13] observed the similar result on the circle). Then, a more appropriate smoothing
(spline) model on the sphere is the smoothing formula (10) in Section 2.3. A general approach for
hyper-spheres can be found in [34].
For surface fitting in multi-dimensional Euclidean spaces, thin-plate splines are often used [6,
23, 37]. Some applications [25] have extended thin-plate splines to the sphere by mapping longitude
and latitude directly on two-dimensional Euclidean spaces. It shall be warned that such practice
is flawed. The radial basis kernel of thin-plate splines are [6, 23]
E(x, y) = E(d∗ (x, y)) = |d∗ (x, y)|2 log |d∗ (x, y)|,
where d∗ (x, y) is the Euclidean distance between x and y. When thin-plate splines are obtained in
Euclidean spaces, E(·) is shown [37] to be conditionally positive definite of order 2. However, when
x, y are on the sphere, the distance between x and y is the spherical distance d(x, y) in (2). To
obtain a valid thin-plate spline on the sphere, E(d(x, y)) needs to be conditionally positive definite
of order 2. Menengatto and Peron [24] showed that the coefficients of E(d(x, y)) in the Legendre
polynomials expansion must be non-negative for all orders higher than or equal to 2. However,
direct computation shows that E(·) yields negative coefficients for such expansion on the sphere.
For example, the fourth-order coefficient is negative. Hence, thin-plate splines are not directly
applicable for surface fitting on the sphere.

14

2.5

Extensions

An IRFκ is characterized by its truncated process, while the mean function falls in nil space
(9). This is a direct extension of IRFs in Euclidean spaces. It is important and noteworthy that
the spherical harmonics are orthogonal, while the monomials used in Euclidean spaces are not
necessarily so. Therefore, in developing the notion of IRFs, the allowable measure can also be
defined to annihilate the spherical harmonics of orders in a set {i1 , . . . , iκ } of integers, instead of
just {1, . . . , κ}. Therefore, paralleling to Theorem 1, a continuous random process in the sense of
quadratic mean is an IRF(i1 , . . . , iκ ) on the sphere if and only if
Z(i1 ,...,iκ ) (x) =

X

l
X

Zl,m Ylm (x)

l∈{i
/ 1 ,...,iκ } m=−l

is stationary. This is a more general result that includes the IRFκ developed in Section 2.1 as a
special case and may offer more flexibility in modeling the mean function of an IRF on the sphere.

References
[1] N Aronszajn. Theory of reproducing kernels. Transactions of the American Mathematical
Society, 68:337–404, 1950.
[2] J. Chilès and P. Delfiner. Geostatistics: Modeling Spatial Uncertainty. Wiley, New York, 2
edition, 2012.
[3] N. Cressie. Geostatistics. The American Statistician, 43:197–202, 1989.
[4] N. Cressie. Reply to “comment on cressie” by G. Wahba. The American Statistician, 44:256–
258, 1990.
[5] N. Cressie. Statistics for Spatial Data, revised ed. Wiley, New York, 1993.
[6] J. Duchon. Splines minimizing rotation-invariant semi-norms in Sobolev spaces. In Constructive Theory of Functions of Several Variables, pages 85–100. Springer-Verlag, Berlin, 1977.
[7] E. Furrer and D. Nychka. A framework to understand the asymptotic properties of kriging
and splines. Journal of the Korean Statistical Society, 36:57–76, 2007.

15

[8] I.M. Gel’fand and N.Y. Vilenkin. Applications of Harmonic Analysis (Generalized Functions),
volume 4. Academica Press, New York, 1964.
[9] T. Gneiting. Strictly and non-strictly positive definite functions on spheres. Bernoulli, 19:1327–
1349, 2013.
[10] M. Heaton, M. Katzfuss, C. Berrett, and D. Nychka. Constructing valid spatial processes on
the sphere using kernel convolutions. Environmetrics, 25:2–15, 2014.
[11] C. Huang, H. Zhang, and S. Robeson. On the validity of commonly used covariance and
variogram functions on the sphere. Mathematical Geosciences, 43:721–733, 2011.
[12] C. Huang, H. Zhang, and S. Robeson. A simplified representation of the covariance structure
of axially symmetric processes on the sphere. Statistics and Probability Letters, 82:1346–1351,
2012.
[13] C. Huang, H. Zhang, and S. Robeson. Intrinsic random functions and universal kriging on the
circle. Statistics and Probability Letters, 108:33–39, 2016.
[14] Levesley J., W. Light, D. Ragozin, and X. Sun. A simple approach to the variational theory for
interpolation on spheres. International Series of Numerical Mathematics, 132:117–143, 1999.
[15] A. H. Jones. Stochastic processes on a sphere. Annals of Mathematical Statistics, 34:213–217,
1963.
[16] J. Kent and K. Mardia. The link between kriging and thin-plate splines. In F.P. Kelly, editor,
Probability, Statistics, and Optimization: A tribute to Peter Whittle, pages 325–339. Wiley,
Chichester, 1994.
[17] G.F. Laslett. Kriging and splines: An empirical comparison of their predictive performance in
some applications. Journal of American Statistics Association, 89:391–400, 1994.
[18] A.C. Lorenc. Analysis methods for numerical weather prediction. Quarterly Journal of the
Royal Meteorological Society, 112:1177–1194, 1986.

16

[19] D. Marinucci and G. Peccati. Random Fields on the Sphere: Representation, Limit Theorems, and Cosmological Applications. London Mathematical Society, Lecture Notes Series
389, Cambridge, 2011.
[20] G. Matheron. The intrinsic random functions and their applications. Advances in Applied
Probability, 5:439–468, 1973.
[21] G. Matheron.

Comment translater les catastrophes. La structure des F.A.I. ǵenérales.

manuscript. Technical Report N-167, Centre de Géostatistique, Fontainebleau, France, 1979.
[22] G. Matheron. Splines and kriging: their formal equivalence. In Syracuse University Geological
Contributions. Syracuse University, 1981.
[23] J. Meiguet. Multivariate interpolation in arbitrary points made simple. Journal of Applied
Mathematics Physics, 30:292–304, 1979.
[24] V. Menegatto and A. Peron. Conditionally positive definite kernels on Euclidean domains.
Journal of Mathematical Analysis and Applications, 294:345–359, 2004.
[25] M. New, D. Lister, M. Hulme, and I. Makin. A high-resolution data set of surface climate over
global land areas. Climate Research, 21:1–25, 2002.
[26] A. M. Obukhov. Statistical homogeneous random fields on a sphere. Uspekhi Matematicheskikh,
2:196–198, 1947.
[27] S.M. Robeson. Spherical methods for spatial interpolation: Review and evaluation. Cartography and Geographic Information Systems, 24:3–20, 1997.
[28] R. Roy. Processus stochastiques sur la sphère. PhD thesis, Université de Montréal, 1969.
[29] R. Roy. Estimation of the covariance function of a homogeneous process on the sphere. Annals
of Statistics, 1:780–785, 1973.
[30] H.L. Royden. Real Analysis. Prentice Hall, Englewood Cliffs, NJ, 3rd edition, 1988.
[31] I.J. Schoenberg. Positive definite functions on spheres. Duke Mathematics Journal, 9:96–108,
1942.

17

[32] M. L. Stein. Spatial variation of total column ozone on a global scale. Annals of Applied
Statistics, 1:191–210, 2007.
[33] M.L. Stein. Statistical Interpolation of Spatial Data: Some Theory for Kriging. Springer, New
York, 1999.
[34] H.J. Taijeron, A.G. Gibson, and C. Chandler. Spline interpolation and smoothing on hyperspheres. SIAM Journal on Scientific Computing, 15:1111–1125, 1994.
[35] G. Wahba. Spline interpolation and smoothing on the sphere. SIAM Journal on Scientific and
Statistical Computing, 2, 1981.
[36] G. Wahba. Comment on Cressie. The American Statistician, 44:255–256, 1990.
[37] G. Wahba. Spline Models for Observational Data. CBMS-NSF regional conference series in
applied mathematics, Philadelphia, 1990.
[38] G.S. Watson. Smoothing and interpolation by kriging with splines. Mathematical Geology,
16:601–615, 1984.
[39] A.M. Yaglom. Correlation theory of processes with random stationary nth increments (in
Russian). Matematicheskii Sbornik, 37:141–196, 1955. English translation in (1958) American Mathematical Society Translations: Series 2, 8, 87-141. American Mathematical Society,
Providence, R. I., 1958.
[40] A.M. Yaglom. Second-order homogeneous random fields. Fourth Berkeley Symposium on
Mathematical Statistics and Probability, 2:593–622, 1961.
[41] H. Zhang and C. Huang. A note on processes with random stationary increments. Statistics
and Probability Letters, 94:153–161, 2014.
[42] Z. Zhu and Y. Wu. An efficient algorithm for estimation and prediction of a class of convolutionbased spatial nonstationary models. Journal of Computational and Graphical Statistics, 19:74–
95, 2010.

18

