Math 441

Interpolation: Lagrange and Newton’s Bases

Interpolation Problem

Given \(n + 1\) data points \((x_i,y_i)\), \(i = 0, 1, \dots, n\), find a function \(f\) that belongs to an \(n+1\)-dimensional vector space \(V\) with basis \(\left\{\phi_0, \phi_1, \dots, \phi_n\right\}\) such that \[f(x_i) = y_i \text{ for } i = 0, 1, \dots, n.\]

This reduces to a system of equations:

\[\sum_{k = 0}^n a_k\class{emph}{\phi_k(x_i)} = \class{emph}{y_i} \text{ for } i = 0, 1, \dots, n.\]

Matrix form:

\[ \class{emph}{\begin{bmatrix} \phi_0(x_0) & \phi_1(x_0) & \phi_2(x_0) & \cdots & \phi_n(x_0)\\ \phi_0(x_1) & \phi_1(x_1) & \phi_2(x_1) & \cdots & \phi_n(x_1)\\ \phi_0(x_2) & \phi_1(x_2) & \phi_2(x_2) & \cdots & \phi_n(x_2)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \phi_0(x_n) & \phi_1(x_n) & \phi_2(x_n) & \cdots & \phi_n(x_n)\\ \end{bmatrix}} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_n \end{bmatrix} = \class{emph}{\begin{bmatrix} y_0\\y_1\\y_2\\\vdots\\y_n \end{bmatrix}} \]

Polynomial Interpolation

Vector space: polynomials of degree at most \(n\)

Monomial basis: \(\{1, x, x^2, x^3, \dots, x^n\}\).

Van der Monde matrix:

\[ \begin{bmatrix} 1 & x_0^1 & x_0^2 & \cdots & x_0^n\\ 1 & x_1^1 & x_1^2 & \cdots & x_1^n\\ 1 & x_2^1 & x_2^2 & \cdots & x_2^n\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n^1 & x_n^2 & \cdots & x_n^n\\ \end{bmatrix} \]

Gauss-Jordan is slow, matrix is ill-conditioned.

Lagrange Basis

Idea: Find the basis function \(\ell_k\), \(k = 0, 1, \dots, n\) such that

Each \(\ell_k\) is a polynomial of degree \(n\).
For each \(k = 0, \dots, n\), for each \(i = 0, \dots, n\), \(\displaystyle\ell_k(x_i) = \delta_{ki} = \begin{cases} 1 & \text{ if } k = i\\0 & \text{ otherwise}\end{cases}\)

If we can find such polynomials, then:

The set \(\{\ell_0, \ell_1, \dots, \ell_n\}\) is linearly independent:
- \(\ell_k(x_k) = 1\) while \(\ell_j(x_k) = 0\) if \(j \neq k\).
- If \(p\) is any linear combination of polynomials \(\ell_j\) for \(j \neq k\), then \(p(x_k) = 0\).

The set \(\{\ell_0, \ell_1, \dots, \ell_n\}\) is a linearly independent set of \(n+1\) polynomials of degree at most \(n\).

This set forms a basis of the vector space of polynomials of degree at most \(n\).

Lagrange Basis Formulas

\[\ell_k(x) = \prod_{\substack{j = 0\\ j\neq k}}^n\frac{x - x_j}{x_k - x_j} = \frac{x - x_0}{x_k - x_0}\cdots \frac{x - x_{k-1}}{x_k - x_{k-1}}\cdot\frac{x-x_{k+1}}{x_k - x_{k+1}}\cdots\frac{x - x_n}{x_k-x_n}\]

\[\ell_k(x_k) = \prod_{\substack{j = 0\\ j\neq k}}^n\frac{x_k - x_j}{x_k - x_j} = 1\]

\[\ell_k(x_i) = \frac{x_i - x_0}{x_k - x_0}\cdots\class{standout}{\overbrace{\frac{x_i - x_i}{x_k - x_i}}^0}\cdots \frac{x - x_{k-1}}{x_k - x_{k-1}}\cdot\frac{x-x_{k+1}}{x_k - x_{k+1}}\cdots\frac{x - x_n}{x_k-x_n} = 0\]

\[\text{Then } \displaystyle p(x) = \sum_{k = 0}^n y_k\ell_k(x) = y_0\ell_0(x) + y_1\ell_1(x) + \cdots + y_n\ell_n(x)\]

Example 1

Find a polynomial \(p\) of degree 2 such that \(p(-1) = 11\), \(p(1) = -1\), and \(p(2) = 2\).

