Lecture AIMS Notes 2006
Peter J. Olver
13. Approximation and Interpolation
We will now apply our minimization results to the interpolation and least squares fitting of data and functions.
13.1. Least Squares.
Linear systems with more equations than unknowns typically do not have solutions. In such situations, the least squares solution to a linear system is one means of getting as close as one can to an actual solution. Definition 13.1. A least squares solution to a linear system of equations Ax = b is a vector x R n that minimizes the Euclidean norm A x - b . If the system (13.1) actually has a solution, then it is automatically the least squares solution. Thus, the concept of least squares solution is new only when the system does not have a solution. To find the least squares solution, we need to minimize the quadratic function Ax - b where
2
(13.1)
= xT AT A x - 2 xT AT b + bT b = xT Kx - 2 xT f + c, f = AT b, c= b
2
= (A x - b)T (A x - b) = (xT AT - bT )(A x - b)
K = AT A,
.
(13.2)
According to Theorem 12.10, the Gram matrix K = AT A is positive definite if and only if ker A = {0}. In this case, Theorem 12.12 supplies us with the solution to this minimization problem. Theorem 13.2. Assume that ker A = {0}. Set K = AT A and f = AT b. Then the least squares solution to the linear system A x = b is the unique solution x to the so-called normal equations Kx = f namely x = (AT A)-1 AT b. 3/15/06 212
c 2006
or, explicitly,
(AT A) x = AT b,
(13.3) (13.4)
Peter J. Olver
The least squares error is Ax - b
2
= b
2
- fT x = b
2
- bT A (AT A)-1 AT b.
(13.5)
Note that the normal equations (13.3) can be simply obtained by multiplying the original system A x = b on both sides by AT . In particular, if A is square and invertible, then (AT A)-1 = A-1 (AT )-1 , and so (13.4) reduces to x = A-1 b, while the two terms in the error formula (13.5) cancel out, producing zero error. In the rectangular case -- when inversion is not allowed -- (13.4) gives a new formula for the solution to a compatible linear system A x = b. Example 13.3. Consider the linear system x1 + 2 x 2 = 1, 3 x1 - x2 + x3 = 0,
- x1 + 2 x2 + x3 = -1, x1 - x2 - 2 x3 = 2, 2 x1 + x2 - x3 = 2, consisting of 5 equations in 3 unknowns. The coefficient matrix and right hand side are 1 3 A = -1 1 2 2 0 -1 1 2 1 , -1 -2 1 -1 1 0 b = -1 . 2 2
A direct application of Gaussian Elimination shows that the system is incompatible -- it has no solution. Of course, to apply the least squares method, we are not required to check this in advance. If the system has a solution, it is the least squares solution too, and the least squares method will find it. To form the normal equations (13.3), we compute 16 -2 -2 8 K = AT A = -2 11 2 , f = AT b = 0 . -2 2 7 -7 Solving the 3 3 system K x = f by Gaussian Elimination, we find x = K -1 f ( .4119, .2482, -.9532 )
T
to be the least squares solution to the system. The least squares error is b - Ax ( -.0917, .0342, .1313, .0701, .0252 )
T
.1799,
which is reasonably small -- indicating that the system is, roughly speaking, not too incompatible. 3/15/06 213
c 2006 Peter J. Olver
13.2. Data Fitting and Interpolation.
One of the most important applications of the least squares minimization process is to the fitting of data points. Suppose we are running an experiment in which we measure a certain time-dependent physical quantity. At time ti we make the measurement yi , and thereby obtain a set of, say, m data points (t1 , y1 ), (t2 , y2 ), ... (tm , ym ). (13.6)
Suppose our theory indicates that the data points are supposed to all lie on a single line y = + t, (13.7)
whose precise form -- meaning its coefficients , -- is to be determined. For example, a police car is interested in clocking the speed of a vehicle by using measurements of its relative distance at several times. Assuming that the vehicle is traveling at constant speed, its position at time t will have the linear form (13.7), with , the velocity, and , the initial position, to be determined. Experimental error will almost inevitably make this impossible to achieve exactly, and so the problem is to find the straight line (13.7) that "best fits" the measured data and then use its slope to estimate the vehicle's velocity. The error between the measured value yi and the sample value predicted by the function (13.7) at t = ti is ei = yi - ( + ti ), i = 1, . . . , m.
We can write this system of equations in matrix form as e = y - A x, where e1 e2 . . . em y1 y2 . . . ym 1 , A=. . . 1 tm 1 t1 t2 . . .
We call e the error vector and y the data vector . The coefficients , of our desired function (13.7) are the unknowns, forming the entries of the column vector x. If we could fit the data exactly, so yi = + ti for all i, then each ei = 0, and we could solve A x = y for the coefficients , . In the language of linear algebra, the data points all lie on a straight line if and only if y rng A. If the data points are not collinear, then we seek the straight line that minimizes the total squared error or Euclidean norm Error = e = e2 + + e 2 . 1 m
e=
,
y=
,
while
x=
.
(13.8)
Pictorially, referring to Figure 13.1, the errors are the vertical distances from the points to the line, and we are seeking to minimize the square root of the sum of the squares of 3/15/06 214
c 2006 Peter J. Olver
Figure 13.1.
Least Squares Approximation of Data by a Straight Line.
the individual errors , hence the term least squares. In other words, we are looking for the T coefficient vector x = ( , ) that minimizes the Euclidean norm of the error vector e = Ax - y . (13.9)
Thus, we recover the problem of characterizing the least squares solution to the linear system A x = y. Theorem 13.2 prescribes the solution to this least squares minimization problem. We form the normal equations (AT A) x = AT y, with solution x = (AT A)-1 AT y. (13.10)
Invertibility of the Gram matrix K = AT A relies on the assumption that the matrix A has linearly independent columns. For the particular matrix in (13.8), linear independence of its two columns requires that not all the ti 's be equal, i.e., we must measure the data at at least two distinct times. Note that this restriction does not preclude measuring some of the data at the same time, e.g., by repeating the experiment. However, choosing all the ti 's to be the same is a silly data fitting problem. (Why?)
This choice of minimization may strike the reader as a little odd. Why not just minimize the sum of the absolute value of the errors, i.e., the 1 norm e 1 = | e1 | + + | en | of the error vector, or minimize the maximal error, i.e., the norm e = max{ | e1 |, . . . , | en | }? Or, even better, why minimize the vertical distance to the line? The perpendicular distance from each data point to the line might strike you as a better measure of error. The answer is that, although each of these alternative minimization criteria is interesting and potentially useful, they all lead to nonlinear minimization problems, and so are much harder to solve! The least squares minimization problem can be solved by linear algebra, and so, purely on the grounds of simplicity, is the method of choice in most applications. Moreover, one needs to fully understand the linear problem before diving into more treacherous nonlinear waters.
3/15/06
215
c 2006
Peter J. Olver
where the overbars, namely 1 t= m
m
Under this assumption, we then compute 1 t 1 1 t2 m ti 1 t 1 1 ... 1 AT A = =m . = 2 . . t1 t2 . . . t m . t t2 ti (ti ) . . 1 tm y 1 y2 y yi 1 1 ... 1 = AT y = , =m . t1 t2 . . . t m . ty t i yi . ym ti ,
i=1
, (13.11)
1 y= m
m
yi ,
i=1
t2
1 = m
m
ti2 ,
i=1
1 ty = m
m
t i yi ,
i=1
(13.12)
denote the average sample values of the indicated variables. Warning: The average of a product is not equal to the product of the averages! In particular, t2 = ( t ) 2 , t y = t y. Substituting (13.11) into the normal equations (13.10), and canceling the common factor of m, we find that we have only to solve a pair of linear equations + t = y, for the coefficients: = y - t , = t2 ty-ty - ( t )2 = (ti - t ) yi . (ti - t )2 (13.13) t + t2 = t y,
Therefore, the best (in the least squares sense) straight line that fits the given data is y = (t - t ) + y, where the line's slope is given in (13.13). Example 13.4. Suppose the data points are given by the table ti yi 0 2 1 3 3 7 6 12 (13.14)
To find the least squares line, we construct 1 0 1 1 1 A= AT = , 1 3 0 1 6 3/15/06 216
1 1 1 3
1 6
,
2 3 y = . 7 12
c 2006 Peter J. Olver
12 10 8 6 4 2 1 2 3 4 5 6
Figure 13.2. Therefore AT A =
Least Squares Line.
4 10 10 46
,
AT y =
24 96
.
The normal equations (13.10) reduce to 4 + 10 = 24, 10 + 46 = 96, so =
12 7 ,
=
12 7 .
Therefore, the best least squares fit to the data is the straight line y=
12 7
+
12 7
t.
Alternatively, one can compute this formula directly from (13.1314). As you can see in Figure 13.2, the least squares line does a fairly good job of approximating the data points. Example 13.5. Suppose we are given a sample of an unknown radioactive isotope. At time ti we measure, using a Geiger counter, the amount mi of radioactive material in the sample. The problem is to determine the initial amount of material along with the isotope's half life. If the measurements were exact, we would have m(t) = m 0 e t , where m0 = m(0) is the initial mass, and < 0 the decay rate. The half-life is given by t = -1 log 2. As it stands this is not a linear least squares problem. But it can be easily converted to the proper form by taking logarithms: y(t) = log m(t) = log m0 + t = + t. We can thus do a linear least squares fit on the logarithms yi = log mi of the radioactive mass data at the measurement times ti to determine the best values for and = log m0 . Polynomial Approximation and Interpolation The basic least squares philosophy has a variety of different extensions, all interesting and all useful. First, we can replace the straight line (13.7) by a parabola defined by a quadratic function y = + t + t2 . (13.15) 3/15/06 217
c 2006 Peter J. Olver
For example, Newton's theory of gravitation says that (in the absence of air resistance) a falling object obeys the parabolic law (13.15), where = h 0 is the initial height, = v0 is the initial velocity, and = - 1 g is minus one half the gravitational constant. Suppose 2 we observe a falling body on a new planet, and measure its height y i at times ti . Then we can approximate its initial height, initial velocity and gravitational acceleration by finding the parabola (13.15) that best fits the data. Again, we characterize the least squares fit by minimizing the sum of the squares of the individual errors ei = yi - y(ti ). The method can evidently be extended to a completely general polynomial function y(t) = 0 + 1 t + + n tn
m
(13.16)
of degree n. The total least squares error between the data and the sample values of the function is equal to e where 1 1 A=. . . t1 t2 . . . t2 1 t2 2 . . . t2 m ... ... .. . ... tn 1 tn 2 . . . tn m 1 x = 2 , . . . n
0 2
=
i=1
yi - y(ti )
2
= y - Ax
2
,
(13.17)
1 tm
The coefficient m (n + 1) coefficient matrix is known as a Vandermonde matrix , named after the eighteenth century French mathematician, scientist and musicologist Alexandre Thophile Vandermonde -- despite the fact that it appears nowhere in his four mathemate ical papers! In particular, if m = n + 1, then A is square, and so, assuming invertibility, we can solve A x = y exactly. In other words, there is no error, and the solution is an interpolating polynomial , meaning that it fits the data exactly. A proof of the following result can be found at the end of this section. Lemma 13.6. If t1 , . . . , tn+1 are distinct, ti = tj , then the (n + 1) (n + 1) Vandermonde interpolation matrix (13.18) is nonsingular. This result immediately implies the basic existence theorem for interpolating polynomials. Theorem 13.7. Let t1 , . . . , tn+1 be distinct sample points. Then, for any prescribed data y1 , . . . , yn+1 , there exists a unique interpolating polynomial of degree n with the prescribed sample values y(ti ) = yi for all i = 1, . . . , n + 1. Thus, two points will determine a unique interpolating line, three points a unique interpolating parabola, four points an interpolating cubic, and so on; see Figure 13.3. The basic ideas of interpolation and least squares fitting of data can be applied to approximate complicated mathematical functions by much simpler polynomials. Such approximation schemes are used in all numerical computations. Your computer or calculator is only able to add, subtract, multiply and divide. Thus, when you ask it to compute t or 3/15/06 218
c 2006 Peter J. Olver
,
y=
y1 y2 . . . ym
.
(13.18)
Linear Figure 13.3.
Quadratic Interpolating Polynomials.
Cubic
et or cos t or any other non-rational function, the program must rely on an approximation scheme based on polynomials . In the "dark ages" before computers, one would consult precomputed tables of values of the function at particular data points. If one needed a value at a nontabulated point, then some form of polynomial interpolation would be used to accurately approximate the intermediate value. Example 13.8. Suppose that we would like to compute reasonably accurate values for the exponential function et for values of t lying in the interval 0 t 1 by approximating it by a quadratic polynomial p(t) = + t + t2 . (13.19)
If we choose 3 points, say t1 = 0, t2 = .5, t3 = 1, then there is a unique quadratic polynomial (13.19) that interpolates et at the data points, i.e., p(ti ) = eti for i = 1, 2, 3.
In this case, the coefficient matrix (13.18), namely 1 0 0 A = 1 .5 .25 , 1 1 1
is nonsingular. Therefore, we can exactly solve the interpolation equations 1. e t1 A x = y, where y = et2 = 1.64872 e t3 2.71828 is the data vector, which we assume we already know. The solution 1. x = = .876603 .841679
Actually, one could also allow interpolation and approximation by rational functions, a subject known as Pad approximation theory, [ 3 ]. e
3/15/06
219
c 2006
Peter J. Olver
3
3
3
2.5
2.5
2.5
2
2
2
1.5
1.5
1.5
1
1
1
0.5
0.5
0.5
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
Figure 13.4.
Quadratic Interpolating Polynomial for et .
yields the interpolating polynomial p(t) = 1 + .876603 t + .841679 t2 . (13.20)
It is the unique quadratic polynomial that agrees with et at the three specified data points. See Figure 13.4 for a comparison of the graphs; the first graph shows e t , the second p(t), and the third lays the two graphs on top of each other. Even with such a primitive interpolation scheme, the two functions are quite close. The maximum error or L norm of the difference is et - p(t)
= max | et - p(t) |
0t1
.01442,
with the largest deviation occurring at t .796. There is, in fact, an explicit formula for the interpolating polynomial that is named after the influential eighteenth century ItaloFrench mathematician JosephLouis Lagrange. Suppose we know the solutions x1 , . . . , xn+1 to the particular interpolation systems A xk = ek , k = 1, . . . , n + 1, (13.21)
where e1 , . . . , en+1 are the standard basis vectors of R n+1 . Then the solution to A x = y = y1 e1 + + yn+1 en+1 is given by the superposition formula x = y1 x1 + + yn+1 xn+1 . The particular interpolation equation (13.21) corresponds to the interpolation data y = e k , meaning that yk = 1, while yi = 0 at all points ti with i = k. If we can find the n + 1 particular interpolating polynomials that realize this very special data, we can use superposition to construct the general interpolating polynomial. Theorem 13.9. Given distinct sample points t1 , . . . , tn+1 , the k th Lagrange interpolating polynomial is given by Lk (t) = (t - t1 ) (t - tk-1 )(t - tk+1 ) (t - tn+1 ) , (tk - t1 ) (tk - tk-1 )(tk - tk+1 ) (tk - tn+1 ) 220 k = 1, . . . , n + 1. (13.22)
c 2006 Peter J. Olver
3/15/06
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
L1 (t) Figure 13.5.
L2 (t)
L3 (t)
Lagrange Interpolating Polynomials for the Points 0, .5, 1.
It is the unique polynomial of degree n that satisfies Lk (ti ) = 1, 0, i = k, i = k, i, k = 1, . . . , n + 1. (13.23)
Proof : The uniqueness of the Lagrange interpolating polynomial is an immediate consequence of Theorem 13.7. To show that (13.22) is the correct formula, we note that when t = ti for any i = k, the factor (t - ti ) in the numerator of Lk (t) vanishes, while the denominator is not zero since the points are distinct. On the other hand, when t = t k , the numerator and denominator are equal, and so Lk (tk ) = 1. Q.E.D. Theorem 13.10. If t1 , . . . , tn+1 are distinct, then the polynomial of degree n that interpolates the associated data y1 , . . . , yn+1 is p(t) = y1 L1 (t) + + yn+1 Ln+1 (t). Proof : We merely compute p(tk ) = y1 L1 (tk ) + + yk Lk (t) + + yn+1 Ln+1 (tk ) = yk , where, according to (13.23), every summand except the k th is zero. Q.E.D. Example 13.11. For example, the three quadratic Lagrange interpolating polynomials for the values t1 = 0, t2 = 1 , t3 = 1 used to interpolate et in Example 13.8 are 2 L1 (t) = L2 (t) = L3 (t) = t- 0-
1 2 1 2 1 2
(13.24)
Thus, we can rewrite the quadratic interpolant (13.20) to et as y(t) = L1 (t) + e1/2 L2 (t) + e L3 (t) = (2 t2 - 3 t + 1) + 1.64872(- 4 t2 + 4 t) + 2.71828(2 t2 - t). 3/15/06 221
c 2006 Peter J. Olver
(t - 0) t - 1 2 (1 - 0) 1 - 1 2
(t - 0)(t - 1) -0 1 -1 2
(t - 1) = 2 t2 - 3 t + 1, (0 - 1) = - 4 t2 + 4 t, = 2 t2 - t. (13.25)
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
-3
-2
-1 -0.2
1
2
3
-3
-2
-1 -0.2
1
2
3
-3
-2
-1 -0.2
1
2
3
Figure 13.6.
Degree 2, 4 and 10 Interpolating Polynomials for 1/(1 + t 2 ).
We stress that this is the same interpolating polynomial -- we have merely rewritten it in the more transparent Lagrange form. You might expect that the higher the degree, the more accurate the interpolating polynomial. This expectation turns out, unfortunately, not to be uniformly valid. While low degree interpolating polynomials are usually reasonable approximants to functions, high degree interpolants are not only more expensive to compute, but can be rather badly behaved, particularly near the ends of the interval. For example, Figure 13.6 displays the degree 2, 4 and 10 interpolating polynomials for the function 1/(1 + t 2 ) on the interval -3 t 3 using equally spaced data points. Note the rather poor approximation of the function near the ends of the interval. Higher degree interpolants fare even worse, although the bad behavior becomes more and more concentrated near the endpoints. As a consequence, high degree polynomial interpolation tends not to be used in practical applications. Better alternatives rely on least squares approximants by low degree polynomials, to be described next, and interpolation by piecewise cubic splines, a topic that will be discussed in depth later. If we have m > n + 1 data points, then, usually, there is no degree n polynomial that fits all the data, and so we must switch over to a least squares approximation. first The requirement is that the associated m (n + 1) interpolation matrix (13.18) has rank n + 1; this follows from Lemma 13.6, provided that at least n + 1 of the values t 1 , . . . , tm are distinct. Thus, given data at m n + 1 different sample points t1 , . . . , tm , we can uniquely determine the best least squares polynomial of degree n that fits the data by solving the normal equations (13.10). Example 13.12. Let us return to the problem of approximating the exponential function et . If we use more than three data points, but still require a quadratic polynomial, then we can no longer interpolate exactly, and must devise a least squares approximant. For instance, using five equally spaced sample points t1 = 0, t2 = .25, t3 = .5, t4 = .75, t5 = 1, the coefficient matrix and sampled data vector (13.18) are 1 1 A = 1 1 1 3/15/06 0 0 .25 .0625 .5 .25 , .75 .5625 1 1 222 1. 1.28403 y = 1.64872 . 2.11700 2.71828
c 2006 Peter J. Olver
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
Figure 13.7.
Quadratic Approximant and Quartic Interpolant for et .
The solution to the normal equations (13.3), with 5. K = AT A = 2.5 1.875 is x = K -1 f = ( 1.00514, .864277, .843538 ) . This leads to the quadratic least squares approximant p2 (t) = 1.00514 + .864277 t + .843538 t2 . On the other hand, the quartic interpolating polynomial p4 (t) = 1 + .998803 t + .509787 t2 + .140276 t3 + .069416 t4 is found directly from the data values as above. The quadratic polynomial has a maximal error of .011 over the interval [ 0, 1 ] -- slightly better than the quadratic interpolant -- while the quartic has a significantly smaller maximal error: .0000527. (In this case, high degree interpolants are not ill behaved.) See Figure 13.7 for a comparison of the graphs. Proof of Lemma 13.6 : We will establish the rather striking L U factorization of the transposed Vandermonde matrix V = AT , which will immediately prove that, when t1 , . . . , tn+1 are distinct, both V and A are nonsingular matrices. The 4 4 case is instructive for understanding the general pattern. Applying regular Gaussian Elimination, 3/15/06 223
c 2006 Peter J. Olver T
2.5 1.875 1.5625
1.875 1.5625 , 1.38281
8.76803 f = AT y = 5.45140 , 4.40153
we find the explicit L U factorization 1 0 1 1 1 1 1 t1 t2 t3 t4 t1 2 2 2 2= 2 t1 t2 t3 t4 t1 t1 + t 2 t3 1 t3 2 t3 3 t3 4 t3 1
0 0 1
0
In the general (n + 1) (n + 1) case, the individual entries of the matrices appearing in factorization V = L U are vij = ti-1 , j
ij
t2 + t 1 t2 + t 2 t1 + t 2 + t 3 1 1 2 1 1 1 t3 - t 1 0 t2 - t1 0 0 (t3 - t1 ) (t3 - t2 ) 0 0 0
0 0
1 t4 - t 1 (t4 - t1 ) (t4 - t2 )
(t4 - t1 ) (t4 - t2 ) (t4 - t3 )
.
i, j = 1, . . . , n + 1,
ii ij
(13.26) = 1, = 0, i = 1, . . . , n + 1, 1 i < j n + 1, j = 1, . . . , n + 1, 1 j < i n + 1.
=
1k1 ki-j j i
tk1 tk2 tki-j , 1 j < i n + 1, 1 < i j n + 1,
uij =
k=1
(tj - tk ),
u1j = 1, uij = 0,
Full details of the proof that V = L U can be found in [20, 41]. (Surprisingly, as far as we know, these are the first places this factorization appears in the literature.) The entries of L lying below the diagonal are known as the complete monomial polynomials since ij is obtained by summing, with unit coefficients, all monomials of degree i-j in the j variables t1 , . . . , tj . The entries of U appearing on or above the diagonal are known as the Newton difference polynomials. In particular, if t1 , . . . , tn are distinct, so ti = tj for i = j, all entries of U lying on or above the diagonal are nonzero. In this case, V has all nonzero pivots, and is a regular, hence nonsingular matrix. Q.E.D. Approximation and Interpolation by General Functions There is nothing special about polynomial functions in the preceding approximation scheme. For example, suppose we were interested in finding the best trigonometric approximation y = 1 cos t + 2 sin t to a given set of data. Again, the least squares error takes the same form y - A x in (13.17), where cos t y1 sin t1 1 y2 cos t2 sin t2 1 , A= x= , y = . . . . . . 2 . . . . ym cos tm sin tm 3/15/06 224
c 2006 2
as
Peter J. Olver
Thus, the columns of A are the sampled values of the functions cos t, sin t. The key is that the unspecified parameters -- in this case 1 , 2 -- occur linearly in the approximating function. Thus, the most general case is to approximate the data (13.6) by a linear combination y(t) = 1 h1 (t) + 2 h2 (t) + + n hn (t) of prescribed functions h1 (x), . . . , hn (x). The least squares error is, as always, given by
m
Error =
i=1
yi - y(ti )
2
=
y - Ax ,
If A is square and nonsingular, then we can find an interpolating function of the prescribed form by solving the linear system A x = y. (13.28) A particularly important case is provided by the 2 n + 1 trigonometric functions 1, cos x, sin x, cos 2 x, sin 2 x, ... cos n x, sin n x.
where the sample matrix A, are h1 (t1 ) h2 (t1 ) h1 (t2 ) h2 (t2 ) A= . . . . . . h1 (tm ) h2 (tm )
the vector of unknown coefficients x, and the data vector y ... ... .. . ... hn (t1 ) 1 2 x = . , . . n y1 y2 . . . ym
hn (t2 ) , . . . hn (tm )
y=
.
(13.27)
Interpolation on 2 n + 1 equally spaced data points on the interval [ 0, 2 ] leads to the Discrete Fourier Transform, used in signal processing, data transmission, and compression. If there are more than n data points, then we cannot, in general, interpolate exactly, and must content ourselves with a least squares approximation that minimizes the error at the sample points as best it can. The least squares solution to the interpolation equations (13.28) is found by solving the associated normal equations K x = f , where the (i, j) entry of K = AT A is m times the average sample value of the product of hi (t) and hj (t), namely
m
kij = m hi (t) hj (t) =
=1
hi (t ) hj (t ),
(13.29)
whereas the ith entry of f = AT y is
m
fi = m hi (t) y =
=1
hi (t ) y .
(13.30)
The one issue is whether the columns of the sample matrix A are linearly independent. This is more subtle than the polynomial case covered by Lemma 13.6. Linear independence 3/15/06 225
c 2006 Peter J. Olver
of the sampled function vectors is, in general, more restrictive than merely requiring the functions themselves to be linearly independent. If the parameters do not occur linearly in the functional formula, then we cannot use linear algebra to effect a least squares approximation. For example, one cannot determine the frequency , the amplitude r, and the phase shift of the general trigonometric approximation y = c1 cos t + c2 sin t = r cos( t + ) that minimizes the least squares error at the sample points. Approximating data by such a function constitutes a nonlinear minimization problem. Weighted Least Squares Another extension to the basic least squares method is to introduce weights in the measurement of the error. Suppose some of the data is known to be more reliable or more significant than others. For example, measurements at an earlier time may be more accurate, or more critical to the data fitting problem, than later measurements. In that situation, we should penalize any errors in the earlier measurements and downplay errors in the later data. In general, this requires the introduction of a positive weight c i > 0 associated to each data point (ti , yi ); the larger the weight, the more vital the error. For a straight line approximation y = + t, the weighted least squares error is defined as
m m
Error =
i=1
ci e2 = i
i=1
ci yi - ( + ti )
2
.
Let us rewrite this formula in matrix form. Let C = diag (c1 , . . . , cm ) denote the diagonal weight matrix . Note that C > 0 is positive definite, since all the weights are positive. The least squares error, Error = eT C e = e , is then the norm of the error vector e with respect to the weighted inner product v , w = vT C w. Since e = y - A x, e
2
= Ax - y
2
where
= xT AT C A x - 2 xT AT C y + yT C y = xT K x - 2 xT f + c, f = AT C y, c = yT C y = y
2
= (A x - y)T C (A x - y)
(13.31)
K = AT C A,
.
Note that K is the weighted Gram matrix derived in (12.10), and so is positive definite provided A has linearly independent columns or, equivalently, has rank n. Theorem 13.13. Suppose A is an m n matrix with linearly independent columns. Suppose C > 0 is any positive definite m m matrix. Then, the quadratic function (13.31) giving the weighted least squares error has a unique minimizer, which is the solution to the weighted normal equations AT C A x = AT C y, 3/15/06 so that 226 x = (AT C A)-1 AT C y.
c 2006
(13.32)
Peter J. Olver
12 10 8 6 4 2 1 2 3 4 5 6
Figure 13.8.
Weighted Least Squares Line.
In brief, the weighted least squares solution is obtained by multiplying both sides of the original system A x = y by the matrix AT C. The derivation of this result allows C > 0 to be any positive definite matrix. In applications, the off-diagonal entries of C can be used to weight cross-correlation terms in the data, although this extra freedom is rarely used in practice. Example 13.14. In Example 13.4, we fit the following data ti yi ci 0 2 3 1 3 2 3 7
1 2
6 12
1 4
with an unweighted least squares line. Now we shall assign the weights listed in the last row of the table for the error at each sample point. Thus, errors in the first two data values carry more weight than the latter two. To find the weighted least squares line y = + t that best fits the data, we compute 1 0 3 0 0 0 23 5 1 1 1 1 0 2 0 01 1 T 4 A CA = , = 1 31 1 3 0 0 2 0 0 1 3 6 5 2 1 6 0 0 0 1 4 3 0 0 0 2 37 1 1 1 1 0 2 0 0 3 2 AT C y = = . 1 69 0 1 3 6 0 0 2 0 7 2 0 0 0 1 12 4 Thus, the weighted normal equations (13.32) reduce to
23 4
+ 5 =
37 2 ,
5 +
31 2
=
69 2 ,
so
= 1.7817,
= 1.6511.
Therefore, the least squares fit to the data under the given weights is y = 1.7817 + 1.6511 t, as plotted in Figure 13.8. 3/15/06 227
c 2006 Peter J. Olver
13.3. Splines.
Polynomials are but one of the options for interpolating data points by smooth functions. In preCAD (computer aided design) draftsmanship, a spline was a long, thin, flexible strip of wood that was used to draw a smooth curve through prescribed points. The points were marked by small pegs, and the spline rested on the pegs. The mathematical theory of splines was first developed in the 1940's by the Romanian mathematician Isaac Schoenberg as an attractive alternative to polynomial interpolation and approximation. Splines have since become ubiquitous in numerical analysis, in geometric modeling, in design and manufacturing, in computer graphics and animation, and in many other applications. We suppose that the spline coincides with the graph of a function y = u(x). The pegs are fixed at the prescribed data points (x0 , y0 ), . . . , (xn , yn ), and this requires u(x) to satisfy the interpolation conditions u(xj ) = yj , j = 0, . . . , n. (13.33)
The mesh points x0 < x1 < x2 < < xn are distinct and labeled in increasing order. The spline is modeled as an elastic beam, [38], and so xj x xj+1 , u(x) = aj + bj (x - xj ) + cj (x - xj )2 + dj (x - xj )3 , (13.34) j = 0, . . . , n - 1, is a piecewise cubic function -- meaning that, between successive mesh points, it is a cubic polynomial, but not necessarily the same cubic on each subinterval. The fact that we write the formula (13.34) in terms of x - xj is merely for computational convenience. Our problem is to determine the coefficients aj , bj , cj , dj , j = 0, . . . , n - 1.
Since there are n subintervals, there are a total of 4 n coefficients, and so we require 4 n equations to uniquely prescribe them. First, we need the spline to satisfy the interpolation conditions (13.33). Since it is defined by a different formula on each side of the mesh point, this results in a total of 2 n conditions: u(x+ ) = aj = yj , j u(x- ) = aj + bj hj + cj h2 + dj h3 = yj+1 , j j j+1 where we abbreviate the length of the j th subinterval by hj = xj+1 - xj . The next step is to require that the spline be as smooth as possible. The interpolation conditions (13.35) guarantee that u(x) is continuous. The condition u(x) C 1 be continuously differentiable requires that u (x) be continuous at the interior mesh points x1 , . . . , xn-1 , which imposes the n - 1 additional conditions bj + 2 cj hj + 3 dj h2 = u (x- ) = u (x+ ) = bj+1 , j j+1 j+1 3/15/06 228 j = 0, . . . , n - 2.
c 2006
j = 0, . . . , n - 1,
(13.35)
(13.36)
Peter J. Olver
To make u C2 , we impose n - 1 further conditions 2 cj + 6 dj hj = u (x- ) = u (x+ ) = 2 cj+1 , j+1 j+1 to ensure that u is continuous at the mesh points. We have now imposed a total of 4 n - 2 conditions, namely (13.3537), on the 4 n coefficients. The two missing constraints will come from boundary conditions at the two endpoints, namely x0 and xn . There are three common types: (i ) Natural boundary conditions: u (x0 ) = u (xn ) = 0, whereby c0 = 0, cn-1 + 3 dn-1 hn-1 = 0. (13.38) Physically, this models a simply supported spline that rests freely on the first and last pegs. (ii ) Clamped boundary conditions: u (x0 ) = , u (xn ) = , where , , which could be 0, are fixed by the user. This requires b0 = , bn-1 + 2 cn-1 hn-1 + 3 dn-1 h2 = . n-1 (13.39) j = 0, . . . , n - 2, (13.37)
This corresponds to clamping the spline at prescribed angles at each end. (iii ) Periodic boundary conditions: u (x0 ) = u (xn ), u (x0 ) = u (xn ), so that b0 = bn-1 + 2 cn-1 hn-1 + 3 dn-1 h2 , n-1 c0 = cn-1 + 3 dn-1 hn-1 . (13.40) If we also require that the end interpolation values agree, u(x0 ) = y0 = yn = u(xn ), (13.41) then the resulting spline will be a periodic C2 function, so u(x+p) = u(x) with p = xn -x0 for all x. The periodic case is used to draw smooth closed curves; see below. Theorem 13.15. Suppose we are given mesh points a = x0 < x1 < < xn = b, and corresponding data values y0 , y1 , . . . , yn , along with one of the three kinds of boundary conditions (13.38), (13.39), or (13.40). Then there exists a unique piecewise cubic spline function u(x) C2 [ a, b ] that interpolates the data, u(x0 ) = y0 , . . . , u(xn ) = yn , and satisfies the boundary conditions. Proof : We first discuss the natural case. The clamped case is left as an exercise for the reader, while the slightly harder periodic case will be treated at the end of the section. The first set of equations in (13.35) says that aj = y j , Next, (13.3738) imply that dj = j = 0, . . . , n - 1. cj+1 - cj . 3 hj (13.42)
(13.43)
This equation also holds for j = n - 1, provided that we make the convention that cn = 0.
This is merely for convenience; there is no cn used in the formula for the spline.
3/15/06
229
c 2006
Peter J. Olver
We now substitute (13.4243) into the second set of equations in (13.35), and then solve the resulting equation for bj = (2 cj + cj+1 ) hj yj+1 - yj - . hj 3 (13.44)
Substituting this result and (13.43) back into (13.36), and simplifying, we find hj cj + 2 (hj + hj+1 ) cj+1 + hj+1 cj+2 = 3 yj+1 - yj yj+2 - yj+1 - hj+1 hj = zj+1 , (13.45)
where we introduce zj+1 as a shorthand for the quantity on the right hand side. In the case of natural boundary conditions, we have c0 = 0, cn = 0,
and so (13.45) constitutes a tridiagonal linear system A c = z, for the unknown coefficients c = c1 , c2 , . . . , cn-1
T
(13.46) , with coefficient matrix
(13.47) and right hand side z = z1 , z2 , . . . , zn-1 . Once (13.47) has been solved, we will then use (13.4244) to reconstruct the other spline coefficients aj , bj , dj . The key observation is that the coefficient matrix A is strictly diagonally dominant, cf. Definition 6.25, because all the hj > 0, and so
T
2 (h + h ) h1 0 1 h1 2 (h1 + h2 ) h2 h2 2 (h2 + h3 ) h3 A= .. .. .. . . . hn-3 2 (hn-3 + hn-2 ) hn-2 hn-2 2 (hn-2 + hn-1 )
2 (hj-1 + hj ) > hj-1 + hj . Theorem 6.26 implies that A is nonsingular, and hence the tridiagonal linear system has a unique solution c. This suffices to prove the theorem in the case of natural boundary conditions. Q.E.D. To actually solve the linear system (13.46), we can apply our tridiagonal solution algorithm (4.47). Let us specialize to the most important case, when the mesh points are equally spaced in the interval [ a, b ], so that xj = a + j h, 3/15/06 where h = hj = 230 b-a , n j = 0, . . . , n - 1.
c 2006 Peter J. Olver
2 1.5 1 0.5 1 -0.5 -1 2 3 4
Figure 13.9.
A Cubic Spline.
In this case, the coefficient matrix A = h B is equal to h times the tridiagonal matrix 4 1 B= 1 4 1 1 4 1 1 4 1 1 4 1 .. .. .. . . .
that first appeared in Example 4.26. Its L U factorization takes on an especially simple form, since most of the entries of L and U are essentially the same decimal numbers. This makes the implementation of the Forward and Back Substitution procedures almost trivial. Figure 13.9 shows a particular example -- a natural spline passing through the data points (0, 0), (1, 2), (2, -1), (3, 1), (4, 0). As with the Green's function for the beam, the human eye is unable to discern the discontinuities in its third derivatives, and so the graph appears completely smooth, even though it is, in fact, only C2 . In the periodic case, we set an+k = an , bn+k = bn , cn+k = cn , dn+k = dn , zn+k = zn .
With this convention, the basic equations (13.4245) are the same. In this case, the coefficient matrix for the linear system A c = z, 3/15/06 with c = c0 , c1 , . . . , cn-1 231
T
,
z = z0 , z1 , . . . , zn-1
c 2006
T
,
Peter J. Olver
Figure 13.10.
Three Sample Spline Letters.
is of circulant tridiagonal form: 2 (h h0 hn-1 n-1 + h0 ) h0 2 (h0 + h1 ) h1 h1 2 (h1 + h2 ) h2 . A= .. .. .. . . . hn-3 2 (hn-3 + hn-2 ) hn-2 hn-1 hn-2 2 (hn-2 + hn-1 )
(13.48) Again A is strictly diagonally dominant, and so there is a unique solution c, from which one reconstructs the spline, proving Theorem 13.15 in the periodic case. One immediate application of splines is curve fitting in computer aided design and T graphics. The basic problem is to draw a smooth parametrized curve u(t) = ( u(t), v(t) ) T that passes through a set of prescribed data points xk = ( xk , yk ) in the plane. We have the freedom to choose the parameter value t = tk when the curve passes through the k th point; the simplest and most common choice is to set tk = k. We then construct the functions x = u(t) and y = v(t) as cubic splines interpolating the x and y coordinates of the data points, so u(tk ) = xk , v(tk ) = yk . For smooth closed curves, we require that both splines be periodic; for curves with ends, either natural or clamped boundary conditions are used. Most computer graphics packages include one or more implementations of parametrized spline curves. The same idea also underlies modern font design for laser printing and typography (including the fonts used in this book). The great advantage of spline fonts over their bitmapped counterparts is that they can be readily scaled. Some sample letter shapes parametrized by periodic splines passing through the indicated data points are plotted in Figure 13.10. Better fits can be easily obtained by increasing the number of data points. Various extensions of the basic spline algorithms to space curves and surfaces are an essential component of modern computer graphics, design, and animation, [15, 44].
3/15/06
232
c 2006
Peter J. Olver