Lecture 14: Interpolation

The linear algebra concepts arising from study of linear mappings between vector spaces $𝒇:U\to V$, $𝒇(\alpha 𝒖+\beta 𝒗)=\alpha f\left(𝒖\right)+\beta 𝒇\left(𝒗\right)$, are widely applicable to nonlinear functions also. The study of nonlinear approximation starts with the simplest case of approximation of a function with scalar values and arguments, $f:ℝ\to ℝ$ through additive corrections.

1.Function spaces

An immediate application of the linear algebra framework is to construct vector spaces of real functions $ℱ=(F,+,\cdot )$, with $F=\left\{f|f:ℝ\to ℝ.\right\}$, and the addition and scaling operations induced from $ℝ$,

Comparing with the real vector space $({ℝ}^m,+,\cdot )$ in which the analogous operation is $\alpha 𝒖+\beta 𝒗,𝒖,𝒗\in {ℝ}^m,\alpha ,\beta \in ℝ$, or componentwise

the key difference that arises is the dimension of the set of vectors. Finite-dimensional vectors within ${ℝ}^m$ can be regarded as functions defined on a finite set $𝒖\iff u:\{1,2,\dots ,m\}\to ℝ$, with $u\left(i\right)=u_i$. The elements of $F$ are however functions defined on $ℝ$, a set with cardinality $𝔠=2^{{\aleph }_0}$, with ${\aleph }_0$ the cardinality of the naturals $\mathbb{N}$. This leads to a review of the concept of a basis for this infinite-dimensional case.

1.1.Infinite dimensional basis set

In the finite dimensional case $𝑩\in {ℝ}^{m\times m}$ constituted a basis if any vector $𝒚\in {ℝ}^m$ could be expressed uniquely as a linear combination of the column vectors of

While the above finite sum is well defined, there is no consistent definition of an infinite sum of vectors. As a simple example, in the vector space of real numbers ${ℛ}_1=(ℝ,+,\cdot )$, any finite sum of reals is well defined, for instance

but the limit $S_{n\to \infty }$ cannot be determined. This leads to the necessity of seeking finite-dimensional linear combinations to span a vector space $𝒱=(V,S,+,\cdot )$. First, define linear independence of an infinite (possibly uncountable) set of vectors $ℬ=\{v_{\gamma }|.\gamma \in \Gamma ,v_{\gamma }\in V\}$, where $\Gamma$ is some indexing set.

Definition. The vector set $ℬ=\{v_{\gamma }|.\gamma \in \Gamma ,v_{\gamma }\in V\},$is linearly independent if the only $n\in \mathbb{N}$ scalars, $x_1,\dots ,x_n\in S$, that satisfy

are $x_1=0$, $x_2=0$,…,$x_n=0$.

The important aspect of the above definition is that all finite vector subsets are linearly independent. The same approach is applied in the definition of a spanning set.

Definition. Vectors within the set $ℬ=\{v_{\gamma }|.\gamma \in \Gamma ,v_{\gamma }\in V\},$span $V$, stated as $V=\mathrm{span}\left(ℬ\right)$, if for any $u\in V$ there exist $n\in \mathbb{N}$ scalars, $x_1,\dots ,x_n\in S$, such that

This now allows a generally-applicable definition of basis and dimension.

Definition. The vector set $ℬ=\{v_{\gamma }|.\gamma \in \Gamma ,v_{\gamma }\in V\}$ is a basis for vector space $𝒱=(V,S,+,\cdot )$ if

Definition. The dimension of a vector space $𝒱=(V,S,+,\cdot )$ is the cardinality of a basis set $ℬ$, $\dim \left(𝒱\right)=\left|ℬ\right|$.

The use of finite sums to define linear independence and bases is not overly restrictive since it can be proven that every vector space has a basis. The proof of this theorem is based on Zorn's lemma from set theory, and asserts exsitence of a basis, but no constructive procedure. The difficulty in practical construction of bases for infinite dimensional vector spaces is ascertained through basic examples.

Example. ${ℛ}_{\infty }$. As a generalization of ${ℛ}_m=({ℝ}^m,ℝ,+,\cdot )$, consider the vector space of real sequences ${\left\{x_n\right\}}_{n\in \mathbb{N}}$ represented as a vectors with a countably infinite number of components $𝒙={\left[\begin{array}{llll}{x}_1& x_2& x_3& \dots \end{array}\right]}^T$. Linear combinations are defined by

