Approximation Theory

In Interpolation, we were given a set of points $\{(x_i,y_i){\}}_{i=1}^n$, $y_i=f(x_i)$, and asked to construct a polynomial $P(x)$ exactly matching these points - more formally,

$$P(x_i)=f(x_i)\forall 1\le i\le n$$

But this exact fit nature of interpolation isn’t always the best! In the real world, error is inevitable, and our interpolation would fit these errors precisely! In practical applications, this can give us really poor predictions.

Below, we examine various ways we can approximate our functions from given datapoints.

Approximating $f$ Given Noisy Data

Problem Defined

Suppose we have noisy data. How can we create a function to best approximate this data?

In general, we can express our datapoints as

$$y_i=f(x_i)+{ϵ}_i$$

Where ${ϵ}_i$ is some sort of error offseting our true value.

Pick $m$ “known” set of functions $\{{\varphi }_i(x){\}}_{i=1}^m$. For constants $c_1,\dots c_m$, we want to find approximating function

$$f(x)=c_1{\varphi }_1(x)+c_2{\varphi }_2(x)+\dots c_m{\varphi }_m(x)=\sum _{i=1}^m{c}_j{\varphi }_j(x)$$

Matching our data for all $x\in I$.

Note that to determine $m$ unknowns, we need at least $m$ datapoints. However, in most situations, we’ll have an overdetermined situation, where we have $>m$ data points, meaning we’ll have more conditions than unknowns!

In such situations, it’s not possible to find a “solution” that perfectly fits our data. Thus, we need some notion of best fit! Given our function $f$, we define the residual error at any datapoint $\{x_i,y_i\}$ to be

$$r_i=f(x_i)-y_i=\sum _{j=1}^m{c}_j{\varphi }_j(x_i)-y_{i}1\le i\le n$$

This gives us a residual vector $r\in {\mathbb{R}}^n$

$$r=\left[\begin{array}{c}{r}_1\\ r_2\\ ⋮\\ r_n\end{array}\right]$$

We want to find constants $\{c_j{\}}_{j=1}^m$ such that the magnitude of this residual vector is minimized (meaning we’ve minimized our residual error), given some norm $||\cdot ||$.

This is known as the Linear Least Squares Problem! The solution is described below.

Linear Least Squares Problem

Let $f$ be our approximating function.

For all datapoints $\{x_i,y_i\},1\le i\le n$, we have approximations given by

$$\begin{array}{rl}{\tilde{y}}_1=& f(x_1)=c_1{\varphi }_1(x_1)+c_2{\varphi }_2(x_1)+\cdots +c_m{\varphi }_m(x_1)\\ {\tilde{y}}_2=& f(x_2)=c_1{\varphi }_1(x_2)+c_2{\varphi }_2(x_2)+\cdots +c_m{\varphi }_m(x_2)\\ & ⋮\end{array}$$

We can represent these approximations in the matrix form $Ac,$ where $A_{ij}={\varphi }_j(x_i)$, and $c$ is our constant vector.

$$A=\left[\begin{array}{cccc}{\varphi }_1(x_1)& {\varphi }_2(x_1)& \dots & {\varphi }_m(x_1)\\ {\varphi }_1(x_2)& {\varphi }_2(x_2)& \dots & {\varphi }_m(x_2)\\ ⋮& ⋮& & ⋮\\ {\varphi }_1(x_n)& {\varphi }_2(x_n)& \dots & {\varphi }_m(x_n)\end{array}\right]c=\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_m\end{array}\right]$$

Now, given our vector of true datapoints $y=(y_1,y_2,\dots y_n{)}^T$, we can define our residual error in the matrix form

$$r=Ac-y$$

As mentioned before, we want to minimize the magnitude of this vector! We’ll use the $l_2$ norm to minimize our error. How can we minimize $||r|{|}_2=\sqrt{{\sum }_{i=1}^n{r}_i^2}$?

Let’s define a function $ϵ(c)$, where $c$ are our constants as

$$\begin{array}{rl}ϵ(c)& =(||r|{|}_2{)}^2=\sum _{i=1}^n(r_i(c){)}^2=(||Ac-y|{|}_2{)}^2\\ & =(Ac-y{)}^T(Ac-y)=c^T{A}^{T}Ac-2y^{T}Ac+y^{T}y\end{array}$$

