Polynomial interpolation

Suppose you have a list of \(n+1\) given point \((x_i, y_i)\) with \(i \in \{0, \dots, n\}\) and \(\forall i,j \in \{0, \dots, n\}: i \neq j \Rightarrow x_i \neq x_j\).

Now you want to find a polynomial \(\displaystyle p(x) = \sum_{i=0}^n a_i \cdot x^i\) that goes through all of the given points. This means:

\(\forall i \in \{0, \dots, n\}: p(x_i) = y_i\)

Existence and uniqueness

Theorem: When you have a list of \(n+1\) points with mutually different \(x_i\) there is exactly one polynom of degree \(\leq n\) that fits those points.

Proof:

You can formulate a linear system of equations:

$$\underbrace{\begin{pmatrix} x_0^0 & \cdots & x_0^n \\ \vdots & \ddots & \vdots \\ x_n^0 & \cdots & x_n^n \end{pmatrix} }_{:= A \in \mathbb{R}^{(n+1) \times (n+1)}} \begin{pmatrix} a_0 \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} y_0 \\ \vdots \\ y_n \end{pmatrix}$$

A matrix of the form of \(A\) is called Vandermonde matrix. When the data points \(x_i\) are mutually different, it is known that the Vandermonde matrix is invertible (source).

This means, the solution \((a_0 \dots a_n)^T\) of this linear equation gives the polynom \(p(x) = \sum_{i=0}^n a_i x^i\)

So the solution exists and is unique \(\blacksquare\)

Straight forward interpolating polynomials

For this algorithm, I'll find the polynomial in its monomial from \(p(x) = \sum_{i=0}^n a_i x^i\). I'll use the matrix \(A\) from section "Uniqueness".

You might want to take a look at my article about Gaussian elimination.

#!/usr/bin/env python
# -*- coding: utf-8 -*-


def pprintGaus(A):
    """ Pretty print a n×n matrix with a result vector n×1. """
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n + 1):
            line += str(A[i][j]) + "\t"
            if j == n - 1:
                line += "| "
        print(line)
    print("")


def pprintPolynomial(A):
    """ Pretty print a polynomial. """
    line = ""
    for i in range(len(x) - 1, -1, -1):
        if x[i] != 0:
            if i == 0:
                line += "+" + str(x[i])
            else:
                if x[i] == 1:
                    line += "+ x^" + str(i) + "\t"
                elif x[i] == -1:
                    line += "- x^" + str(i) + "\t"
                else:
                    line += "+" + str(x[i]) + "·x^" + str(i) + "\t"
    print(line)


def gauss(A):
    """ Solve a linear sysem of equations given by a n×n matrix
        with a result vector n×1. """
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i + 1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = A[k][i]
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n + 1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i + 1, n):
            c = -A[k][i] / A[i][i]
            for j in range(i, n + 1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n - 1, -1, -1):
        x[i] = A[i][n] / A[i][i]
        for k in range(i - 1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x


def setGauss(points):
    """ Create a system of equations for gaussian elimination from
        a set of points. """
    n = len(points) - 1
    A = [[0 for i in range(n + 2)] for j in range(n + 1)]
    for i in range(n + 1):
        x = points[i]["x"]
        for j in range(n + 1):
            A[i][j] = x ** j
        A[i][n + 1] = points[i]["y"]
    return A


if __name__ == "__main__":
    from fractions import Fraction

    # Read input data
    points = []
    points.append({"x": Fraction(-1), "y": Fraction(1)})
    points.append({"x": Fraction(1), "y": Fraction(1)})
    points.append({"x": Fraction(2), "y": Fraction(2)})

    A = setGauss(points)

    # Print input
    pprintGaus(A)

    # Calculate solution
    x = gauss(A)

    # Print result
    pprintPolynomial(x)

It is also interesting to get the value of \(p(x)\) at any given point \(x \in \mathbb{R}\):

def evaluatePolynomial(p, x):
    y = 0
    xi = 1
    for i, a in enumerate(p):
        y += a * xi
        xi *= x
    return y

Time complexity to get the polynomial: \(\frac{1}{3} n^3 + \mathcal{O}(n^2)\) (where \(n\) is the number of points) Space complexity for the polynomial: \(\mathcal{O}(n^2)\) (where \(n\) is the number of points)

Time complexity to evaluate the value of \(p(x)\) for any \(x \in \mathbb{R}\): \(\mathcal{O}(n)\).

Polynomials

\(\displaystyle \mathbb{R}_n[X] := \left \{p:\mathbb{R} \rightarrow \mathbb{R} | p(x) = \sum_{i=0}^n a_i \cdot x^i \text{ with } a_0, \dots, a_n \in \mathbb{R} \right \}\)

So \(\mathbb{R}_n[X]\) are all polynomials with real coefficients and degree not higher than latex]n\(. $\mathbb{R}_n[X]\) forms a vector space for \(n \in \mathbb{N}_0\). The degree of that vector space is \(n+1\).