Let ${𝒆}_i$ denote the vector of all zeros except the $i^{\mathrm{th}}$ position. In ${ℝ}^m$, the identity matrix $𝑰=\left[\begin{array}{lll}{𝒆}_1& \dots & {𝒆}_m\end{array}\right]$ was a basis, but this does not generalize to ${ℝ}^{\infty }$; for example the vector $𝒗={\left[\begin{array}{llll}1& 1& 1& \dots \end{array}\right]}^T$ cannot be obtained by finite linear combination of the ${𝒆}_i$ vectors. In fact, there is no countable set of vectors that spans ${ℝ}^{\infty }$.

Example. $P\left(ℝ\right)$. The vector space of polynomials $P\left(ℝ\right)=\{p|p\left(t\right)=c_n{t}^n+c_{n-1}{t}^{n-1}+\cdots +c_0.,n\in \mathbb{N},c_i\in ℝ,i=0,1,\dots ,N\}$ on the real line has an easily constructed basis, namely the set of the monomials

an infinite set with the cardinality as the naturals $\left|ℬ\right|=\left|\mathbb{N}\right|={\aleph }_0$.

1.2.Alternatives to the concept of a basis

The difficulty in ascribing significance to an infinite sum of vectors ${\sum }_{i=1}^{\infty }{𝒗}_i$ can be resolved by endowing the vector space with additional structure, in particular a way to define convergence of the partial sums

to a limit ${lim}_{n\to \infty }{𝒔}_n=𝒗$.

Fourier series.

and norm $||f||={(f,f)}^{1/2}$. An orthonormal spanning set for $C[-\pi ,\pi ]$ is given by

Vector spaces with an inner product are known as Hilbert spaces.

Taylor series.

for which a spanning set is given by the monomials $\{1,t,t^2,\dots \}$, and the infinite exapnsion

is convergent, with coefficients given by the Taylor series

Note that orthogonality of the spanning set cannot be established, absent an inner product.

1.3.Common function spaces

Several function spaces find widespread application in scientific computation. An overview is provided in Table 1.

2.Interpolation

The interpolation problem seeks the representation of a function $f$ known only through a sample data set $𝒟=\left\{(x_i,y_i=f\left(x_i\right))|i=0,\dots ,m.\right\}\subset ℝ\times ℝ$, by an approximant $p\left(t\right)$, obtained through combination of elements from some family of computable functions, $ℬ=\{b_0,\dots ,b_n|b_k:ℝ\to ℝ.\}$. The approximant $p\left(t\right)$ is an interpolant of $𝒟$ if

or $p\left(t\right)$ passes through the known poles $(x_i,y_i)$ of the function $f$. The objective is to use $p\left(t\right)$ thus determined to approximate the function $f$ at other points. Assuming $x_0<x_1<\cdots <x_m$, evaluation of $p\left(t\right)$ at $t\in (x_0,x_m)$ is an interpolation, while evaluation at $t<x_0$ or $t>x_m$, is an extrapolation. The basic problems arising in interpolation are:

• choice of the family from which to build the approximant $p\left(t\right)$;

• choice of the combination technique;

• estimation of the error of the approximation given some knowledge of $f$.

git config --global user.name "SorinMitran"

git config --global user.email "mitran@unc.edu"

∴ 

home=unsafe_string(ccall(:getenv, Cstring, (Cstring,), "HOME"));

∴ 

pkgdir = home * "/courses/MATH661/packages";

∴ 

∴ 

using PkgTemplates

∴ 

tmpl=Template(;
  user="SorinMitran", dir=pkgdir,
  plugins = [ License(; name="MPL"),
              Git(; manifest=true, ssh=true),
              GitHubActions(; x86=true),
              Codecov(),
              Documenter{GitHubActions}(),
              Develop(),
             ],
);

∴ 

tmpl("Interpolation661");

∴ 

2.1.Additive corrections

As to be expected, a widely used combination technique is linear combination,

The idea is to capture the nonlinearity of $f\left(t\right)$ through the functions $b_0\left(t\right),\dots ,b_n\left(t\right)$, while maintaining the framework of linear combinations. Sampling of $b_j\left(t\right)$ at the poles $x_i$ of a data set $𝒟$, constructs the vectors