Note that when minimizing these distances, it’s often the case we square the term to get rid of the square root.

We wish to find minimizing $c$, and to do this, we can find when our derivative is 0! Differentiating with respect to any $c_j$, we find

$$\begin{array}{rl}\frac{\partial ϵ(c)}{\partial c_j}& =\sum _{i=1}^n\frac{\partial }{\partial c_j}(r_i(c{)}^2)=\sum _{i=1}^{n}2r_i(c)\frac{\partial r_i(c)}{\partial c_j}\\ & =\sum _{i=1}^{n}2r_i{\varphi }_j(x_i)=2(A_{1j}{r}_1+A_{2j}{r}_2+\dots A_{nj}{r}_n)\end{array}$$

Continuing this for all $c_j$, we find that we want to find $c$ such that

$$A^{T}r=0$$

Otherwise known as the normal equations.

We plug this back into our original equation to obtain

$$\begin{array}{r}r=Ac-y\\ A^{T}r=0=A^T(Ac-y)\\ A^{T}Ac=A^{T}r+A^{T}y\\ A^{T}Ac=(0)A^{T}y\end{array}$$

This is a linear system that we can solve for $c$!

Thus, our algorithm for solving the linear least squares problem is as follows:

1. First, construct $A$ from our $x_i$‘s and ${\varphi }_j$‘s given.
2. Then, construct $M=A^{T}A$ and $b=A^{T}y$.
3. Finally, solve $Mc=b$.

Typically, we want to solve the normal equations ourselves, as solving them on a computer can lead to a highly unstable algorithm!

Theorem: Uniqueness of Solutions

If the columns of $A$ are linearly independent, then

1. The normal equations have a unique solution for $c$.
2. If $c$ is this unique solution, then it minimizes the residual error - that is, $||Ac-y||\le ||A{c}^{\ast }-y||$ for all $c\ne c^{\ast }$.

Example: Deriving $f$

Let $m=2$, where

$${\varphi }_1(x)=1{\varphi }_2(x)=x$$

Say we have $n$ points $\{(x_i,y_i){\}}_{i=1}^n$. We want function

$$f(x)=c_1+c_{2}x$$

We are essentially finding a line that minimizes the square of the vertical distances to each point.

To find our constants, we first find $A$ as

$$A=\left[\begin{array}{cc}1& x_1\\ 1& x_2\\ ⋮& \\ 1& x_n\end{array}\right]$$

We use this to find $A^{T}A$, $A^{T}y$ as

$$A^{T}A=\left[\begin{array}{cc}n& {\sum }_{i=1}^n{x}_i\\ {\sum }_{i=1}^n{x}_i& {\sum }_{i=1}^n{x}_i^2\end{array}\right]A^{T}y=\left[\begin{array}{c}{\sum }_{i=1}^n{y}_i\\ {\sum }_{i=1}^n{x}_i{y}_i\end{array}\right]$$

We can use this to solve for $c$!

Linear Least Squares Pros

Note that Linear Least Squares only works well if we know the following about the noise:

• The average noise has no bias; as in, the average is 0.
• The noise is uncorrelated.
• The noise’s variance is approximately the same.

Approximating $f$ Given True $f$

Problem Defined

Suppose we have a known function $f$, but it is really complicated. Thus, in practical applications, it may be worth approximating it with a cleaner solution.

How can we create a polynomial of degree $\le N$ that best approximates our function?

In other words, given functions $\{{\varphi }_i(x){\}}_{i=0}^N$, we want to approximate $f$ as

$$f(x)=c_1{\varphi }_1(x)+c_2{\varphi }_2(x)+\dots c_N{\varphi }_N(x)$$

How can we do this?

This is known as the weighted least squares problem.

Weighted Least Squares Problem

Linear Independence

Consider functions $\{{\varphi }_i(x){\}}_{i=0}^N$ that we want to approximate $f$ with. These functions are linearly independent if

$$c_0{\varphi }_0(x)+c_1{\varphi }_1(x)+\cdots +c_N{\varphi }_N(x)=0$$

is only possible for $c_0=c_1=\cdots =c_N=0$.

Theorem: Independence of Distinct Degrees

