Journal of Math-for-Industry, Vol. 3 (2011A-6), pp. 63–71

RBF interpolation and Gaussian process regression
through an RKHS formulation
D案

Ken Anjyo and J. P. Lewis
Received on February 11, 2011
Abstract.
Radial Basis Function (RBF) interpolation is a common approach to scattered data interpolation.
Gaussian Process regression is also a common approach to estimating statistical data. Both techniques play a central role, for example, in statistical or machine learning, and recently they havebeen increasingly applied in other ﬁelds such as computer graphics. In this survey we describe theformulation of both techniques as instances of functional regression in a Reproducing Kernel Hilbert
Space. We then show that the RBF and Gaussian Process techniques can in some cases be reduced
to an identical formulation, diﬀering primarily in their assumptions on when the data locations andvalues are known, as well as in their (respectively) deterministic and stochastic perspectives. The
scope and eﬀectiveness of the RBF and Gaussian process techniques are illustrated through several
applications in computer graphics.
Keywords. Reproducing Kernel Hilbert Space, Radial Basis Function, scattered data interpolation,
Gaussian process regression,

1.

Both RBF interpolation and Gaussian process regression involve the inverse problem of estimating an unknown
function from a ﬁnite number of data samples. We describe how an RKHS is the appropriate theoretical basis
for both techniques. In particular, for our practical purposes of solving interpolation/regularization problems inan RKHS, we utilize a powerful feature of RKHSs, known
as the Representer Theorem [16]. This theorem provides
eﬀective solutions to function estimation problems in theform of weighted combinations of a kernel function. Givenﬁnite data samples, seeking an unknown function of thissimple form is much easier to justify than a search within
the unrestricted (inﬁnite-dimensional) Hilbert space. Appropriate choices of the kernel function can summarize theknown or assumed characteristics of the particular classof functions, thus leading to sensible regularization of theinverse problem. The RBF and GPR techniques can beseen as reﬂecting diﬀerent approaches to specifying theseassumptions.

Introduction

Radial Basis Function (RBF) interpolation is a common
approach to scattered data interpolation. Gaussian Process regression (GPR) can also be used for scattered datainterpolation and extrapolation. In this paper we interpret both of these techniques in the common framework of
Reproducing Kernel Hilbert Spaces (RKHSs). The RKHS
interpretation helps reveal the commonalities of these twotechniques as well as clarify their diﬀerences. The practicalrelevance of the techniques will also be illustrated throughseveral examples in computer graphics.
RBF Interpolation was popularized in the late 1980s,
initially in the machine learning community [4]. More recently RBF interpolation has had numerous applicationsin other ﬁelds. For example, in computer graphics RBF interpolation has been used for reconstructing models fromscanned data [5], for deﬁning the skin motion of characters [10, 8], for interpolating volume medical scans [15], for
editing lighting and shading on anime characters [17], andother purposes.

In addition, linear operations within the RKHS can beinterpreted as nonlinear operations in the input space (thisis known as the “kernel trick” in the machine learning community). This enables problems to be solved as linear systems while still having the some of the power of a nonlinear
approach.

Gaussian process regression was initially explored byKolmogorov and Wiener in the 1940s. It was later developed in the geostatistics community where it goes by thename Kriging [6]. Gaussian process techniques have been
widely applied, for example in computer graphics Gaussian
processes have been used for modeling the motion of human characters [7, 14], for synthesizing terrains [9] and forgenerating novel variations of human drawings [2].

The RBF and GPR techniques and the notion of RKHS
are deﬁned in the next section. The RKHS formulationof RBF interpolation is presented in section 3. Section
4 derives GPR from functional regression in an RKHS.
63

E案

64

Journal of Mathematics for Industry, Vol. 3 (2011A-6)

Section 5 applies the preceding discussion to establish an
equivalence between RBF interpolation and GPR.

2.
2.1.

Definitions and Theory

Reproducing Kernel

Let E be an abstract set, and H be a Hilbert space consisting of the (real-valued1 ) functions deﬁned on E, with
the inner product ⟨, ⟩.
Deﬁnition 1. The function K : E × E → R is called a
reproducing kernel of H, if it satisﬁes the following conditions2 :1. For any ﬁxed y ∈ E, K(x, y) belongs to H as a function of x on E.
2. For any f ∈ H, we have f (y) = ⟨f (x), K(x, y)⟩x .
Deﬁnition 2. If Hilbert space H has the above kernel K,then H is referred to as a reproducing kernel Hilbert space(RKHS).
The following proposition will be used in characterizing
the reproducing kernel in the next section.
Proposition 1. For the reproducing kernel K, we have:(1)

K(y, z) = ⟨K(x, y), K(x, z)⟩x