which gathered together into a matrix leads to the formulation of the interpolation problem as

Before choosing some specific function set $ℬ$, some general observations are useful.

1. The function values $y_i=f\left(x_i\right),i=0,\dots ,m$, are directly incorporated into the interpolation problem (3). Any estimate of the error at other points requires additional information on $f$. Such information can be furnished by bounds on the function values, or knowledge of its derivatives for example.

2. A solution to (3) exists if $𝒚\in C\left(𝑩\right)$. Economical interpolations would use $n<m$ functions in the set $ℬ$, hopefully $n\ll m$.

2.2.Polynomial interpolation

Monomial form of interpolating polynomial.

With ${ℳ}_{n+1}\left(t\right)$ denoting the first $n+1$ monomials

a polynomial of degree $n$ is the linear combination

Let $𝑴\in {ℝ}^{(m+1)\times (n+1)}$ denote the matrix obtained from evaluation of the first $n+1$ monomials at the sample points $𝒙={\left[\begin{array}{llll}{x}_0& x_1& \dots & x_m\end{array}\right]}^T$,

The above notation conveys that a finite-dimensional matrix $𝑴\in {ℝ}^{(m+1)\times (n+1)}$ is obtained from evaluation of the row vector of the monomial basis function $ℳ\left(x\right):ℝ\to {ℝ}^{n+1}$, at the column vector of sample points $𝒙\in {ℝ}^{m+1}$. The interpolation condition $p\left(𝒙\right)=𝒚$ leads to the linear system

For a solution to exist for arbitrary $𝒚$, $𝑴$ must be of full rank, hence $m=n$, in which case $𝑴$ becomes the Vandermonde matrix

known to be ill-conditioned. Since $𝑴$ is square and of full rank, (4) has a unique solution.

Finding the polynomial coefficients by solving the above linear system requires $𝒪(n^3/3)$ operations. Evaluation of the monomial form is economically accomplished in $𝒪\left(n\right)$ operations through Horner's scheme

$\circ$

Figure 1. Monomial basis over interval $[-\pi ,\pi ]$

∴	n=6; x=range(-π,π,length=n+1); y=sin.(x);

∴	t=range(-π,π,length=10*n);

∴	M=ones(length(t),length(x));

∴	for j=2:n+1 M[:,j] = t.^(j-1) end

∴	figure(1); clf();

∴	for j=1:n+1 plot(t,M[:,j]) end

∴	xlabel("t"); ylabel("m(j,t)"); grid("on"); title("Monomial basis");

∴	imdir=home*"/courses/MATH661/images/";

∴	savefig(imdir*"L18MonomialBasis.eps")

∴

∴ 

function Horner(t,a)
   n=length(a)-1; p=a[n+1]
   for k=n:-1:1
     p=a[k]+p*t
   end
   return p
end;

∴ 

p2(t)=3*t^2-2*t+1; ap2=[1 -2 3];

∴ 

t=-3:3; [p2.(t) Horner.(t,Ref(ap2))]

∴ 

Lagrange form of interpolating polynomial.

Let ${ℒ}_{n+1}\left(t\right)=\left[\begin{array}{llll}{\ell }_0\left(t\right)& {\ell }_1\left(t\right)& \dots & {\ell }_n\left(t\right)\end{array}\right]$ denote another set of basis functions that evaluates to the identity matrix at the sample points $𝒙$, such that ${ℒ}_{n+1}\left(𝒙\right)=𝑰$,

For arbitrary $t$, the relationship

describes a linear mapping between the monomials ${ℳ}_{n+1}\left(t\right)$ and the ${ℒ}_{n+1}\left(t\right)$ functions, a mapping which is invertible since $𝑴=𝑳𝑼$ is of full rank

Note that organization of bases as row vectors of functions leads to linear mappings expressed through right factors.

In[37]:= 

V[n_]:=Table[Subscript[x, i]^j,{i, 0, n},{j, 0, n}];
n=3; n1=n+1; M=V[n]

In[13]:= 

LU=Simplify[LUDecomposition[M]][[1]]

In[14]:= 

L=Table[ If[i<j,LU[[j,i]], If[i==j, 1, 0]],{j,1,n1},{i,1,n1}]