If for $0\le i\le N$, ${\varphi }_i(x)$ is a polynomial of degree $=i$, then $\{{\varphi }_i(x){\}}_{i=0}^N$ is linearly independent on any $[a,b]$.

We furthermore define ${\prod }_{i=0}^N[a,b]$ as the set of all polynomials $P(x)$ on $[a,b]$ with degree $\le N$. We can use this to define what polynomials our functions can span.

Theorem: Linear Independence and Span

Let $\{{\varphi }_i{\}}_{i=0}^N$ be a set of linearly independent polynomials in ${\prod }_N[a,b]$. Then, any polynomial in ${\prod }_N[a,b]$ can be uniquely written as a linear combination of $\{{\varphi }_i{\}}_{i=0}^N$.

Weight Functions

We define $w(x)$ to be a weight function on $[a,b]$ if it is continuous and non-negative in $(a,b)$, and

$${\int }_a^{b}w(x)dx>0$$

In other words, $w(x)$ has a positive max.

Given any weight function $w$, the weighted $L^2$ norm of some function $f$ on $[a,b]$ is defined as

$$||f|{|}_{2,w}=\sqrt{{\int }_a^b{f}^2(x)\ast w(x)}$$

Note that if $w(x)=1$, then our weighted $L^2$ norm is the same as a normal $L^2$ norm!

Problem Solution

We can our weighted $L^2$ norm to get a measure of error / distance!

Let $w(x)$ be a weight function, $f$ be the function we want to approximate. We want to construct a polynomial $P_N^{\ast }$ minimizing

$$(||f-P_N|{|}_{2,w}{)}^2={\int }_a^b(f(x)-P_N(x){)}^{2}w(x)dx$$

Or in other words, the weighted $L^2$ distance to $f$.

Let us have a set of polynomials $\{{\varphi }_i(x){\}}_{i=0}^N$ with $\mathrm{deg}({\varphi }_i)=i$. Then, given our earlier theorems, we guarantee that any $P_N(x)\in {\prod }_N[a,b]$ can be uniquely written as

$$P_N(x)=c_0{\varphi }_0(x)+\cdots +c_N{\varphi }_N(x)=\sum _{j=0}^N{c}_j{\varphi }_j(x)$$

Meaning we can use this set of polynomials to construct our $P_N^{\ast }$ mentioned earlier.

Thus, we want to find constants $c^{\ast }=[c_0^{\ast },\dots ,c_N^{\ast }]$ minimizing

$$ϵ(c)={\int }_a^b{\left(f(x)-\sum _{i=0}^N{c}_i{\varphi }_i(x)\right)}^{2}w(x)dx$$

To find this minimum, we’ll find where our derivative is equal to 0. For any $c_i$, we differentiate with respect to $c_i$ to get

$$\begin{array}{rl}0& =\frac{\partial ϵ(c)}{\partial c_i}{|}_{c=c^{\ast }}=-2{\int }_a^b\left(\frac{\partial ϵ(c)}{\partial c_i}f(x)-\frac{\partial ϵ(c)}{\partial c_i}\sum _{j=0}^N{c}_j^{\ast }{\varphi }_j(x)\right){\varphi }_i(x)w(x)dx\\ 0& =\sum _{j=0}^N\left[c_j^{\ast }{\int }_a^b{\varphi }_i(x){\varphi }_j(x)w(x)dx\right]\end{array}$$

Note that we applied the chain rule here!

Continuing this for all $c_i$, we get a linear system of equations!

$$M{c}^{\ast }=d$$

where

$$\begin{array}{r}{M}_{ij}={\int }_a^b{\varphi }_i(x){\varphi }_j(x)w(x)dx\\ d_i={\int }_a^{b}f(x){\varphi }_i(x)w(x)\end{array}$$

Naive Solutions

Making the simplest choice of $w(x)=1$ and ${\varphi }_i=x^i$, we can simplify this system of equations to

$$\begin{array}{r}{M}_{ij}={\int }_a^b{x}^{i+j}dx\\ d_i={\int }_a^{b}f(x)x^i\end{array}$$

If $[a,b]=[0,1]$, we get matrix

$$M=\left[\begin{array}{cccc}1/1& 1/2& \dots & 1/N\\ 1/2& 1/3& & ⋮\\ ⋮& & & \\ 1/(N+1)& \dots & & 1/(2N+1)\end{array}\right]$$