What does that mean?

This means you can use different bases for polynomials. The base we usually use for \(\mathbb{R}_n[X]\) is \(\{x^0, x^1, x^2, x^3, \dots\}\), but you can switch the base.

Lagrange interpolation

Define \(n+1\) polynomials like this:

\(\displaystyle L_{i}(x) := \prod_{j=0 \atop j \neq i}^n \frac{x-x_j}{x_i - x_j}\)

This means:

\(\displaystyle L_{i}(x_i) := \prod_{j=0 \atop j \neq i}^n \frac{x_i-x_j}{x_i - x_j} = 1\)

and

\(\displaystyle L_{i}(x_p) := \prod_{j=0 \atop j \neq i}^n \frac{x_p-x_j}{x_i - x_j} = 0 \;\; p \in \{0, \dots, n\} \setminus \{i\}\)

So the polynomial \(y_i \cdot L_{i}(x)\) fits the point \((x_i, y_i)\) and is zero for all other points. This implies that \(\displaystyle p(x) = \sum_{i=0}^n y_i \cdot L_i(x)\) is an interpolation of our data points. The degree of \(p(x)\) is not higher than \(n\). You can see that easily when you take a look at the definition of \(p(x)\) and \(L_i(x)\).

The polynomials \(L_i(x)\) form another base for \(\mathbb{R}_n[X]\).

Lagranges way to interpolate polynomials can be implemented like this:

def lagrangeInterpolation(points):
    p = []
    for i in range(len(points)):
        Li = {"y": points[i]["y"], "polynomial": []}
        for j in range(len(points)):
            if j == i:
                continue
            Li["polynomial"].append(
                {"sub": points[j]["x"], "divisor": points[i]["x"] - points[j]["x"]}
            )
        p.append(Li)
    return p


def evaluateLagrangePolynomial(p, x):
    y = 0
    for Li in p:
        prod = 1
        for term in Li["polynomial"]:
            prod *= (x - term["sub"]) / term["divisor"]
        y += Li["y"] * prod
    return y

Time complexity to get the polynomial: \(n^2 + \mathcal{O}(n)\) (where \(n\) is the number of points) Space complexity for the polynomial: \(\mathcal{O}(n^2)\) (where \(n\) is the number of points)

Time complexity to evaluate the value of \(p(x)\) for any \(x \in \mathbb{R}\): \(\mathcal{O}(n^2)\).

Newton interpolation

Define

\(\displaystyle N_i(x) := \prod_{j=0}^{i-1} (x-x_j)\)

So you know that

\(N_i(x) = 0 \Leftrightarrow \exists p \in \{1, \dots, i\}: x = x_{i-p}\)

So your polynomial is

\(\displaystyle p(x) = \sum_{i=0}^n c_i \cdot N_i(x)\)

for the correct \(c_i\). You can do this by solving the following system of equations. Please note that you don't have to store that matrix to get those \(c_i\), divided differences are better.

$$\begin{pmatrix} 1 & & & & 0 \\ 1 & (x_1 - x_0) & & & \\ 1 & (x_2 - x_0) & (x_2 - x_0)(x_2 - x_1) & & \\ \vdots & \vdots & & \ddots & \\ 1 & (x_n - x_0) & \cdots & & \prod_{i=0}^{n-1}(x_n - x_i) \\ \end{pmatrix} \cdot \begin{pmatrix} c_0 \\ \vdots \\ c_n \end{pmatrix} = \begin{pmatrix} f_0 \\ \vdots \\ f_n \end{pmatrix}$$

You can get this lower triangular matrix like this:

def getGaussSystemForNewton(points):
    n = len(points) - 1
    A = [[0 for i in range(n + 2)] for j in range(n + 1)]
    for j in range(0, n + 2):
        for i in range(j, n + 1):
            if j == 0:
                A[i][j] = 1
            else:
                A[i][j] = A[i][j - 1] * (points[i]["x"] - points[j - 1]["x"])
        if j == n + 1:
            for i in range(0, n):
                A[i][j] = points[i]["y"]
    return A

From my previous posts about solving equations of upper triangular matrices and lower unitriangular matrices you know that the space complexity of this is in \(\Theta(n^2)\).

According to Wikipedia, you can use Horner's method to evaluate this Polynom in \(\mathcal{O}(n)\). But I really don't want to implement this today.

Interactive example

  • Click to add points.
  • Ctrl+Click to remove point.
  • You can enter a function that is valid JavaScript in the text field.

Published


by

Category

Code

Tags

Contact