In[15]:= 

U=Table[ If[i>=j,LU[[j,i]], 0],{j,1,n1},{i,1,n1}]

In[18]:= 

Simplify[L.U]

In[38]:= 

Factor[Det[M]]

In[39]:= 

In[19]:= 

Linv=Simplify[Inverse[L]]

In[20]:= 

Uinv=Simplify[Inverse[U]]

In[21]:= 

In[23]:= 

Simplify[{1,t,t^2,t^3}.Uinv.Linv]

In[24]:= 

By construction, through the condition ${ℒ}_{n+1}\left(𝒙\right)=𝑰$, a Lagrange basis function evaluated at a sample point is

A polynomial of degree $n$ is expressed as a linear combinations of the Lagrange basis functions by

The interpolant of data $\{(x_i,y_i=f\left(x_i\right)),i=0,1,\dots ,n\}$ is determined through the conditions

i.e., the linear combination coefficients are simply the sampled function values $c_i=y_i=f\left(x_i\right)$.

Determining the linear combination coefficients may be without cost, but evaluation of the Lagrange form (7) of the interpolating polynomial requires $𝒪\left(n^2\right)$ operations, significantly more costly than the $𝒪\left(n\right)$ operations required by Horner's scheme (5)

∴ 

function Lagrange(t,x,y)
  n=length(x)-1; p=0
  for i=1:n+1
    w=1
    for j=1:n+1
      if (i!=j) w=w*(t-x[j])/(x[i]-x[j]); end
    end
    p = p + w*y[i]
  end
  return p
end;

∴ 

p2(t)=3*t^2-2*t+1;

∴ 

x=[-2 0 2]; y=p2.(x);

∴ 

t=-3:3; [p2.(t) Lagrange.(t,Ref(x),Ref(y))]

∴ 

$•$

Figure 2. Lagrange basis for $n=6$ for $\sin \left(x\right)$ over interval $[-\pi ,\pi ]$

∴	n=6; x=range(-π,π,length=n+1); y=sin.(x);

∴	t=range(-π,π,length=10*n);

∴	M=ones(length(t),length(x));

∴	function lagrange(j,t,x,y) n=length(x)-1; l=1 for i=1:n+1 if (i!=j) l = l*(t-x[i])/(x[j]-x[i]); end end return l end;

∴

∴	for j=1:n+1 M[:,j] = lagrange.(j,t,Ref(x),Ref(y)) end

∴	figure(1); clf();

∴	for j=1:n+1 plot(t,M[:,j]) end

∴	xlabel("t"); ylabel("l(i,t)"); grid("on"); title("Lagrange basis");

∴	imdir=home*"/courses/MATH661/images/";

∴	savefig(imdir*"L18LagrangeBasis.eps")

∴

A reformulation of the Lagrange basis can however reduce the operation count. Let $\ell \left(t\right)={\prod }_{k=0}^n(t-x_k)$, and rewrite ${\ell }_i\left(t\right)$ as

with the weights

depending only on the function sample arguments $x_i$, but not on the function values $y_i$. The interpolating polynomial is now

Interpolation of the function $g\left(t\right)=1$ would give

and taking the ratio yields

known as the barycentric Lagrange formula (by analogy to computation of a center of mass). Evaluation of the weights $w_i$ costs $𝒪\left(n^2\right)$ operations, but can be done once for any set of $x_i$. The evaluation of $p\left(t\right)$ now becomes an $𝒪\left(2n\right)$ process, comparable in cost to Horner's scheme.

∴ 

function BaryLagrange(t,x,y)
  n=length(x)-1; w=ones(size(x));
  for i=1:n+1
    w[i]=1
    for j=1:n+1
      if (i!=j) w[i]=w[i]/(x[i]-x[j]); end
    end
  end
  q=r=0
  for i=1:n+1
    d=t-x[i]
    if d≈0 return y[i]; end
    s=w[i]/d; q=q+y[i]*s; r=r+s
  end
  return q/r
end;

∴ 

p2(t)=3*t^2-2*t+1;

∴ 

x=[-2 0 2]; y=p2.(x);

∴ 

t=-3:3; [p2.(t) BaryLagrange.(t,Ref(x),Ref(y))]

∴ 

Newton form of interpolating polynomial.