Interpolation and Least Squares Curve Fitting

Interpolation

Given n+1 points, {(Subscript[x, 0]), (Subscript[y, 0])}, {(Subscript[x, 1]), (Subscript[y, 1])}, {(Subscript[x, 2]), (Subscript[y, 2])},...{(Subscript[x, n]), (Subscript[y, n])}, find the polynomial of degree n that passes through these points.  The polynomial is of the form p(x)=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)+...+(Subscript[a, n])(x^n).  From this polynomial form and from the given points, you get the system of equations:

(Subscript[y, 0])=(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 0])+(Subscript[a, 2])((Subscript[x, 0])^2)+...+(Subscript[a, n])((Subscript[x, 0])^n)
(Subscript[y, 1])=(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 1])+(Subscript[a, 2])((Subscript[x, 1])^2)+...+(Subscript[a, n])((Subscript[x, 1])^n)
(Subscript[y, 2])=(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, 2])+(Subscript[a, 2])((Subscript[x, 2])^2)+...+(Subscript[a, n])((Subscript[x, 2])^n)
...
(Subscript[y, n])=(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, n])+(Subscript[a, 2])((Subscript[x, n])^2)+...+(Subscript[a, n])((Subscript[x, n])^n)

This is a system of n+1 equations with n+1 unknowns (Subscript[a, 0]),(Subscript[a, 1]),...,(Subscript[a, n]), and it can be solved as long as the determinant of the coefficient matrix is not zero.  However, if the points (Subscript[x, 0]),(Subscript[x, 1]),...,(Subscript[x, n]) are distinct this determinant is not zero.  This is illustrated in the following special case in which the determinant is computed and simplified by Mathematica.

Simplify[Det[{{1,x0,x0^2,x0^3},{1,x1,x1^2,x1^3},{1,x2,x2^2,x2^3},
{1,x3,x3^2,x3^3}}]]

All determinants of this form are nonzero because they reduce to a product of terms each of the form xi-xj, and as long as xi is not equal to xj, the system of equations has a unique solution.

Mathematica has a special function for producing these interpolating polynomials.  It is called InterpolatingPolynomial and has the following definition:

For example, to find the interpolating 3rd degree polynomial through the points {1,3},{2,4},{3,5}, and {5,2} use

pts={{1,3},{2,4},{3,5},{5,2}};

ip=InterpolatingPolynomial[pts,x]

Expand[%]

Next, the following three graphs show the points, the interpolating polynomial, and the points and polynomial on the same graph.

g1=ListPlot[pts,PlotStyle->PointSize[0.05]];

g2=Plot[ip,{x,0,6}];

Show[g1,g2];

The Approximation Theorem

