An Introduction to Scattered Data Approximation

Greg Coombe
October 31, 2006

1 Scattered Data Approximation

Scattered Data Approximation is the name given to a set of techniques which attempt to fill a domain with
reasonable data values given samples of the data in arbitrary locations. More formally, given a set of sites
V = {ν0 , . . . , νm } ∈ Rd and function values at these sites fν ∈ R, we seek an approximant f : Rd → R which
as close as possible our data.
In general, there are two types of functions; interpolating and approximating. They are differentiated by
their behavior at the sites ν. The residual at the sites is

 = kf (ν) − fν k for ν ∈ V

For interpolating functions,  = 0 for all ν ∈ V. For approximating functions, the residual is minimized for
the given set of basis functions.

1.1 Piecewise Polynomial

The most common scattered data approximation technique is to triangulate the data sites. This defines the
regions where the polynomials lie as well as how they are joined at common vertices and edges. It is usually
used for two-dimensional data, as higher dimensions are too difficult to triangulate. The quality of the spline
approximation depends strongly on triangulation, so it is very important to get a good triangulation; long
thin triangles will destroy the accuracy of the approximation.
The most commonly used triangulation is the Delaunay triangulation (). The Delaunay triangulation is
the dual of the Voronoi diagram. The Voronoi diagram is a partitioning of the space into cells such that
every point within the cell is closer to one site than any other site. This is written as:

Ty = {x ∈ R2 s.t. kx − yk = min kx − νk, y and ν ∈ V}

ν

The Delaunay triangulation is constructed from the Voronoi diagram by using the sites as vertices and
connecting neighboring sites with edges.
A piece-wise linear interpolation over Delaunay triangulation can is defined by linearly interpolating the
function values from the vertices. This can be extended to a piece-wise quadratic function by using additional
information such as the mid-points of the edges. In general, extending this to higher-order continuity requires
subdividing the triangulation.

1.2 Shepard’s Method

Shepard’s Method is a non-polynomial global multivariate interpolation scheme for scattered data. The
advantage of global methods is that they define a single function f (x) over all of Rd , instead of a (potentially
large) set of local piece-wise polynomial functions. Shepard’s Method has the form:
P
fν φ(kx − νk)
f (x) = Pν∈V
ν∈V φ(kx − νk)