known as a Hilbert Matrix. While this matrix is simple, it’s really terrible for error, as it has a condition number that grows extremely fast!

Is there another matrix we could use that’s a better choice?

Note that for any linear system, the easiest matrix we could use to solve it is a diagonal matrix.

If we wanted this for our system, we’d need to define a set of functions $\{{\varphi }_j{\}}_{j=0}^N$ that are orthogonal on $[a,b]$ with respect to $w(x)$. In other words,

$${\int }_a^b{\varphi }_i(x){\varphi }_j(x)w(x)dx=\{\begin{array}{ll}0& i\ne j\\ {\alpha }_i>0& i=j\end{array}$$

We can furthermore say that these functions are orthonormal if the same conditions hold, but ${\alpha }_i=1$. Note that we can convert any orthogonal set into an orthonormal set by taking

$${\hat{\varphi }}_j(x)=\frac{{\varphi }_j(x)}{\sqrt{{\alpha }_i}}=\frac{{\varphi }_j(x)}{||{\varphi }_j|{|}_{2,w}}$$

We can find such a set using the Gram-Schmidt process of orthonormalization! See the below example.

Gram-Schmidt Orthonormalization (Example)

Suppose we have linearly independent set of functions $\{{\varphi }_0(x),{\varphi }_1(x),{\varphi }_2(x),\dots {\varphi }_N(x)\}$.

Find the corresponding orthonormal set of functions, given a weight function $w(x)$.

1. Start with ${\hat{\varphi }}_0(x)$.
2. Let ${\hat{\varphi }}_1(x)={\varphi }_1(x)-{\beta }_1^0{\hat{\varphi }}_0(x)$. Find ${\beta }_1^0$ such that $${\int }_a^b{\hat{\varphi }}_0(x){\hat{\varphi }}_1(x)w(x)=0$$
3. Let ${\hat{\varphi }}_2(x)={\varphi }_2(x)-{\beta }_2^0{\hat{\varphi }}_0(x)-{\beta }_2^1{\hat{\varphi }}_1(x)$. Find ${\beta }_2^0$ and ${\beta }_1^0$ such that $$\begin{array}{r}{\int }_a^b{\hat{\varphi }}_2(x){\hat{\varphi }}_0(x)w(x)dx=0\\ {\int }_a^b{\hat{\varphi }}_2(x){\hat{\varphi }}_1(x)w(x)dx=0\end{array}$$
4. …

Continue this for all ${\hat{\varphi }}_i(x)$ to obtain an orthogonal set. Finally, orthonormalize each function to finish!

If we choose a set of functions $\{{\varphi }_j{\}}_{j=0}^N$ to be orthonormal, then $M$ becomes the identity matrix, giving us

$$c_i^{\ast }=d_i={\int }_a^b{\varphi }_i(x)f(x)w(x)dx$$

Which makes it super easy to find an approximation for $f$!

Some important examples of these orthonormal polynomials are given as follows:

• The Legendre Polynomials: Given $w(x)=1$, and applying Grahm Schmidt on $[-1,1]$, we have

  $$\begin{gathered} {\varphi }_0(x)=1 \\ {\varphi }_1(x)=x \\ {\varphi }_{n+1}(x)=\frac{(2n+1)x{\varphi }_n(x)-n{\varphi }_{n-1}(x)}{n+1}\forall n\ge 1 \end{gathered}$$

  with $L^2$ norm given as

  $$||{\varphi }_n|{|}_{2,w}=\sqrt{\frac{2}{2n+1}}$$

• The Chebyshev Polynomials: Given $w(x)=\frac{1}{\sqrt{1-x^2}}$ on $[-1,1]$, we get

  $$\begin{gathered} {\varphi }_0(x)=1 \\ {\varphi }_1(x)=x \\ {\varphi }_{n+1}=2x{\varphi }_n(x)-{\varphi }_{n-1}(x)\forall n\ge 1 \end{gathered}$$

  With $L^2$ norm given as

  $$||{\varphi }_n|{|}_{2,w}=\{\begin{array}{ll}\sqrt{\pi }& n=0\\ \sqrt{\frac{\pi }{2}}& n\ne 0\end{array}$$