Proof. From condition 1 of the reproducing kernel, we
have K(y, z) ∈ H for a ﬁxed z. Then, by putting f (y) =K(y, z) in condition 2, we haveK(y, z) = ⟨K(x, z), K(x, y)⟩x = ⟨K(x, y), K(x, z)⟩x .
The next theorem is well known as a classical result in the
reproducing kernel theory (N. Aronszajn [1], S. Bergman[3]).
Theorem 1. Hilbert space H(E) has a reproducing kernel,
if and only if the following condition is satisﬁed:
For any y ∈ E, there exists a positive constant C = Cy ,such that(2)|f (y)| ≤ Cy ∥f ∥, for any f ∈ H.

Proof. [only if part] The assertion follows from Schwarz’
inequality. [if part] For a ﬁxed y ∈ E, let us consider
the linear functional δy : H(E) → R, which is deﬁned
as δy (f ) := f (y), for f ∈ H(E). Condition (2) thenmeans that δy is continuous. The assertion thus follows
from Riesz’ Theorem [21]. That is, for a ﬁxed y ∈ E, there
exists a unique function Ky (x) ∈ H(E) such that δy (f ) =⟨f, Ky ⟩ for any f ∈ H(E). Then we put K(x, y) = Ky (x).
Relating to the above theorem, we note that the reproducing kernel K is uniquely determined for an RKHS H(E)(also see Theorem 3 in the next section).
1 This article treats only real-valued functions, while the extension
to complex-valued functions is obvious.
2 In condition 2, the inner product ⟨, ⟩ means that we get thexinner product value of the two functions with variable x.

RKHS Examples
1. R: Let E be {1}. Here we identify H(E) with R. An
arbitrary element of H(E) is a map f : E ≡ {1} →R. Specifying f ∈ H(E) therefore means specifying areal number x with f (1) = x. Let the inner product⟨, ⟩ for H(E) be the ordinary multiplication in R :⟨x, y⟩ = x · y. We deﬁne K : E × E → R as K(1, 1) :=1. It is then easy to see K satisﬁes condition 1 in thedeﬁnition of the reproducing kernel. As for condition
2, we have: ⟨f (1), K(1, 1)⟩ = f (1)1 = f (1).
2
2. l2 : Let a = (aj )∞j=1 ∈ l . Then a deﬁnes a map α :N → R with α(i) := ai (i ∈ N). We thus identify a ↔
α. ∑By setting E = N, we have H(N) ≡ {α : N →∞R| i=1 |α(i)|2 < ∞} ∼= l2 with its kernel function K
as being K(i, j) = δij .

3. Let A be an n-th order, symmetric, and positive semideﬁnite matrix. Then A(Rn ) is RKHS and its reproducing kernel is A (see the discussions in the next
section):A(Rn ) ≡ {Ax ∈ Rn | x ∈ Rn }.
2.2.

Fundamental Properties

Theorem 2. Let K : E × E → R be the kernel function of
RKHS H(E). Then K satisﬁes the following properties:
1. K is symmetric: K(x, y) = K(y, x) for any x, y ∈ E.
2. K is positive semi-deﬁnite: For any (x1 , · · · , xn )TnTn∈∑nE and (a1 , a2 , · · · , an ) ∈ R , we always havei,j=1 ai aj K(xi , xj ) ≥ 0.
Proof. (1) in Proposition 1 says that K is symmetric,
since the inner product itself is symmetric. For (x1 , · · · , xn )T∈ E n and (a1 , a2 , · · · , an )T ∈ Rn , it follows from (1) thatn∑ai aj K(xi , xj ) =

i,j=1

n
∑

ai aj ⟨K(x, xi ), K(x, xj )⟩x

i,j=1

=

⟨ n
∑

aj K(x, xi ),

n
∑

i=1

j=1

k=1

x

⟩
aj K(x, xj )

° n
°2
°∑
°
°
°
=°
ak K(x, xk )° ≥ 0.
°°x

The following theorem says that RKHS is constructed by
specifying a symmetric, positive semi-deﬁnite function. Itshould also be noted that the proof is constructive so thatit might be useful even in our practical situations.
Theorem 3. Suppose that K is a symmetric, positive semideﬁnite function on E×E. Then there exists a Hilbert space
H that has K as its reproducing kernel.
∑lSketch of the proof. We put F := { i=1 αi K(x, xi )| l ∈N, αi ∈ R, xi ∈ E}. By deﬁning addition and multiplication by constant as usual, we can make F a vector space.

65

Ken Anjyo and J. P. Lewis