1
The weight functions φ(r) are chosen to be large at the sites and to fall off away from the points. A common
weight function is φ(r) = r−µ . These radially-symmetric functions decay as the distance from the site
increases, where the µ parameter controls the rate of decay. One problem with these functions is that they
are not compact, so their support can extend globally. Consequently, another popular choice for weight
function is the bounded exponential decay function:
exp[−t2 /(t2 −r)2 ]
(
exp[r 2 /h2 ]−1 if r ≤ t
φ(r) =
0 otherwise
The user parameters consist of the support radius t and the decay rate h.
Shepard’s Method is based on the Partition of Unity (?), which is a general technique to weight a set of
basis functions such that they sum to 1. The form of the Partition of Unity is:
P
b(x)φ(x)
f (x) = P
φ(x)

The Partition of Unity is a non-linear combination of the basis functions, which we will talk about further
in Section 2.2. The disadvantage of the Partition of Unity is that it does not preserve the differentiability
or continuity of the basis functions. For example, a major disadvantage of Shepard’s Method is that it has
stationary points (vanishing gradients) at all data sites.

1.3 Radial Basis Functions

Another common scattered data approximation technique is the Radial Basis Function(RBF), which consists
of a finite linear combination of translated basis functions. These basis functions are radially-symmetric
functions of the form φ(k · k), where k · k is the Euclidean norm. The general form of f is:
X
f (x) = λk φ(kx − ck k)
k

where λk ∈ R are the coefficients of the basis functions and ck ∈ Rd are known as the centers. For
interpolating functions, the centers are defined to be the original sites V. For approximating functions, the
centers are optimized to minimize the residual.
One advantage of RBFs is that they are independent of dimension. The basis functions φ(r) take Eulidean
distance as input, so they can be trivially extended to arbitrary dimensions. This is in contrast to piece-wise
polynomial schemes, which usually rely on the tensor product of the dimensions.
The basis functions are chosen to emphasize certain desirable qualities such as continuity, hole-filling, or
differentiability. Commonly used basis functions include:

φ(r) = r bi-harmonic
φ(r) =√r2 log r thin-plate splines
φ(r) = r2 + c2 multiquadrics
2
φ(r) = e−αr Gaussian
Once the basis functions are chosen, the coefficients can be computed by solving the linear system:

Aλ = f

where A = {φ(kxi − xj k}, and xi and xj are the centers. A is known as the interpolation matrix and is
square, symmetric, and non-singular (?). This is one of the primary advantages of RBFs; unlike high-degree
spline approximation or polynomial interpolation, RBF interpolation is always uniquely solvable.
Note that the biharmonic, thin-plate splines, and the multiquadrics all increase as the distance from
the point increases, while the Gaussian falls off to zero. This behavior effects the form of the interpolation
matrix. For the biharmonic function, the matrix will be very dense with a zero diagonal. Gaussian basis

2
functions, which have local support, generate a matrix which is banded and sparse. This local support also
means that adding a point only effects a small set of interpolation coefficients. Adding a point into a global
RBF requires recomputing the entire set of interpolation coefficients.
The combination of local support and matrix form would seem to imply that Gaussians are a much
better choice for basis function. However, it order to effectively use a local function such as a Gaussian, the
sampling density must be fairly uniform, which is not often the case with many datasets. Global functions
such as thin-plate splines and the multiquadrics are very good at filling holes in the sampling domain. These
occur, for example, when reconstructing surfaces from laser-scanned data. Thus global basis functions tend
to be more widely used, and research is focused on fast methods for solving the interpolation matrix (?).

2 Least Squares
Least Squares methods are a set of linear approximation techniques for scattered data. As before, we are
given a set of N scalar samples fν ∈ R at points ν ∈ Rd , and we want a globally-defined function f (x) that
best approximates the fν samples. For the approximation problem, the goal is to generate this function f (x)
such that the distance between the scalar data values fν and the function evaluated at the data points f (νi )
is as small as possible. This is written as:
X
min{ kf (νi ) − fν k}
i
(Note: this discussion follows Nealen (?)). Typically, f (x) is a polynomial of degree m in d spatial dimensions.
Thus f (x) can be written as

f (x) = Bc
T
where c = [c1 . . . ck ] is the unknown coefficient vector. The B matrix is the polynomial basis matrix, and
is composed of a set of basis functions b(~x). These functions are chosen based on the properties of the data
and the dimensionality. In our case, we use either the 2D quadratic basis set

b(x, y) = c1 + c2 x + c3 y + c4 x2 + c5 xy + c6 y 2
or the 2D linear basis set

b(x, y) = c1 + c2 x + c3 y
The polynomial basis matrix B is creating by stacking all of the k polynomial basis equations:
   
b1 (x, y) 1 x1 y1
 b2 (x, y)   1 x2 y2 
B= = 
 ···   ··· 
bk (x, y) 1 xk yk
The data points are also written as a vector f = [f1 (x) . . . fk (x)]T . The fitting problem can now be stated
as finding the coefficient vector c that minimizes the error between the polynomial approximation Bc and
the data values f :
min{||Bc − f k}
c
This can be solved using the Method of Normal Equations (?):
Bc = f
B Bc = B T f
T

c = [B T B]−1 B T f
Solving the Normal Equation requires inverting the matrix B T B. Since this inversion must be performed
for every Least Squares fit, it is worthwhile to look into how this can be performed efficiently.

3
2.1 Solving the Normal Equations
The matrix inversion is possible because B T B is square and non-singular (?). The size of the matrix depends
upon the dimensionality d of the data and the degree k of the polynomial basis. For small matrices, it can be
inverted directly. For larger matrices, I will discuss several efficient algorithms, including QR Decomposition
and SVD. These are often implemented in matrix inversion packages such as BLAS (?) and TGT (?).

2.1.1 QR Decomposition
One problem with the Normal Equations is that the matrix B T B can be ill-conditioned, and thus subject
to numerical precision errors when inverting. This can be addressed by using a technique called QR De-
composition. This technique decomposes a matrix B into an upper-trapezoidal matrix R and an orthogonal
matrix Q such that B = QR. By substituting this into the above Normal Equations and using the fact that
QT = Q−1 :

c = [B T B]−1 B T f
= [(QR)T QR]−1 (QR)T f
= [RT QT QR]−1 RT QT f
= R−1 (RT )−1 RT QT f
= R−1 QT f

This formulation is much less prone to numerical precision errors. The QR Decomposition can be com-
puted using three different techniques; Gram-Schmidt orthogonalization, Householder transformations, or
Givens rotations. Householder transformations are most commonly used, but there are advantages and dis-
advantages to each of these. For more information on these techniques, see (?)[note: or any textbook on
Linear Algebra! ].

2.1.2 Singular Value Decomposition

Another technique is to invert the matrix using the SVD . Recall from Section ?? that the SVD decomposes
a matrix A = U SV T . U and V are orthogonal and S is a diagonal matrix with the elements σ1 , σ2 , . . . , σn
along the diagonal. This leads to a simple method for computing A−1 :

A−1 = [U SV T ]−1
= [V T ]−1 S −1 U −1
= V S −1 U T

where S −1 is a diagonal matrix with the elements 1/σ1 , 1/σ2 , . . . , 1/σn along the diagonal. When A is a
square matrix, this computes the inverse of the matrix. When A is not square a true matrix inverse does
not exist, and this is known as a pseudo-inverse. The pseudo-inverse A+ is a generalization of the inverse
which satisfies the linear equation x = A+ b.

2.2 Weighted Least Squares

One of the limitations of Least Squares fitting is that the solution encompasses the entire domain. This
global complexity makes it difficult to handle large data sets or data sets with local high frequencies. We
prefer a method that considers samples that are nearby as more important than samples that are far away.
This can be accomplished by adding a distance-weighting term φ(r) to the Least Squares minimization. We
are now trying to minimize the function

4
Figure 1: A diagram of the weighted least squares approach to function representation. Each center con-
structs a low-degree polynomial approximation based on samples in their neighborhood. These local approx-
imations are then combined to form a global approximation.

X
min{ φ(kx − xi k)kf (xi ) − fi k}
i

Common choices for the distance-weighting basis functions φ(r) are the Wendland function (?)
4r r
φ(r) = (1 + )(1 − )4
h h
and the interpolation function
φ(r) = 1 + 2(r/h)3 − 3(r/h)2
which are both 1.0 at d = 0, and falls off to zero at the edges of the support radius h. These functions
are illustrated in Figure 2. Both functions have compact support, so each sample only affects a small
neighborhood.
It is interesting to note the connection with the Radial Basis Functions discussed in Section 1.3. In fact,
we could use some of the other weighting functions that were discussed such as multiquadrics and thin-plate
splines. However, these functions have infinite support and thus must be solved globally. The advantage
of Weighted Least Squares is that instead of evaluating a single global approximation for all of the data
samples, we create a set of local approximations. We then blend these local approximations together to get
a global approximation. This provides us with an local approximation behavior through the Least Squares
fitting, which is useful because the data that we get is often noisy.
Each of the approximations is associated with a point x̄, which is known as thecenters (the concept
is analogous to the centers of Radial Basis Functions). At each of these centers, a low-degree polynomial
approximation is computed using the distance-weighted samples xi in the local neighborhood. This means
that the coefficients are now a function of x̄, and only defined locally around these centers. We can define a
distance weighting vector Φx̄ with elements φi = φ(kx̄ − xi k).
The system is solved using Normal Equations:

cx̄ = [Φx̄ B T B]−1 Φx̄ B T f

5
 Distance Weighting Functions with Support = 0.5
1.6
Interpolation
Wendland
1.4

1.2

1
Weight

0.8

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Distance

Figure 2: Two possible distance-weighting functions, with a variable support. Both functions are 1 at the
center, and taper off to zero as the distance increases. While they are well-behaved within the support
interval, they must be clamped outside this interval to avoid spurious weighting.

6
 Expanding this out in equation form:
X X
c(x̄) = [ φ(kx̄ − xi k)b(xi )b(xi )T ]−1 φ(kx̄ − xi k)b(xi )fi
i i

We now have a set of local approximations at each center. During the reconstruction step, we need to
combine these local approximations to form a global approximation. Since this global function is a weighted
combination of the basis functions, it has the same continuity properties.
The first step is to determine the m nearby local approximations that overlap this point and combine
them using a weight based on distance. However, the functions cannot just be added together, since the
weights may not sum to 1.0. To get the proper weighting of the local approximations, we use a technique
known as the Partition of Unity (?), which allows us to extend the local approximations to cover the entire
domain. A new set of weights Θj are computed by considering all of the m local approximations that overlap
this point

φj (r)
Θj (r) = Pm
k=1 φk (r)
The global approximation of this function is computed by summing the weighted local approximations:
m
X
f (x) = Θj (x)b(x)T c(x̄j )
j=1

This representation allows us to place the centers x̄ at the most effective locations. We discuss several
strategies for center placement in Section ??.
Note that Weighted Least Squares is very similar to Shepard’s Method, which we discussed in Section 1.2.
Recall that Shepard’s Method has the form:
P
fν φ(kx − νk)
f (x) = P
φ(kx − νk)
There are several differences between the two methods. Weighted Least Squares uses a distance-weighted
Least Squares approximation to the points, while Shepard’s Method uses the distance-weighted points di-
rectly. [note: What advantage does this give us? Since Shepard’s Method has vanishing points at sites, does
WLS? ]

2.3 Moving Least Squares

Moving Least Squares (?; ?) is another technique for approximating scattered data using Least Squares.
Instead of establishing a set of centers and computing Weighted Least Squares approximations for these
centers, the Moving Least Squares technique treats each evaluation point individually as an center. The
center is “moved” over the parameter domain and a new approximation is computed everywhere. This
generates a continuously differentiable global function f (x) if and only if the weighting function φ(r) is
continuously differentiable (?). Moving Least Squares consists of a set of local functions fx (x):
X
f (x) = fx (x), min{ φ(kx − xi k)kfx (xi ) − fi k}
fx
i

The functions fx (x) are constructed and evaluated at each point in the domain. The parameters of the basis
functions can be adjusted to smooth the data, which makes it useful for reconstructing functions from noisy
data (?).