Interpolation and Splines

1
Lecture 9

• Interpolation and Splines

2
Lingo

• Interpolation filling in gaps in data
• Find a function f(x) that
• 1) goes through all your data points
• 2) does something sensible in between

3
Lingo

• Splines a broad class of ways of performing
  interpolation
• (well get to the details, eventually)

4

• Find a function f(x) that
• 1) goes through all your data points
• (observations)
• 2) does something sensible in between
• (prior information)
• Why not just use least-squares?

5

• Remember this?
• mest mA M dobs GmA
• where M GTCd-1G Cm-1-1 GT Cd-1

6

• m a vector of all the points at which you want
  to estimate the function, including the points
  for which you have observations
• d a vector of just those points where you have
  observations
• So the equation Gmd is very simple, a model
  parameter equals the data when the corresponding
  observation is available

0 0 1 0 0
mi
dj

Just a single 1 per row
7

• You then implement a smoothness constraint
  through minimizing Dm2, where D is some measure
  of the non-smoothness of m
• Thus mA 0 and Cm-1 e2DTD and CdI

0 1 -2 1 0
D
One possibility is to use the finite-difference
approximation of the second derivative
8

• First derivative
• dm/dxi ? (1/Dx) mi mi-1
• ? mi mi-1
• Second derivative
• d2m/dx2i ? dm/dxi1 - dm/dxi
• mi1 mi mi mi-1
• mi1 2mi mi-1

9
example

• 101 equally spaced along the x-axis
• So 101 values of the function f(x)
• 40 of these values measured (the data, d)
• the rest are unknown
• Two prior information
• minimize 2nd derivative for interior 99 xs
• minimize 1st derivative at left and right xs
• (nice to have the same number
• of priors as unknowns, but not
• required)

10
e 10-6
data
f(x)
result
x
11
e can be chosen by trail and errorbut usually
the result fairly insensitive to e, as long as
its small
12
e varying over six orders of magnitude
f(x)
log10(Total Error)
log10(e)
x
13
A purist might say that this is not really
interpolation, because the curve goes through the
data only in the limit e?0but for small esthe
error is extremely small
14
an aside

• Construct an equation F m h as follows

d 0
G eD
m
then note FTF-1FT h GTGe2 DTD-1GT d so if
you want, you can just append eD to the bottom of
G and solve by simple least-squares
15
solved via solved viaFTF-1FT
h GTGe2 DTD-1GT d
exactly the same!
16
another reason to work with
d 0
G eD

• F m h

both G and D, and therefore F, too, are mostly
zero (that is, theyre sparse matrices) very
efficient algorithms are available for solving
Fmh in the least-squares sense when F is a
sparse matrix (note GTG and DTD are not as sparse
as G or D and GTG and DTD-1 is not sparse at
all)
17
2D Example(here a sparse solver would really be
useful, for the number of unknowns is very large)
18
44 observed data
21 unknowns
21?21441 unknowns
21 unknowns
19

• Prior information
• ?2f d2f/dx2 d2f/dy2 0 in interior of the
  box
• n??f 0 on edges of box
• a generalization of the 1D case

20
results
21
comparison
22
one limitation of this method is that it is
discreteit only gives the unknown function at
specific, prescribed values of xione might
prefer to have an analytic formula for the value
of the function at any x
23
LINGOan analytic formula that gives the value
of the function at any x is called an interpolant
24
high order polynomialsomething that sound like
a good idea but isnteven though an N-1
polynomal can computed to pass through any N
points
25
example 10th order polynomial fit to 11 points
Big swings not what we hoped for
26
solution

• simple function
• e.g. a low order polynomial
• that is valid in some interval near xi
• obviously, we need many such polynomials to over
  the whole x-axis
• This approach is called a spline

27
weve all used one already -linear splines
in this interval y(x) yi (yi1-yi)?(x-xi)/(xi
1-xi)
y
yi1
yi
x
xi
xi1
28
cubic splines somewhat more complicatedbut a
lot nicer
y
cubic abxcx2dx3 in this interval
yi1
a different cubic in this interval
yi
x
xi
xi1
29
counting up unknowns
unknowns N data N-1 intervals 4 coefficients per
interval 4(N-1)4N-4 coefficients total4N-4
unknowns
four coefficients a, b, c, d in every interval
y
constraints curve goes thru point at end of its
interval 2(N-1)2N-2 dy/dx match at interior
points N-2 constraints d2y/dx2 match at interior
points N-2 constraints d2y/dx2 0 at end
points 2 constraints total 4N-4 constraints
yi1
yi
x
xi
xi1
30
formulating the cubic spline problem in an
efficient manner