The “nodes”: \(x_0 = -1\), \(x_1 = 1\), \(x_2 = 2\).

\(\displaystyle\ell_0 = \frac{x - 1}{-1 -1}\frac{x - 2}{-1 - 2} = \frac{x - 1}{-2}\frac{x - 2}{-3} = \frac{1}{6}(x-1)(x-2)\)
\(\displaystyle\ell_1 = \frac{x - (-1)}{1 - (-1)}\frac{x - 2}{1 - 2} = \frac{x + 1}{2}\frac{x - 2}{-1} = -\frac{1}{2}(x+1)(x-2)\)
\(\displaystyle\ell_2 = \frac{x - (-1)}{2 - (-1)}\frac{x - 1}{2 - 1} = \frac{x + 1}{3}\frac{x - 1}{1} = \frac{1}{3}(x+1)(x-1)\)
Then \[p(x) = 11\ell_0(x) - \ell_1(x) + 2\ell_2(x) = \frac{11}{6}(x-1)(x-2) + \frac{1}{2}(x+1)(x-2) + \frac{2}{3}(x+1)(x-1)\]
Simplified \(\displaystyle p(x) = \frac{11}{6}x^2 - \frac{11}{2}x + \frac{11}{3} + \frac{1}{2}x^2 - \frac{1}{2}x - 1 + \frac{2}{3}x^2 - \frac{2}{3} = 3x^2 - 6x + 2\)

Lagrange Polynomials

What we have so far:

Monomial basis:
- Very hard to find interpolating polynomial
- The polynomial is in a very simple easy to use form
Lagrange basis:
- Very easy to find interpolating polynomial
- The polynomial is in a complicated hard to use form

A Compromise

The Newton basis:

\[ \begin{align} \pi_0(x) &= 1\\ \pi_k(x) &= \prod_{j=0}^{k-1}(x-x_j) \text{ for } k = 1, 2, \dots, n \end{align} \]

Degree of \(\pi_k\) is \(k\), for \(k = 0, 1, \dots, n\).
When adding polynomial of different degrees, the degree of the sum is always the highest of the degrees.
These polynomials are linerly independent.
There are \(n+1\) of them, with the highest degree being \(n\).
They form a basis of the space of polynomial of degree at most \(n\).
Important property: \(\pi_k(x_i) = 0\) for \(i < k\).

The System

\[ \class{emph}{\begin{bmatrix} \pi_0(x_0) & \pi_1(x_0) & \pi_2(x_0) & \cdots & \pi_{n-1}(x_0) & \pi_n(x_0)\\ \pi_0(x_1) & \pi_1(x_1) & \pi_2(x_1) & \cdots & \pi_{n-1}(x_1) & \pi_n(x_1)\\ \pi_0(x_2) & \pi_1(x_2) & \pi_2(x_2) & \cdots & \pi_{n-1}(x_2) & \pi_n(x_2)\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \pi_0(x_{n-1}) & \pi_1(x_{n-1}) & \pi_2(x_{n-1}) & \cdots & \pi_{n-1}(x_{n-1}) & \pi_n(x_{n-1})\\ \pi_0(x_n) & \pi_1(x_n) & \pi_2(x_n) & \cdots & \pi_{n-1}(x_n) & \pi_n(x_n)\\ \end{bmatrix}} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_{n-1}\\a_n \end{bmatrix} = \class{emph}{\begin{bmatrix} y_0\\y_1\\y_2\\\vdots\\y_{n-1}\\y_n \end{bmatrix}} \]