Also we can∑mintroduce the inner product∑n for the expressions f = i=1 αi K(x, xi ) and g = j=1 βj K(x, xj ) ∈ F∑m ∑n
as follows: ⟨f, g⟩ := i=1 j=1 αi βj K(xi , xj ) ∈ R. Since
K is positive semi-deﬁnite, it follows that ⟨f, f ⟩ ≥ 0. It iseasy to see that F is a pre-Hilbert space. Next we set H as
the completion3 of F . Then, with g(x) ∑≡ gy (x) = K(x, y),
m
we
have
⟨(f
(x),
K(x,
y)⟩
=
⟨f,
g
⟩
=
x
y
i=1 αi K(xi , y) =
∑m
α
K(y,
x
)
=
f
(y),
for
any
f
∈
F
.
This also holdsi
i−1 i
for any f of H, because, for any Cauchy sequence {fm } of
F , we have(3)

|fm (y) − fn (y)| ≤ |⟨(fm − fn )(x), K(x, y)⟩x |≤ ∥fm − fn ∥x · ∥K(x, y)∥x√
= ∥fm − fn ∥ · K(y, y).

This means that {fm (y)} ⊂ R converges at any y. H
therefore contains f = limn→∞ fn , because f also satisﬁes
condition 2 in the deﬁnition of reproducing kernel.
One more remark is about an RKHS that has the complete orthonormal system (CONS). Let us describe the reproducing kernel K with CONS ≡ {φj }∞j=1 . From condition 1 of the reproducing kernel, we ﬁrst getK(x, y) =∞∑αi (y)φi (x).

i=1Then, taking f (y) = φi (y) in condition 2, we haveφi (y) = ⟨φi , K(·, y)⟩ = ⟨φi ,
=

∞
∑

∞
∑

αk (y)φk ⟩

k=1
∞
∑

αk (y)⟨φi , φk ⟩ =

k=1

E×E

This yields that κ(f ) ∈ L2 and that κ is a continuous
linear operator, known as Hilbert-Schmidt integral operator.
According to Mercer’s theorem,
we then have the eigen∑
decomposition: K(x, x′ ) = ν≥1 λν ϕν (x)ϕν (x′ ) where ν ∈N, and λν , ϕν are eigen value and∑eigen functions of κ,
respectively. Assumption (4) yields ν≥1 λ2ν < ∞, so that
we have limk→∞ λk = 0. We now assume that λ1 ≥ λ2 ≥· · · ≥ λn ≥ · · · > 0. We then know that {ϕn }∞n=1 is a CONSof L2 , which consequently gives the following result:
Theorem 4. Let Hκ be the totality of the functions f ∈∑fk2
L2 , satisfyingk≥1 λk < ∞, where fk = (f, ϕk )L2 , referred to as the Fourier coeﬃcient of f w.r.t. {ϕn }∞n=1 .
We then have:
1. H
space with the inner product: ⟨f, g⟩κ =∑κ is afkHilbert
gk<∞,
where f, g ∈ Hκ , with fk = (f, ϕk )L2
k≥1 λkand gk = (g, ϕk )L2 .
2. For g ∈ Hκ , we have(g, ϕν )L2,λν
1
1⟨ϕν , ϕν ⟩κ =(i.e., ∥ϕν ∥κ = √ ) (ν = 1, 2, · · · ).
λν
λν⟨g, ϕν ⟩κ =αk (y)δik = αi (y).

k=1

We therefore haveK(x, y) =∞∑3. K is the reproducing kernel of Hκ : f (x) = ⟨f, K(·, x)⟩κfor any f ∈ Hκ .
φi (x)φi (y).

i=1

2.3.

Then κ can be considered as a linear operator: L2 → L2 ,because, using Schwarz’ inequality, we have
)2
∫
∫ (∫
2
′
|κ(f )(y)| dy =K(x, x )f (x)dx dx′
E
E
E
)
∫ (∫
∫
′ 2
2
≤
|K(x, x )| dx|f (y)| dy dx′
E
E
E
∫∫
∫
′ 2
′
=
|K(x, x )| dxdx
|f (y)|2 dy
E×E
E
∫∫
=
|K(x, x′ )|2 dxdx′ · ∥f ∥2L2 < ∞.

RKHS in L2 space

We brieﬂy describe an inﬁnite-dimensional RKHS in L2 (E),where E is assumed to be a domain in Rn . We refer to L2
instead of L2 (E) for simplicity. The∫ inner product (, )L2
in L2 is then deﬁned as: (f, g)L2 = E f (x)g(x)dx for any
f, g ∈ L2 (E). Supposing that K is a symmetric, positive
semi-deﬁnite function on E × E, we ﬁrst
∫ deﬁne the realvalued function κ(f ) on E : κ(f )(y) := E K(x, y)f (x)dx,
for any y ∈ E. Let us further suppose that∫∫(4)

We skipped the mathematical details and the rigorous
proof of the above theorem. Instead, we should keep inmind the relation between Hκ and L2 through the CONS
derived from the Hilbert-Schmidt operator κ. We also notethat a similar result is obtained in a ﬁnite-dimesional case,where K simply means an n-th order symmetric, positive
semi-deﬁnite matrix and L2 ∼= Rn .
2.4.

Radial Basis Functions

The form(5)fˆ(x) =N∑

wk G(x − xk )k|K(x, x′ )|2 dxdx′ < ∞.

E×E
3 The completion simply means the space that consists of all thelimits of Cauchy sequences of F .

is an example of radial basis interpolation of a function
f (x) known at a set of points xk . The function G() is assumed to be a radially symmetric function of its argument,
i.e. G(r) := ϕ(|r|).

66

Journal of Mathematics for Industry, Vol. 3 (2011A-6)

The linear system (6) of course has a solution for a positivedeﬁnite G().
In other cases the matrix Gi,j ≡ G(xi − xj ) has a signiﬁcant nullspace. A prominent example is spline interpolation of scattered data [13, 18], to be discussed further insection 3. In the spline case the interpolation has the form
(7)fˆ(x) =N∑

wk G(x − xk )+pm−1 (x)kwhere pm−1 (x) ∈ Pm−1 . Here we denote the space of polynomials of degree at most k by Pk .
Figure 1: Schematic diagram of RBF interpolation for posespace deformation [10]. The coeﬃcients of each vertex areinterpolated from example shapes in diﬀerent poses. In theﬁgure, the position of a particular vertex on the fore-arm
is interpolated as a function of several examples situatedin a “pose space” described by the elbow and fore-arm rotations.

Figure 2: An animated human hand created with WPSD
(reproduced from [8]). After obtaining example hand poses
from medical scans, the resulting WPSD model can produce a realistic hand in any pose.
Given known values yk ≡ f (xk ), the weights wk can besolved as a linear system(6)   G(x1 − x1 ) G(x1 − x2 ) G(x1 − x3 ) · · ·w1y1 G(x2 − x1 ) G(x2 − x2 ) w2  y2 ···    G(x3 − x2 ) w3  = y3 ···   
..
..
..
.
.
.
Since the argument to the kernel function G() depends
only on the distance between points, the interpolation takes
the form of a weighted sum of radially symmetric functions.
Common choices of the kernel function G() with ϕ(r)are:
• Gaussian

ϕ(r) = exp(−(r/c)2 )√• Hardy multiquadratic ϕ(r) = r2 + c2 , c > 0
Deﬁnition 3. [19] A function G : Rn → R is positive semideﬁnite if, for any N ∈ N, all sets of N distinct points
xk ∈ Rn , and any w = (w1 , w2 , · · · , wN )T ∈ RN , the∑N ∑Nquadratic form i=1 j=1 wi wj G(xi −xj ) is non-negative.
The function G is called positive deﬁnite if the quadratic
form is strictly positive, for any w ̸= 0 ∈ RN .

RBF Application Examples
1. Pose space deformation (PSD)[10]: This technique isan example-based skinning method. The word “skinning” refers to the skin (surface) deformation of a 3D
character such as human as the underlying skeleton
moves (Figure 1). PSD uses skin surface examples ofvarying posture and interpolates them during animation as a function of the vector of joint angles of theskeleton. The skin is prescribed with its vertices, sothat RBF interpolation is applied to the coeﬃcientsof those vertices. In the notation of (5), x is the vector of joint angles or other parameters describing thecurrent pose of the skeleton, xk is the pose of the kth
example, and fˆ(x) is the interpolated value of one ofthe vertex coeﬃcients. The weights wk are diﬀerent
for each vertex coeﬃcient. Kurihara and Miyata [8]introduced the improved weighted PSD (WPSD) algorithm. In this approach the distance between poses∥xj − xk ∥ is diﬀerent at each vertex, with individual
terms in the distance being weighted by the inﬂuenceof that parameter on the particular vertex. This effectively replaces the single pose space in PSD withsmaller coupled pose spaces at each vertex, thereby
greatly reducing the number of needed training examples. Both PSD and WPSD use RBFs of Gaussian
type (and its nomalized version in [8]).
2. Learning doodles by example: Baxter and Anjyo [2]proposed the concept of a latent doodle space, a lowdimensional space derived from a set of input doodles, or simple line drawings in R2 . This approachﬁrst gives a heuristic algorithm for ﬁnding stroke correspondences between the input drawings, and thenproposes a few latent variable techniques to automatically extract a low-dimensional latent doodle spacefrom the inputs. By the stroke correspondence algorithm, each of the input drawings is represented asa feature vector by combining all the x and y coordinates of each point on each stroke into one vector.
One of the latent variable techniques ﬁrst employs
PCA (Principal Component Analysis) applied to thefeature vectors. The ﬁrst two principal components
are then used to constitute the 2-dimensional latent
doodle space, and, using thin plate spline RBF, we
can obtain new drawings in this 2-d space (Figure 3)

67

Ken Anjyo and J. P. Lewis(a) cartoon face
Figure 4: Synthetic ocean waves generated by Gaussian
process regression with a speciﬁed covariance [9].
is optimal when the error is orthogonal in expectation to
the data:
E [(û(x) − u(x)) yk ] = 0 ∀ k(b) jellyﬁsh drawing
Figure 3: “Doodle” synthesis. The left images were drawn
by an artist; the right images were then synthesized using
the technique in [2].

2.5.

Gaussian Process regression

A Gaussian process is a function from which any ﬁnite
number of samples have a jointly Gaussian distribution.
A Gaussian process u(x) is thus completely speciﬁed byits mean E[u(x)] and covariance function Σu (x1 , x2 ) =E[u(x1 )u(x2 )]. Gaussian process regression (GPR) estimates the values of a Gaussian process given observed values at other locations.
A value u(x) of a Gaussian process at location x can be
estimated as(8)û(x) = σ T Σ−1u ywhere y ∈ Rm is the vector of known data values, andσ = E[u(x), u(xk )] is the covariance between u(x) (thefunction at the point being evaluated) and the function
at the data points xk . If it is assumed that the data y
are observed in the presence of an independent identically
distributed Gaussian observation noise ϵ then (8) takes the
form
(9)

û(x) = σ T (Σu + λI)−1 y

with λ being the variance of the observation noise.
Expression (8) is derived by assuming a linear estimate(10)û(x) = wT y

The weights w for the optimal linear estimate are foundby applying the orthogonality principle, i.e., the estimate

Although the true value u(x) is unknown, substituting (10)and applying the expectations E[yj yk ] = Σ(xj , xk ) andE[u(x)yk ] = Σ(x, xk ) results in (8) with w = Σ−1u σ. Notethat the linear estimate (10) is in fact the optimal estimator (among all linear and non-linear estimators) if theunderlying function is in fact Gaussian [20].
Another useful expression is the variance[ of the regres- ]sion estimate. This can be obtained from E (û(x) − u(x))2
by substituting (10), resulting in
[
]
[
]
(11)
E (û(x) − u(x))2 = E u(x)2 − σ T Σ−1u σ
GPR Application Examples
1. Kriging:
Gaussian process regression is form of Kriging, a term
used in geostatistics to refer to estimation of quantities based on spatial samples and covariance (orvariogram) information [11]. More speciﬁcally, (8) is
identical to simple Kriging, whereas the more commonly used ordinary Kriging adds an additional constraint that the weights w sum to one. In addition
to its use in geostatistical interpolation, Kriging has
recently been applied in computer graphics for interpolating human body poses from motion capture[14].
2. Terrain Synthesis:
Gaussian process regression is also called Wiener interpolation in the signal estimation literature. [9]utilized this interpolation technique in a hierarchical subdivision scheme to incrementally synthesize
random textures and terrains. Whereas the fractal approach to terrain synthesis is restricted to the
1/f p family of power spectra, in Wiener interpolation
an arbitary covariance is speciﬁed, allowing the synthesis of non-fractal landscapes such as ocean waves(Figure 4). Note that the covariance function andpower spectrum are a Fourier transform pair by the68

Journal of Mathematics for Industry, Vol. 3 (2011A-6)

Wiener-Khinchine theorem. The goal in texture andterrain synthesis is to obtain a representative sample from the underlying random process. To achievethis, in the synthesis each variable is ﬁrst estimatedwith (10), but then this estimate is displaced by ani.i.d. (independent identically distributed) Gaussianpseudo-random variable with variance (11).

3.
3.1.

n
It would therefore be nice, if we could have Bm
as then
RKHS in the above theorem. However, Jm cannot be thensquared norm for Bm, as described next.
n
The functional Jm has the following properties:n
1. Jm(f ) = 0 ⇔ f ∈ Pm−1 .
n2. Jm(f ) = (−1)m ⟨f, ∆m f ⟩L2 , where ∆ is Laplacian:∑n ∂ 2∆:= i=1 ∂x
2.

RKHS and Radial Basis Functions
Regularization problem in RKHS

For simplicity, let E = Rn . Suppose that (xi , yi ) are given
as sample points, where yi ∈ R, xi ∈ Rn (i = 1, 2, · · · , N ).
Let us consider the following regularization problem: Finda function f deﬁned on Rn such that
{N
}
∑
2
n
(12)
min
(yi − f (xi )) + λJm (f ) ,fin
Since the null space of Jm
is equal to Pm−1 , we ﬁrst
nrepresent Bm as being the direct sum of the two functionnn⊕ Pm−1 . Then let us solve the regular= Hmspaces: Bmnization problem on Hm. As a result of [12], we thereforehave:nis an RKHS with:
Theorem 5. If m > n2 , then Hm

(18)
n :=
⟨f, g⟩Hm

i=1

where
(13)

∑

m!
⟨Dα f, Dα g⟩L2
α
!α
!
·
·
·
α
!
1
2
n
α +α +···+α =m1n(f ) :=Jm

2n
The regularization term λJm(f ) in (12) prescribes smooth-ness of a solution. Where can we ﬁnd a solution of thisnproblem (12)? According to the deﬁnition of Jm(f ), weshould ﬁnd a solution to (12) in the following space:n(14) Bm:= {f : Rn → R ∪ {±∞}|Dα f ∈ L2 (Rn ),for any α(|α| = m)}.

Now recall the thin plate spline case. We then start with
B22 in (14) to minimize the energy functional
(15)
¯ )
¯ 2
¯
¯
∫ ∫ (¯ 2 ¯2¯∂ φ¯
¯ ∂ φ ¯2 ¯ ∂ 2 φ ¯ 2
¯
¯
¯
¯
¯
¯
F (φ) :=
¯ ∂x2 ¯ + 2 ¯ ∂x1 ∂x2 ¯ + ¯ ∂x2 ¯ dx1 dx2 .
R2
1
2
In the following we want to solve the regularization problem like (12) in RKHS. The main reason for this is the following nice property of RKHS:
Representer Theorem 4 Let H be an RKHS, with its
reproducing kernel K and norm ∥ ∥H . Consider the regularization problem of the following form: Find f ∈ H such
that

min

{N
∑

}
(yi − f (xi ))2+ λ∥f ∥2H.

i=1

The solution f can then be found in the form:(17)f (x) =N∑

With the above theorem, the regularization problem (12)is restated as:
}
{N
∑
2
n
(yi − f (xi )) + λJm (f )
minn
f ∈Bm

(19)
=

i=1

Please

min

i=1

{N
∑

n , p∈Pg∈Hmm−1}{yi − (g(xi ) + p(xi ))}2+ λ∥g∥2Hm
n

.

i=1

We thus know the solution to (19) is represented in theform of(20)

f (x) =N∑

αi K(x, xi ) + pm−1 (x),i=1where pm−1 (x) ∈ Pm−1 . Here we denote the space of polynomials of degree at most k by Pk .
3.2.

RBF as Green’s function

To describe RBF in the above RKHS framework, we ﬁrst
consider the function Ky (x) = K(x, y) 5 , where K(x, y)nis the kernel of Hm. In [12], it is shown that Ky gives a6
weak solution of a certain type of non-homogeneous linear
diﬀerential equation of m-th iterated Laplacian ∆m . Next
let G be a Green’s function in the sense that(21)

αi K(x, xi ).

4 This is one of the variations of the representer theorem.

refer to [16].

nThis also means ∥f ∥2Hm
n = Jm (f ).

n∂mf
Dα f :=α1αn .
2∂x1 ∂xα
2 · · · ∂xn

f ∈H

n

= ⟨(−1)m ∆m f, g⟩L2 .

m!
∥Dα f ∥2L2 ,
α
!α
!
·
·
·
α
!
1
2
n
α +α +···+α =m
1

(16)

2

∑

∆m G(x) = δ(x),5 The function K (x) is the same as K (x) appeared in the proofyyof Theorem 1.
6 The weak solution is a solution in the sense of distributions [21].

69

Ken Anjyo and J. P. Lewis

where δ is the Dirac delta (generalized) function. Then it
is further shown in [12] that we can replace equation (20)by(22)

f (x) =

N∑

We then know that Z is a pre-Hilbert space having theinner product∫⟨Z1 , Z2 ⟩ =Z1 (ω)Z2 (ω)dP (ω) = E[Z1 , Z2 ]Ω

αi G(x − xi ) + p(x),where Z1 , Z2 ∈ Z. The completion of Z is therefore a
Hilbert space.
where p is a polynomial ∈ Pm−1 . This brings us to our
The covariance function Σu (x, y) = E[u(x)u(y)] inducesfamiliar class of radial basis functions:a Kernel function on E ×E, so that Theorem 3 asserts that(23)
there exists RKHS H associated with this kernel function.
{
2m−n
It’s
easy to see that the completion of Z is isometrically
βmn |x|log |x| if 2m − n is an even integer,
G(x) =
2m−nisomorphic
to H [18].
γmn |x|otherwise,
Now again let us assume that the data y = (y1 , y2 , ..., yN )T
are observed in the presence of an independent identicallywhere βmn and γmn are constants.
For example, for the thin plate spline, we have m = n = 2 distributed observation noise. Consider the regularization
so that G(x) = |x|2 log |x| and the polynomial p in (22) is problem (16) in the above H. According to the Represenof the ﬁrst order (linear). Another familiar case is where ter Theorem, the solution f is represented in the form (17).
2m = 2 and n = 3. Then we have G(x) = |x| and a linear Considering the reproducing property (1), the norm ∥f ∥H
can therefore be written in the form
polynomial for (22).
i=1Regularization with RKHS norm Let us consider annother regularization problem,
where, instead of Jm(f ) in∑n(12) and (13), we take m≥0 am Jm (f ) with am being constants and a0 ̸= 0. According to the scenario used in establishing the above theorems, we have the following resultin which the regularization problem is directly solved in an
RKHS.
1. We can ﬁnd a Green’s function for the operator∑(−1)m ∆m .
m≥0∑N
The solution is then given by f (x) = k=1 ck G(x −xk ). Note that this time we do not need a polynomial
term like (22).
2. Gaussian RBF. In a particular case, where am =σ 2m
m!2m (σ > 0), we have the Green’s function as G(x) =c exp(− ∥x∥2σ 2 ).

∥f ∥2H = αT Kαwith αi gathered into a column vector α and K(xi , xj )gathered into a matrix K. We will rewrite the discrete
regularization problem (16) asmin ∥y − Kα∥2 + λαT Kαfwith y being a vector containing data values yi .
Then expanding the square and taking the derivativewith respect to α, we haved
= 0 = −2KT y + 2KT Kα + 2λKα.
dα
Because K is symmetric, it follows that
KK + λKα = Ky
⇔ (K + λI)α = y⇔ α = (K + λI)−1 y.

24.

From RKHS to Gaussian Processes

GPR (8) and (9) can be derived using the RKHS duality[18] using an appropriate choice of kernel. Here we derive(9) through the regularization problem (16)7 .
Let {P, Ω, B} be a probability space. Suppose then thata zero-mean Gaussian process u(x) ≡ u(x, ω) is given on
the probability space, with x ∈ E and ω ∈ Ω. Consider the
space Z consisting of stochastic variables, which is inducedby the Gaussian process:k∑Z={
aj u(xj ); aj ∈ R, xj ∈ E, k ∈ N}.
j=1
7 As for (8), see [18].

Here we assume that K and K + λI are invertible. Now
substitute this expression for α into f = Kα,
f = K(K + λI)−1 y.
This is an expression to synthesize the Gaussian process at
all the locations corresponding to y. For a particular data
point yi = f (xi ) this specializes toyi = [K(xi , x1 ), K(xi , x2 ), · · · , K(xi , xn )](K + λI)−1 y
Then generalizing from location xi to an arbitrary location x and choosing K = Σ−1 we have (9). It is interesting to note that the regularization parameter λ in RKHS
functional regression takes on the role of the observationnoise variance in (9). This is somewhat intuitive in thatthe functional regression requires increased regularizationin the presence of increased noise.

*

70

Journal of Mathematics for Industry, Vol. 3 (2011A-6)

5.

An equivalence between Radial
Basis and Gaussian Process models

The RBF and GPR models initially seem quite diﬀerent.
GPR is a weighted sum of the data, whereas RBF is aweighted sum of a kernel indexed by the distances between
the data. The RBF weights can be computed in advance,
whereas the GPR weights seem to need to be computed at
run time. However, since the problem of ﬁtting functionsin an RKHS to data yields both Radial Basis function andGaussian process models, there should be a direct and simple relationship between the RBF and GPR models. In thissection we state this relationship.
The solution of the system (6) for the RBF weights canbe written
w = G−1 ywith G being the matrix containing G(xi − xj ) in the i, jthposition. Substituting this expression for the weights, the
RBF interpolation (5) can be written in vector-matrix notation asŷ0 = fˆ(x) = rT w = rT G−1 ywith rT = [G(x − x1 ), G(x − x2 ), G(x − x3 ), · · · ]. If wetake the data covariance as the RBF kernel this is exactly
the same expression as (8). This equivalence is reasonable
in some cases. In particular, the Gaussian distribution isboth a prototypical choice of covariance and a standard
choice for the RBF kernel.
With this equivalence in place, the main operational difference between the RBF and GPR techniques is the assumption on “what is known when”. In the RBF techniquethe data values yi are assumed known in advance, whereas
the locations to be interpolated (or extrapolated) are only
known when the interpolation is evaluated. GPR has the
opposite assumption, that the data covariance is known inadvance (at arbitrary locations), but the surrounding data
values yi are not known until run-time.
On the other hand RBF and GPR models are not always
equivalent. There are many choices of covariance that are
not radially symmetric, and many choices of radial basisthat cannot be covariances. A common instance of the latter was discussed in section 3.2, i.e., the radial basis functions generated as Green’s functions of a regularizer (13).
These functions are zero at the origin, which precludes considering them as covariances.

6.

Conclusion

RBF interpolation and Gaussian process regression canboth be expressed in terms of a common underlying RKHS
formulation. Speciﬁcally, both techniques are instances offunctional regression in an RKHS. This common formulation illuminates the role and choice of the kernel, andallows us to identify an exact formal equivalence between
the two techniques under appropriate choices of this kernel. While we wish to make this point of equivalence, thetwo techniques nevertheless have signiﬁcant diﬀerences inmotivation and practice. The common formulation clariﬁes
both similarities and diﬀerences.

Acknowledgments
We thank Kenji Fukumizu, Geoﬀrey Irving and Jean-Marie
Aubry for discussions of several points. This work is partially supported by Japan Science and Technology Agency,
CREST Project.

References[1] N. Aronszajn, Theory of reproducing kernels, Trans.
Amer. Math. Soc. 68 (1950), 337–404.
[2] W. Baxter and K. Anjyo, Latent doodle space, Computer Graphics Forum 25 (2006), 477–485.
[3] S. Bergman, The kernel function and the conformal
mapping, Mathematical Surveys 5 (1950).
[4] D.S. Broomhead and D. Lowe, Multivariable functional interpolation and adaptive networks, Complex
Systems 2 (1988), 321–355.
[5] J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell,
W. R. Fright, B. C. McCallum, and T. R. Evans, Reconstruction and representation of 3d objects with radial basis functions, SIGGRAPH ’01: Proceedings ofthe 28th annual conference on Computer graphics andinteractive techniques (New York, NY, USA), ACM,
2001, pp. 67–76.
[6] N. Cressie, Statistics for spatial data, Wiley, 1993.
[7] Keith Grochow, Steven L. Martin, Aaron Hertzmann,and Zoran Popovic, Style-based inverse kinematics,
ACM Trans. Graph. 23 (2004), no. 3, 522–531.
[8] Tsuneya Kurihara and Natsuki Miyata, Modeling deformable human hands from medical images, Proceedings of the 2004 ACM SIGGRAPH Symposium on
Computer Animation (SCA-04), 2004, pp. 357–366.
[9] J. P. Lewis, Generalized stochastic subdivision, ACM
Transactions on Graphics 6 (1987), no. 3, 167–190.
[10] J. P. Lewis, Matt Cordner, and Nickson Fong, Pose
space deformation: a uniﬁed approach to shape interpolation and skeleton-driven deformation, SIGGRAPH (New York, NY, USA), ACM Press/AddisonWesley Publishing Co., 2000, pp. 165–172.
[11] Christopher Lloyd, Local models for spatial analysis,
CRC Press, 2007.
[12] J. Meinguet, Multivariate interpolation at arbitrary
points made simple, Z. Angew. Math. 30 (1979), 292–
304.

71

Ken Anjyo and J. P. Lewis[13] Charles Micchelli, Interpolation of scattered data: Distance matrices and conditionally positive deﬁnite functions, Constructive Approximation (1986), no. 2, 11–22.
[14] Tomohiko Mukai and Shigeru Kuriyama, Geostatistical motion interpolation, ACM Trans. Graph. 24(2005), no. 3, 1062–1070.
[15] Taehyun Rhee, J.P. Lewis, Ulrich Neumann, and Krishna Nayak, Scan-based volume animation driven
by locally adaptive articulated registrations, IEEE
Trans. Visualization and Computer Graphics 17(2011), no. 3.
[16] B. Schölkopf and A. J. Smola, Learning with kernels,
MIT Press, 2002.
[17] H. Todo, K. Anjyo, W. Baxter, and T. Igarashi, Locally controllable stylized shading, ACM Trans. Graph.
26 (2007), no. 3, 17.
[18] Grace Wahba, Spline models for observational data,
SIAM, 1990.
[19] Holger Wendland,
Cambridge, 2005.

Scattered data approximation,[20] A. Yaglom, An introduction to the theory of stationary
random functions, Dover, 1973.
[21] K. Yosida, Functional analysis (3rd ed.), Springer,
1971.
Ken Anjyo
OLM Digital, Inc.
Dens Kono Bldg. 302 ViPlus, Tokyo, Japan
E-mail: anjyo(at)acm.org
J. P. Lewis
Weta Digital and Victoria University, Wellington, New ZealandE-mail: zilla(at)computer.org