If W is a subspace of vector space V then every u in V has a unique expression as u=w1+w2 where w1 is in W and w2 is orthogonal to W.  The projection of u onto W is given by the sum of (u.(Subscript[v, i]))(Subscript[v, i]) for each of the basis vectors (Subscript[v, i]) in an orthonormal basis for W.
If W is a subspace of V and u is in V, then the projection of u onto W is the closest approximation in W to u.  This can be stated as ||u-Proj(u onto W)||<=||u-w|| for any other w in W.  The proof of this statement follows:
u-w=u-Proj(u onto W)+Proj(u onto W)-w, and since the first term in the sum is orthogonal to W while the second term in the sum is in W,
(((||u - w||))^2)=(((||u - Proj(u onto W)||))^2)+(((||Proj(u onto W) - w||))^2).  But then
(((||u - w||))^2) >=\(TraditionalForm\`\(\(\ \)\(\(||\)\(u - 
          Proj(u\ onto\ W)\)\( || \^2\)\)\)\) and the result follows by taking square roots.

Least Squares Line of Best Fit

What straight line best fits the points from the example shown above?

pts={{1,3},{2,4},{3,5},{5,2}};

ListPlot[pts,PlotRange->{{0,6},{0,10}},PlotStyle->PointSize[0.05]];

As is evident from the graph no single straight line passes through all of these points.  When we ask what line best fits them, we ask which line y=a+bx best fits in the sense of minimizing the square root of the sum of squares of the distances between the given points and points on the line.  That is, find the a and b for which

f[x_]:=a + b x

Apply[Plus,Map[#^2&,Map[f[#]&,Transpose[pts][[1]]]-
Transpose[pts][[2]]]]

is a minimum.  This is accomplished by the following statement.

FindMinimum[%,{a,0},{b,0}]

However, there is a linear algebra solution to this problem.  The least squares line of best fit to a set of points is found by the following Mathematica statement.

fl=Fit[pts,{1,x},x]

g3=Plot[fl,{x,0,6}];

Show[g1,g3];

The solution is identical to that found by minimizing the sum of squares.  But what is happening here?

Suppose that the given points are {(Subscript[x, 1]),(Subscript[y, 1])}, {(Subscript[x, 2]),(Subscript[y, 2])}, {(Subscript[x, 3]),(Subscript[y, 3])}, and {(Subscript[x, 4]),(Subscript[y, 4])}.  Then the solution of our problem is to find the solution for the system of equations Y=AX or Y-AX=0.  What are Y, A, and X?

A={{1,x1},{1,x2},{1,x3},{1,x4}}

Y={y1,y2,y3,y4};

X={a,b};

Y-A.X

As a and b vary, the linear combinations of the columns of A form a plane.  If the equations have a solution, vector Y will lie in this plane.  If they don't have a solution, vector Y will not lie in the plane produced by linear combinations of the columns of A.  However, the best alternative to an exact solution is to find the vector in this plane that is closest to Y.  This is the projection of Y onto the vector space formed by the columns of A.  If (Subscript[X, 0]) is the vector that is closest to Y, then Y-A(Subscript[X, 0]) will be orthogonal to any vector in the column space of A.  In other words (Y-A(Subscript[X, 0])).A(Subscript[X, 0])=0.  But this can be written 
Transpose[ A.(Subscript[X, 0])](Y-A.(Subscript[X, 0]))=0 or Transpose[(Subscript[X, 0])]Transpose[A](Y-A.(Subscript[X, 0]))=0 or
Transpose[(Subscript[X, 0])](Transpose[A].Y-Transpose[A].A.(Subscript[X, 0]))=0, and this must hold for all vectors (Subscript[X, 0]).  This means that Transpose[A].Y-Transpose[A].A.(Subscript[X, 0]) is orthogonal to every vector X in s-space, and must therefore be the zero vector.  So Transpose[A].Y-Transpose[A].A.(Subscript[X, 0])=0, or solving for (Subscript[X, 0]) gives 
Transpose[A].Y=Transpose[A].A.(Subscript[X, 0]), and finally 
(Subscript[X, 0])=Inverse[Transpose[A].A].Transpose[A].Y

Transpose[A].A

In order that Transpose[A].A have an inverse, the determinant must be nonzero.  If the x's are not identical, this is the case.  Applying this matrix solution method to the set of points given at the beginning results in the same solution as shown there.

A={{1,1},{1,2},{1,3},{1,5}};

Inverse[Transpose[A] . A]

Y={3,4,5,2};

Inverse[Transpose[A].A].Transpose[A].Y

Next, suppose that we would like to fit a polynomial of degree 2 to the given points.  The only change is in matrix A--everything else proceeds in the same way as before and the solution is given by the same formula.

A={{1,1,1},{1,2,4},{1,3,9},{1,5,25}};

Y={3,4,5,2};

Transpose[A].A

Inverse[Transpose[A].A].Transpose[A].Y

fq=Fit[pts,{1,x,x^2},x]

g4=Plot[fq,{x,0,6}];

Show[g1,g4];

Functions Used in This Notebook

norm[vec_]:=Sqrt[vec.vec]

proj[vec1_,vec2_]:=((vec1.vec2)/(vec2.vec2))vec2

CrossProduct[{a_,b_,c_},{e_,f_,g_}]:={b g - c f,c e - a g,a f - b e}

$Post=If[MatrixQ[#],MatrixForm[#] ,#]&;

gramschmidt[B_]:=
     Block[{i,k},
          V=Table[1,{i,1,Length[B]}];
          V[[1]]=B[[1]];
          Do[V[[i]]=B[[i]]-Sum[proj[B[[i]],V[[k]]],{k,1,i-1}],
               {i,2,Length[B]}];
          Do[V[[i]]=V[[i]]/Sqrt[V[[i]].V[[i]]],
               {i,1,Length[V]}];
          V ]