\[ \class{emph}{\begin{bmatrix} \pi_0(x_0) & 0 & 0 & \cdots & 0 & 0\\ \pi_0(x_1) & \pi_1(x_1) & 0 & \cdots & 0 & 0\\ \pi_0(x_2) & \pi_1(x_2) & \pi_2(x_2) & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \pi_0(x_{n-1}) & \pi_1(x_{n-1}) & \pi_2(x_{n-1}) & \cdots & \pi_{n-1}(x_{n-1}) & \mathrlap{\quad\;\,0}\phantom{\pi_n(x_{n-1}}{\;}\\ \pi_0(x_n) & \pi_1(x_n) & \pi_2(x_n) & \cdots & \pi_{n-1}(x_n) & \pi_n(x_n)\\ \end{bmatrix}} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_{n-1}\\a_n \end{bmatrix} = \class{emph}{\begin{bmatrix} y_0\\y_1\\y_2\\\vdots\\y_{n-1}\\y_n \end{bmatrix}} \]

\(\pi_k(x_i) = 0\) for \(i < k\).

\[ \class{emph}{\begin{bmatrix} 1 & 0 & 0 & \cdots & 0 & 0\\ 1 & (x_1-x_0) & 0 & \cdots & 0 & 0\\ 1 & (x_2-x_0) & (x_2-x_0)(x_2-x_1) & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & (x_{n-1}-x_0) & (x_{n-1}-x_0)(x_{n-1}-x_1) & \cdots & \prod_{j-0}^{n-2}(x_{n-1}-x_j) & 0\\ 1 & (x_n-x_0) & (x_n-x_0)(x_n-x_1) & \cdots & \prod_{j=0}^{n-2}(x_n - x_j) & \prod_{j=0}^{n-1}(x_n-x_j)\\ \end{bmatrix}} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_{n-1}\\a_n \end{bmatrix} = \class{emph}{\begin{bmatrix} y_0\\y_1\\y_2\\\vdots\\y_{n-1}\\y_n \end{bmatrix}} \]

Example 2

Use the Newton’s basis to find a polynomial \(p\) of degree 2 such that \(p(-1) = 11\), \(p(1) = -1\), and \(p(2) = 2\).

Nested Multiplication

\[ \begin{aligned} p(x) = a_0&\pi_0(x) + a_1\pi_1(x) + a_2\pi_2(x) + a_3\pi_3(x) + \cdots + a_{n-1}\pi_{n-1}(x) + a_n\pi_n(x)\\ = a_0&\\ {}+ a_1&{\color{green}(x-x_0)} \\ {}+ a_2&{\color{green}(x-x_0)}{\color{blue}(x-x_1)} \\ {}+ a_3&{\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)} \\ &\vdots \\ {}+ a_{n-1}&{\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})} \\ {}+ a_{n}&{\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}{\color{red}(x-x_{n-1})} \end{aligned} \]

\[ \begin{aligned} p(x) &= a_0\pi_0(x) + a_1\pi_1(x) + a_2\pi_2(x) + a_3\pi_3(x) + \cdots + a_{n-1}\pi_{n-1}(x) + a_n\pi_n(x)\\ &= a_0\\ &+ {\color{green}(x-x_0)}a_1 \\ &+ {\color{green}(x-x_0)}{\color{blue}(x-x_1)}a_2 \\ &+ {\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}a_3 \\ &\qquad\vdots \\ &+ {\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}a_{n-1} \\ &+ {\color{green}(x-x_0)}{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}{\color{red}(x-x_{n-1})}a_{n} \end{aligned} \]

\[ \begin{aligned} p(x) &= a_0\pi_0(x) + a_1\pi_1(x) + a_2\pi_2(x) + a_3\pi_3(x) + \cdots + a_{n-1}\pi_{n-1}(x) + a_n\pi_n(x)\\ &= a_0 + {\color{green}(x-x_0)} \begin{aligned}[t] \left(a_1\right. &+ {\color{blue}(x-x_1)}a_2 \\ &+ {\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}a_3 \\ &\qquad\vdots \\ &+ {\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}a_{n-1} \\ &+ \left.{\color{blue}(x-x_1)}{\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}{\color{red}(x-x_{n-1})}a_{n}\right) \end{aligned} \end{aligned} \]

\[ \begin{aligned} p(x) &= a_0\pi_0(x) + a_1\pi_1(x) + a_2\pi_2(x) + a_3\pi_3(x) + \cdots + a_{n-1}\pi_{n-1}(x) + a_n\pi_n(x)\\ &= a_0 + {\color{green}(x-x_0)} \begin{aligned}[t] {\color{green}\left(\strut\right.}a_1 &+ {\color{blue}(x-x_1)} \begin{aligned}[t] {\color{blue}\left(\strut\right.}a_2 &+ {\color{magenta}(x-x_2)}a_3 \\ &\qquad\vdots \\ &+ {\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}a_{n-1} \\ &+ {\color{magenta}(x-x_2)}\cdots{\color{orange}(x-x_{n-2})}{\color{red}(x-x_{n-1})}a_{n}{\color{blue}\left.\strut\right)}{\color{green}\left.\strut\right)} \end{aligned} \end{aligned} \end{aligned} \]