• f(x) with N observations (xi, fi)
• let hi Dxi xi1- xi
• and
• Dfi fi1-fi
• Si(x) are cubic polynomials, one for each
  interval

31

• Let the 2nd derivatives have values
• d2Si/dx2yi at the left hand end of its
  interval
• But since the second derivative is presumed
  continuous across intervals, d2Si/dx2yi1 on
  the right hand side of its interval too.
• since Si(x) is a cubic, its second derivative is
  a linear function
• So within an interval d2Si/dx2 varies linearly
• d2Si/dx2 yi (xj1-x)/hj yi1 (x-xj)/hj

well wind up solving for these yis
32
now integrate twice to get Si(x) d2Si/dx2 yi
(xj1-x)/hj yi1 (x-xj)/hj Si(x) yi
(xj1-x)3/(6hj) yi1 (x-xj)3/(6hj)
ci(x-xi) di(xi1-x) where ci and di are
integration constants
33
now choose the integration constants ci and di
such that the the cubic goes through the data
points. That is, Si(xi)fi and Si(xi1)fi1
this leads to ci fi1/hi yi1hi/6 di
fi/hi yihi/6
34
so Si(x) yi (xj1-x)3/(6hj) yi1
(x-xj)3/(6hj) fj1/hi yi1hi/6(x-xi)
fi/hi yihi/6(xi1-x) where ci and di are
integration constants
but we still havent implemented the continuity
of dS/dx condition
35
so we compute the derivative Si(x) dSi/dx
yi (xi1-x)2/(2hi) yi1
(x-xi)2/(2hi) fi1/hi yi1hi/6 -
fi/hi yihi/6 ½yi (xi1-x)2hi
½yi1 (x-xi)2/hi Dfi1/hi (yi1-yi)hi/6
36
now require the first derivative match across
intervals Si-1(xi)Si(xi)this leads to an
equation for the unknown yihi-1yi-1
2(hihi-1)yi hiyi1 bjwith bj 6Dfi/hi
6Dfi-1/hi-1
37
the equation for the unknown yihi-1yi-1
2(hihi-1)yi hiyi1 biis just a matrix
equation
i1 h0y0 2(h1h0)y1
h1y2 b1 i2 h1y1 2(h2h1)y2
h2y2 b2 i2 h2y2
2(h3h2)y3 h3y3 b3 iN-1
hN-2yN-2 2(hN-1hN-2)yN-1 hN-1yN
bN-1 iN hN-1yN-1 2(hNhN-1)yN
hNyN1 bN
well discuss the issue raise by y0 and yN1
in a moment
38
the matrix equation, with gi2(hihi-1), is
b0b1b2b3bNbN1
h0 g1 h1 h1 g2 h2 h2
g3 h3
hN-2 gN-1 hN-1
hN-1 gN hN
y0y1y2y3yNyN1

A Tridiagonal Matrix, by the way. Very fast
solvers are available
39
the matrix equation, with gi2(hihi-1), is
b0b1b2b3bNbN1
N2
h0 g1 h1 h1 g2 h2 h2
g3 h3
hN-2 gN-1 hN-1
hN-1 gN hN
y0y1y2y3yNyN1
N

Ive written this as if there were two extra
points, one to the left of the first point and
one to the right of the last point. Of course,
there arent. The way to handle this is to
prescribe y0 and yN1 and move them to the
r.h.s. of the equation.
40
moving over these two now-specified unknowns
b1 h0y0b2b3bN hN1yN1
N
g1 h1 h1 g2 h2 h2
g3 h3
hN-1 gN hN hN
gN
y1y2y3yN
N

we can set y0 and yN to whatever we want. A
simple choice is zero, in which case the splines
are called natural cubic splines
41
example of cubic spline interpolation
42
very easy in MatLab

• new_y spline(x,y,new_x)