\[ \begin{aligned} p(x) &= a_0\pi_0(x) + a_1\pi_1(x) + a_2\pi_2(x) + a_3\pi_3(x) + \cdots + a_{n-1}\pi_{n-1}(x) + a_n\pi_n(x)\\ &= a_0 + {\color{green}(x-x_0)} \begin{aligned}[t] {\color{green}\left(\strut\right.}a_1 &+ {\color{blue}(x-x_1)} \begin{aligned}[t] {\color{blue}\left(\strut\right.}a_2 &+ {\color{magenta}(x-x_2)} \begin{aligned}[t] {\color{magenta}\left(\strut\right.}a_3 &+ \cdots {\color{magenta}\left.\strut\right)}{\color{blue}\left.\strut\right)}{\color{green}\left.\strut\right)} \end{aligned} \end{aligned} \end{aligned} \end{aligned} \]

\[ p(x) = a_0 + {\color{green}(x-x_0)\left(\strut\right.}a_1 + {\color{blue}(x-x_1)\left(\strut\right.}a_2 + {\color{magenta}(x-x_2)\left(\strut\right.}a_3 + \cdots + {\color{orange}(x-x_{n-2})\left(\strut\right.}a_{n-1} + {\color{red}(x-x_{n-1})\left(\strut\right.}a_{n} {\color{red}\left.\strut\right)} {\color{orange}\left.\strut\right)} \cdots {\color{magenta}\left.\strut\right)} {\color{blue}\left.\strut\right)} {\color{green}\left.\strut\right)} \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\left(\strut\right.} {\left(a_{n}\right)}\,(x-x_{n-1}) + \left.a_{n-1}\strut\right)(x-x_{n-2}) + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\left(\strut\right.} {\color{red}\left(a_{n}\right)}\,(x-x_{n-1}) + \left.a_{n-1}\strut\right)(x-x_{n-2}) + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\left(\strut\right.} {\color{red}\left(a_{n}\right)(x-x_{n-1})} + \left.a_{n-1}\strut\right)(x-x_{n-2}) + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\color{red}{\left(\strut\right.} {\left(a_{n}\right)}\,(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2}) + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\color{red}{\left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2})} + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} {\,\color{red}{\cdots \left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2}) + \cdots} + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\color{red}{{\left(\strut\right.} \cdots {\left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2})} + \cdots + \left.a_3\strut\right)}\,(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\left(\strut\right.} {\color{red}{{\left(\strut\right.} \cdots {\left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2})} + \cdots + \left.a_3\strut\right)(x-x_2)} + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\left(\strut\right.} {\color{red}{{\left(\strut\right.} {\left(\strut\right.} \cdots {\left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2})} + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)}\,(x-x_1) + \left.a_1\strut\right)(x-x_0) + a_0 \]

\[ p(x) = {\color{red}{{\left(\strut\right.} {\left(\strut\right.} {\left(\strut\right.} \cdots {\left(\strut\right.} {\left(a_{n}\right)(x-x_{n-1}) + \left.a_{n-1}\strut\right)}\,(x-x_{n-2})} + \cdots + \left.a_3\strut\right)(x-x_2) + \left.a_2\strut\right)(x-x_1) + \left.a_1\strut\right)}\,(x-x_0) + a_0 \]

Start with \(a_n\).

Multiply by \((x - x_{n-1})\).

Add \(a_{n-1}\).

Multiply by \((x - x_{n-2})\).

Add, then multiply, and so on…

Add \(a_3\).

Multiply by \((x - x_2)\).

Add \(a_2\).

Multiply by \((x - x_1)\).

Add \(a_1\).

Multiply by \((x - x_0)\).

Add \(a_0\).

Multiply, add, multiply, add, multiply, add, …

\(n\) multiplications, \(n\) additions.

This is also known as Horner’s algorithm.

Horner’s Algorithm

Evaluate the polynomial

\[ \begin{align} p(x_0, x_1, \dots, x_{n-1}) &= a_0 + a_1x_0 + a_2x_0x_1 + a_3x_0x_1x_2 + \cdots a_nx_0x_1x_2\cdots x_{n-1}\\ &= a_0 + x_0\left(a_1 + x_1\left(a_2 + x_2\left(a_3 + \cdots x_{n-1}\left(a_n\right)\right)\right)\right) \end{align} \]

Start with \(a_n\).
Multiply by \(x_{n-1}\).
Add \(a_{n-1}\).
\(\vdots\)
Multiply by \(x_k\).
Add \(a_k\).
\(\vdots\)
Multiply by \(x_0\).
Add \(a_0\).

Sequence of multiply-add steps.
A.k.a. synthetic division
Pro: very fast
Con: hard to parallelize

Example 3

Use the Horner’s algorithm to evaluate the polynomial

\[p(x) = 2 - 3(x+2) + 5(x+2)(x-1) - 2(x+2)(x-1)(x-3)\]

at \(x = 2\).